Line data Source code
1 : ! ***********************************************************************
2 : !
3 : ! Copyright (C) 2020 The MESA Team
4 : !
5 : ! This program is free software: you can redistribute it and/or modify
6 : ! it under the terms of the GNU Lesser General Public License
7 : ! as published by the Free Software Foundation,
8 : ! either version 3 of the License, or (at your option) any later version.
9 : !
10 : ! This program is distributed in the hope that it will be useful,
11 : ! but WITHOUT ANY WARRANTY; without even the implied warranty of
12 : ! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
13 : ! See the GNU Lesser General Public License for more details.
14 : !
15 : ! You should have received a copy of the GNU Lesser General Public License
16 : ! along with this program. If not, see <https://www.gnu.org/licenses/>.
17 : !
18 : ! ***********************************************************************
19 :
20 : module create_FSCRliq8_table
21 :
22 : use const_def, only: dp
23 : use chem_def
24 : use utils_lib, only: is_bad
25 : use math_lib
26 :
27 : implicit none
28 :
29 : public :: do_create_FSCRliq8_table
30 : private
31 :
32 :
33 : contains
34 :
35 :
36 0 : subroutine do_create_FSCRliq8_table(fname,iZion)
37 : character (len=*), intent(in) :: fname
38 : integer, intent(in) :: iZion
39 : real(dp) :: logRS, logRS_min, logRS_max, dlogRS
40 0 : real(dp) :: logGAME, logGAME_min, logGAME_max, dlogGAME
41 : integer :: nlogRS, nlogGAME, io_unit, i, j, ierr
42 0 : real(dp) :: Zion, RS, GAME, FSCR, USCR, PSCR, CVSCR, PDTSCR, PDRSCR
43 : include 'formats'
44 :
45 0 : Zion = iZion
46 :
47 0 : logRS_min = -3.5d0
48 0 : logRS_max = 0.0d0
49 0 : dlogRS = 1d-2
50 0 : nlogRS = (logRS_max - logRS_min)/dlogRS + 1
51 :
52 0 : logGAME_min = -2d0
53 0 : logGAME_max = 4d0
54 0 : dlogGAME = 1d-2
55 0 : nlogGAME = (logGAME_max - logGAME_min)/dlogGAME + 1
56 :
57 0 : io_unit = 40
58 :
59 : !write(*,'(a)') 'create ' // trim(fname)
60 :
61 0 : open(unit=io_unit,file=trim(fname))
62 :
63 0 : write(io_unit,'(99(a14))') 'num logRS', 'logRS min', 'logRS max', 'del logRS', &
64 0 : 'num logGAME', 'logGAME min', 'logGAME max', 'del logGAME', 'Zion'
65 :
66 : write(io_unit,'(2(i10,4x,3(f14.4)),i14)') &
67 0 : nlogRS, logRS_min, logRS_max, dlogRS, &
68 0 : nlogGAME, logGAME_min, logGAME_max, dlogGAME, iZion
69 :
70 0 : do i = 1, nlogRS
71 0 : logRS = logRS_min + (i-1) * dlogRS
72 0 : RS = exp10(logRS)
73 0 : write(io_unit,'(/,7x,a)') 'logRS'
74 0 : write(io_unit,'(2x,f14.6/)') logRS
75 : write(io_unit,'(99(a26))') &
76 0 : 'logGAME', 'FSCR', 'USCR', 'PSCR', 'CVSCR', 'PDTSCR', 'PDRSCR', 'RS', 'GAME'
77 0 : do j = 1, nlogGAME
78 0 : logGAME = logGAME_min + (j-1) * dlogGAME
79 0 : GAME = exp10(logGAME)
80 0 : if (is_bad(GAME)) then
81 0 : write(*,3) 'GAME', i, j, GAME, logGAME, logGAME_min, dlogGAME
82 : end if
83 : ierr = 0
84 0 : call FSCRliq8(RS, GAME, Zion, FSCR, USCR, PSCR, CVSCR, PDTSCR, PDRSCR, ierr)
85 0 : if (ierr /= 0) then
86 0 : write(*,*) 'FSCRliq8 error'
87 0 : write(*,1) 'RS', RS
88 0 : write(*,1) 'GAME', GAME
89 0 : write(*,1) 'Zion', Zion
90 0 : call mesa_error(__FILE__,__LINE__,'do_create_FSCRliq8_table')
91 : end if
92 : write(io_unit,'(99(1pd26.16))') &
93 0 : logGAME, FSCR, USCR, PSCR, CVSCR, PDTSCR, PDRSCR, RS, GAME
94 : end do
95 0 : write(io_unit,*)
96 : end do
97 :
98 0 : write(io_unit,*)
99 0 : write(io_unit,*)
100 0 : close(io_unit)
101 :
102 0 : end subroutine do_create_FSCRliq8_table
103 :
104 :
105 0 : subroutine FSCRliq8(RS,GAME,Zion, &
106 : FSCR,USCR,PSCR,CVSCR,PDTSCR,PDRSCR,ierr) ! fit to the el.-ion scr.
107 : ! Version 11.09.08
108 : ! cleaned 16.06.09
109 : ! Stems from FSCRliq7 v. 09.06.07. Included a check for RS=0.
110 : ! INPUT: RS - density parameter, GAME - electron Coulomb parameter,
111 : ! Zion - ion charge number,
112 : ! OUTPUT: FSCR - screening (e-i) free energy per kT per 1 ion,
113 : ! USCR - internal energy per kT per 1 ion (screen.contrib.)
114 : ! PSCR - pressure divided by (n_i kT) (screen.contrib.)
115 : ! CVSCR - heat capacity per 1 ion (screen.contrib.)
116 : ! PDTSCR,PDRSCR = PSCR + d PSCR / d ln(T,\rho)
117 : real(dp), intent(in) :: RS,GAME
118 : real(dp), intent(in) :: Zion
119 : real(dp), intent(out) :: FSCR,USCR,PSCR,CVSCR,PDTSCR,PDRSCR
120 : integer, intent(out) :: ierr
121 :
122 0 : real(dp) :: SQG, SQR, SQZ1, SQZ, CDH0, CDH, ZLN, Z13
123 0 : real(dp) :: X, CTF, P01, P03, PTX
124 0 : real(dp) :: TX, TXDG, TXDGG, TY1, TY1DG, TY1DGG
125 0 : real(dp) :: TY2, TY2DX, TY2DXX, TY, TYX, TYDX, TYDG, P1
126 0 : real(dp) :: COR1, COR1DX, COR1DG, COR1DXX, COR1DGG, COR1DXG
127 0 : real(dp) :: U0, U0DX, U0DG, U0DXX, U0DGG, U0DXG
128 0 : real(dp) :: D0DG, D0, D0DX, D0DXX
129 0 : real(dp) :: COR0, COR0DX, COR0DG, COR0DXX, COR0DGG, COR0DXG
130 0 : real(dp) :: RELE, Q1, Q2, H1U, H1D, H1, H1X, H1DX, H1DXX
131 0 : real(dp) :: UP, UPDX, UPDG, UPDXX, UPDGG, UPDXG
132 0 : real(dp) :: DN1, DN1DX, DN1DG, DN1DXX, DN1DGG, DN1DXG
133 0 : real(dp) :: DN, DNDX, DNDG, DNDXX, DNDGG, DNDXG
134 0 : real(dp) :: FX, FXDG, FDX, FG, FDG, FDGDH, FDXX, FDGG, FDXG
135 :
136 : real(dp), parameter :: XRS=.0140047d0
137 : real(dp), parameter :: TINY=1.d-19
138 :
139 0 : ierr = 0
140 0 : if (RS<0d0) then
141 0 : ierr = -1
142 0 : return
143 : !call mesa_error(__FILE__,__LINE__,'FSCRliq8: RS < 0')
144 : end if
145 0 : if (RS<TINY) then
146 0 : FSCR=0.d0
147 0 : USCR=0.d0
148 0 : PSCR=0.d0
149 0 : CVSCR=0.d0
150 0 : PDTSCR=0.d0
151 0 : PDRSCR=0.d0
152 0 : return
153 : end if
154 0 : SQG=sqrt(GAME)
155 0 : SQR=sqrt(RS)
156 0 : SQZ1=sqrt(1d0+Zion)
157 0 : SQZ=sqrt(Zion)
158 0 : CDH0=Zion/1.73205d0 ! 1.73205=sqrt(3.)
159 0 : CDH=CDH0*(SQZ1*SQZ1*SQZ1-SQZ*SQZ*SQZ-1d0)
160 0 : SQG=sqrt(GAME)
161 0 : ZLN=log(Zion)
162 0 : Z13=exp(ZLN/3.d0) ! Zion**(1./3.)
163 0 : X=XRS/RS ! relativity parameter
164 0 : CTF=Zion*Zion*.2513d0*(Z13-1d0+.2d0/sqrt(Z13))
165 : ! Thomas-Fermi constant; .2513=(18/175)(12/\pi)^{2/3}
166 0 : P01=1.11d0*exp(0.475d0*ZLN)
167 0 : P03=0.2d0+0.078d0*ZLN*ZLN
168 0 : PTX=1.16d0+0.08d0*ZLN
169 0 : TX=pow(GAME,PTX)
170 0 : TXDG=PTX*TX/GAME
171 0 : TXDGG=(PTX-1.d0)*TXDG/GAME
172 0 : TY1=1d0/(1.d-3*Zion*Zion+2d0*GAME)
173 0 : TY1DG=-2d0*TY1*TY1
174 0 : TY1DGG=-4d0*TY1*TY1DG
175 0 : TY2=1d0+6d0*RS*RS
176 0 : TY2DX=-12d0*RS*RS/X
177 0 : TY2DXX=-3d0*TY2DX/X
178 0 : TY=RS*RS*RS/TY2*(1d0+TY1)
179 0 : TYX=3d0/X+TY2DX/TY2
180 0 : TYDX=-TY*TYX
181 0 : TYDG=RS*RS*RS*TY1DG/TY2
182 0 : P1=(Zion-1d0)/9d0
183 0 : COR1=1d0+P1*TY
184 0 : COR1DX=P1*TYDX
185 0 : COR1DG=P1*TYDG
186 0 : COR1DXX=P1*(TY*(3d0/(X*X)+(TY2DX/TY2)*(TY2DX/TY2)-TY2DXX/TY2)-TYDX*TYX)
187 0 : COR1DGG=P1*RS*RS*RS*TY1DGG/TY2
188 0 : COR1DXG=-P1*TYDG*TYX
189 0 : U0=0.78d0*sqrt(GAME/Zion)*RS*RS*RS
190 0 : U0DX=-3d0*U0/X
191 0 : U0DG=.5d0*U0/GAME
192 0 : U0DXX=-4.d0*U0DX/X
193 0 : U0DGG=-.5d0*U0DG/GAME
194 0 : U0DXG=-3.d0*U0DG/X
195 0 : D0DG=Zion*Zion*Zion
196 0 : D0=GAME*D0DG+21d0*RS*RS*RS
197 0 : D0DX=-63d0*RS*RS*RS/X
198 0 : D0DXX=252d0*RS*RS*RS/(X*X)
199 0 : COR0=1d0+U0/D0
200 0 : COR0DX=(U0DX-U0*D0DX/D0)/D0
201 0 : COR0DG=(U0DG-U0*D0DG/D0)/D0
202 0 : COR0DXX=(U0DXX-(2d0*U0DX*D0DX+U0*D0DXX)/D0+2d0*(D0DX/D0)*(D0DX/D0))/D0
203 0 : COR0DGG=(U0DGG-2d0*U0DG*D0DG/D0+2d0*U0*(D0DG/D0)*(D0DG/D0))/D0
204 0 : COR0DXG=(U0DXG-(U0DX*D0DG+U0DG*D0DX)/D0+2d0*U0*D0DX*D0DG/(D0*D0))/D0
205 : ! Relativism:
206 0 : RELE=sqrt(1.d0+X*X)
207 0 : Q1=0.18d0/sqrt(sqrt(Zion))
208 0 : Q2=0.2d0+0.37d0/sqrt(Zion)
209 0 : H1U=1d0+X*X/5.d0
210 0 : H1D=1d0+Q1*X+Q2*X*X
211 0 : H1=H1U/H1D
212 0 : H1X=0.4d0*X/H1U-(Q1+2d0*Q2*X)/H1D
213 0 : H1DX=H1*H1X
214 : H1DXX=H1DX*H1X &
215 0 : + H1*(0.4d0/H1U-(0.4d0*X/H1U)*(0.4d0*X/H1U)-2d0*Q2/H1D+pow2((Q1+2d0*Q2*X)/H1D))
216 0 : UP=CDH*SQG+P01*CTF*TX*COR0*H1
217 0 : UPDX=P01*CTF*TX*(COR0DX*H1+COR0*H1DX)
218 0 : UPDG=.5d0*CDH/SQG+P01*CTF*(TXDG*COR0+TX*COR0DG)*H1
219 0 : UPDXX=P01*CTF*TX*(COR0DXX*H1+2d0*COR0DX*H1DX+COR0*H1DXX)
220 : UPDGG=-.25d0*CDH/(SQG*GAME) &
221 0 : + P01*CTF*(TXDGG*COR0+2d0*TXDG*COR0DG+TX*COR0DGG)*H1
222 : UPDXG=P01*CTF*(TXDG*(COR0DX*H1+COR0*H1DX) &
223 0 : + TX*(COR0DXG*H1+COR0DG*H1DX))
224 0 : DN1=P03*SQG+P01/RS*TX*COR1
225 0 : DN1DX=P01*TX*(COR1/XRS+COR1DX/RS)
226 0 : DN1DG=.5d0*P03/SQG+P01/RS*(TXDG*COR1+TX*COR1DG)
227 0 : DN1DXX=P01*TX/XRS*(2d0*COR1DX+X*COR1DXX)
228 : DN1DGG=-.25d0*P03/(GAME*SQG) &
229 0 : + P01/RS*(TXDGG*COR1+2d0*TXDG*COR1DG+TX*COR1DGG)
230 0 : DN1DXG=P01*(TXDG*(COR1/XRS+COR1DX/RS)+TX*(COR1DG/XRS+COR1DXG/RS))
231 0 : DN=1d0+DN1/RELE
232 0 : DNDX=DN1DX/RELE-X*DN1/(RELE*RELE*RELE)
233 0 : DNDXX=(DN1DXX-((2d0*X*DN1DX+DN1)-3.d0*X*X*DN1/(RELE*RELE))/(RELE*RELE))/RELE
234 0 : DNDG=DN1DG/RELE
235 0 : DNDGG=DN1DGG/RELE
236 0 : DNDXG=DN1DXG/RELE-X*DN1DG/(RELE*RELE*RELE)
237 0 : FSCR=-UP/DN*GAME
238 0 : FX=(UP*DNDX/DN-UPDX)/DN
239 0 : FXDG=((UPDG*DNDX+UPDX*DNDG+UP*DNDXG-2d0*UP*DNDX*DNDG/DN)/DN-UPDXG)/DN
240 0 : FDX=FX*GAME
241 0 : FG=(UP*DNDG/DN-UPDG)/DN
242 0 : FDG=FG*GAME-UP/DN
243 0 : FDGDH=SQG*DNDG/(DN*DN) ! d FDG / d CDH
244 0 : FDXX=((UP*DNDXX+2d0*(UPDX*DNDX-UP*DNDX*DNDX/DN))/DN-UPDXX)/DN*GAME
245 0 : FDGG=2d0*FG+GAME*((2d0*DNDG*(UPDG-UP*DNDG/DN)+UP*DNDGG)/DN-UPDGG)/DN
246 0 : FDXG=FX+GAME*FXDG
247 0 : USCR=GAME*FDG
248 0 : CVSCR=-GAME*GAME*FDGG
249 0 : PSCR=(X*FDX+GAME*FDG)/3.d0
250 0 : PDTSCR=-GAME*GAME*(X*FXDG+FDGG)/3.d0
251 0 : PDRSCR=(12d0*PSCR+X*X*FDXX+2d0*X*GAME*FDXG+GAME*GAME*FDGG)/9d0
252 0 : return
253 : end subroutine FSCRliq8
254 :
255 :
256 : end module create_FSCRliq8_table
|
