Line data Source code
1 : ! ***********************************************************************
2 : !
3 : ! Copyright (C) 2013-2021 Josiah Schwab & 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 eval_coulomb
21 :
22 : use const_def, only: dp
23 : use math_lib
24 : use auto_diff
25 :
26 : implicit none
27 :
28 : integer, parameter :: None = 0
29 :
30 : ! select the option for the ion chemical potential
31 : integer, parameter :: DGC1973 = 1
32 : integer, parameter :: I1993 = 2
33 : integer, parameter :: PCR2009 = 3
34 :
35 : ! select the option for the electron screening
36 : integer, parameter :: ThomasFermi = 1
37 : integer, parameter :: Itoh2002 = 2
38 :
39 : contains
40 :
41 26 : function do_mui_coulomb(cc, Z) result(mu)
42 :
43 : use rates_def, only: Coulomb_Info, which_mui_coulomb
44 :
45 :
46 : type(Coulomb_Info), intent(in) :: cc
47 : real(dp), intent(in) :: Z ! nuclear charge
48 :
49 : type(auto_diff_real_2var_order1) :: mu ! the chemical potential in units of kT
50 :
51 8 : select case (which_mui_coulomb)
52 : case (None)
53 8 : mu = 0
54 : case (DGC1973)
55 0 : mu = mui_coulomb_DGC1973(cc, Z)
56 : case (I1993)
57 0 : mu = mui_coulomb_I1993(cc, Z)
58 : case (PCR2009)
59 12 : mu = mui_coulomb_PCR2009(cc, Z)
60 : case default
61 20 : stop "STOP (eval_coulomb.f: do_mui_coulomb): incorrect value of which_mui_coulomb"
62 : end select
63 :
64 20 : return
65 :
66 20 : end function do_mui_coulomb
67 :
68 :
69 0 : function mui_coulomb_DGC1973(cc, Z) result(mu)
70 20 : use math_lib
71 : use rates_def, only: Coulomb_Info
72 :
73 :
74 : type(Coulomb_Info), intent(in) :: cc
75 : real(dp), intent(in) :: Z ! nuclear charge
76 :
77 : type(auto_diff_real_2var_order1) :: mu ! the chemical potential in units of kT
78 :
79 : ! calculate the chemical potential of an ion
80 :
81 : type(auto_diff_real_2var_order1) :: gamma
82 0 : real(dp) :: zr, zr_m1o3
83 :
84 : ! define parameters
85 : real(dp), parameter :: c0 = 0.9d0
86 : real(dp), parameter :: c1 = 0.2843d0
87 : real(dp), parameter :: c2 = -0.054d0
88 : real(dp), parameter :: d0 = -9d0/ 16d0
89 : real(dp), parameter :: d1 = 0.460d0
90 :
91 : ! from Dewitt, H. E., Graboske, H. C., & Cooper, M. S. 1973, ApJ, 181, 439
92 : ! via Couch & Loumos ApJ, vol. 194, Dec. 1, 1974, pt. 1, p. 385-392
93 :
94 : ! coulomb coupling parameter
95 0 : gamma = pow(Z,2d0/3d0) * cc% zbar * cc% gamma_e
96 :
97 : ! it appears that Gutierrez et al. use something like
98 : ! despite citing the above references
99 : ! gamma = pow(cc% zbar,5d0/3d0) * cc% gamma_e
100 :
101 : ! ratios for convenience
102 0 : zr = Z/cc% zbar
103 0 : zr_m1o3 = pow(zr,-1d0/3d0) ! zr to the minus 1/3 power
104 :
105 : ! evaluate mu(Z); C&L eq. (11)
106 : mu = -zr * (gamma * (c0 + (c1 + c2 * zr_m1o3) * zr_m1o3) + &
107 0 : (d0 + d1 * zr_m1o3))
108 :
109 0 : return
110 :
111 0 : end function mui_coulomb_DGC1973
112 :
113 :
114 0 : function mui_coulomb_I1993(cc, Z) result(mu)
115 0 : use math_lib
116 : use rates_def, only: Coulomb_Info
117 :
118 :
119 : type(Coulomb_Info), intent(in) :: cc
120 : real(dp), intent(in) :: Z ! nuclear charge
121 :
122 : type(auto_diff_real_2var_order1) :: mu ! the chemical potential in units of kT
123 :
124 : ! calculate the chemical potential of an ion
125 :
126 : ! form from W.L. Slattery, G.D. Doolen, H.E. DeWitt, Phys. Rev. A 26 (1982) 2255.
127 : ! values from Ichimaru, S. 1993, Reviews of Modern Physics, 65, 255
128 : ! via Juodagalvis et al. (2010)
129 :
130 : type(auto_diff_real_2var_order1) :: gamma, gamma14
131 :
132 : ! define parameters
133 : real(dp), parameter :: a = -0.898004d0
134 : real(dp), parameter :: b = 0.96786d0
135 : real(dp), parameter :: c = 0.220703d0
136 : real(dp), parameter :: d = -0.86097d0
137 : real(dp), parameter :: e = -2.52692d0
138 :
139 : ! coulomb coupling parameter for species with charge Z
140 0 : gamma = cc% gamma_e * pow(Z,5d0/3d0)
141 :
142 : ! ratios for convenience
143 0 : gamma14 = pow(gamma, 1d0/4d0)
144 :
145 0 : mu = a * gamma + 4d0 * b * gamma14 - 4d0 * c / gamma14 + d * log(gamma) + e
146 :
147 0 : return
148 :
149 0 : end function mui_coulomb_I1993
150 :
151 :
152 12 : function mui_coulomb_PCR2009(cc, Z) result(mu)
153 0 : use math_lib
154 : use rates_def, only: Coulomb_Info
155 :
156 : ! calculate the chemical potential of an ion
157 :
158 :
159 : type(Coulomb_Info), intent(in) :: cc
160 : real(dp), intent(in) :: Z ! nuclear charge
161 :
162 : type(auto_diff_real_2var_order1) :: mu ! the chemical potential in units of kT
163 :
164 : type(auto_diff_real_2var_order1) :: gamma
165 :
166 : ! expressions taken from
167 : ! [CP98]: Chabrier, G., \& Potekhin, A.~Y.\ 1998, \pre, 58, 4941
168 : ! [PC00]: Potekhin, A.~Y., \& Chabrier, G.\ 2000, \pre, 62, 8554
169 : ! [P+09a]: Potekhin, A.~Y., Chabrier, G., \& Rogers, F.~J.\ 2009, \pre, 79, 016411
170 : ! [P+09b]: Potekhin, A.~Y., Chabrier, G., Chugunov, A.~I., Dewitt, H.~E., \& Rogers, F.~J.\ 2009, \pre, 80, 047401
171 :
172 : ! coulomb coupling parameter for species with charge Z
173 12 : gamma = cc% gamma_e * pow(Z, 5d0/3d0)
174 :
175 : ! first, calculate the shift from the linear mixing relation (LMR)
176 : ! see equation 2 from [P+09b]
177 :
178 12 : mu = fii(gamma) + fie(cc% rs, cc% gamma_e, Z)
179 :
180 : ! skip LMR corrections
181 :
182 12 : return
183 :
184 : contains
185 :
186 12 : function fii(gamma) result(f)
187 :
188 :
189 : type(auto_diff_real_2var_order1), intent(in) :: gamma
190 : type(auto_diff_real_2var_order1) :: f
191 :
192 : ! define parameters
193 : real(dp), parameter :: A1 = -0.907347d0
194 : real(dp), parameter :: A2 = 0.62849d0
195 : real(dp) :: A3
196 : real(dp), parameter :: B1 = 0.0045d0
197 : real(dp), parameter :: B2 = 170d0
198 : real(dp), parameter :: B3 = -8.4d-5
199 : real(dp), parameter :: B4 = 0.0037d0
200 :
201 12 : A3 = -sqrt(3d0)/2d0 - A1/sqrt(A2)
202 :
203 : ! [CP00] Eq. (28)
204 : f = A1 * (sqrt(gamma * (A2 + gamma)) - A2 * log(sqrt(gamma/A2) + sqrt(1+gamma/A2))) + &
205 12 : 2 * A3 * (sqrt(gamma) - atan(sqrt(gamma)))
206 :
207 12 : return
208 :
209 12 : end function fii
210 :
211 12 : function fie(rs, gamma_e, Z) result(f)
212 :
213 :
214 : type(auto_diff_real_2var_order1), intent(in) :: rs
215 : type(auto_diff_real_2var_order1), intent(in) :: gamma_e
216 : real(dp), intent(in) :: Z
217 :
218 : type(auto_diff_real_2var_order1) :: f
219 :
220 12 : real(dp) :: logZ
221 12 : real(dp) :: c_DH, c_inf, c_TF
222 : real(dp) :: a, b, nu
223 :
224 : type(auto_diff_real_2var_order1) :: x, g1, g2, h1, h2
225 :
226 : ! [CP00], between eq 31 and 32
227 12 : logZ = log(Z)
228 12 : a = 1.11_dp * pow(Z, 0.475_dp)
229 12 : b = 0.2_dp + 0.078_dp * logZ * logZ
230 12 : nu = 1.16_dp + 0.08_dp * logZ
231 :
232 : ! [CP00], eq 30 and eq 31
233 12 : c_DH = Z / sqrt(3d0) * (pow(1+Z,3d0/2d0) - 1d0 - pow(Z,3d0/2d0))
234 12 : c_inf = 0.2513_dp
235 12 : c_TF = c_inf * pow(Z,7d0/3d0) * (1d0 - pow(Z,-1d0/3d0) + 0.2d0 * pow(Z,-1d0/2d0))
236 :
237 : ! [CP00], eq 32 and eq 33
238 12 : g1 = 1 + (0.78d0 / (21d0 + gamma_e * pow(Z/rs, 3d0))) * sqrt(gamma_e/Z)
239 12 : g2 = 1 + (Z-1d0)/9d0 * (1d0 + 1d0/(0.001d0*Z*Z + 2*gamma_e)) * pow(rs,3d0) / (1d0 + 6d0 * rs * rs)
240 :
241 : ! [CP00] eq 4
242 12 : x = 0.014005_dp / rs
243 :
244 : ! [CP98], below eq 33
245 12 : h1 = 1d0 / (1d0 + pow(x/sqrt(1+x*x),6d0) * pow(Z,-1d0/3d0))
246 12 : h2 = 1d0 / sqrt(1+x*x)
247 :
248 : f = - gamma_e * (c_DH * sqrt(gamma_e) + c_TF * a * pow(gamma_e, nu) * g1 * h1) / &
249 12 : (1d0 + (b * sqrt(gamma_e) + a * g2 * pow(gamma_e, nu)/rs) * h2)
250 :
251 :
252 12 : end function fie
253 :
254 : end function mui_coulomb_PCR2009
255 :
256 :
257 :
258 10 : function do_Vs_coulomb(cc, Z) result(Vs)
259 :
260 : use rates_def, only: Coulomb_Info, which_Vs_coulomb
261 :
262 :
263 : type(Coulomb_Info), intent(in) :: cc
264 : real(dp), intent(in) :: Z ! nuclear charge
265 :
266 : type(auto_diff_real_2var_order1) :: Vs ! the screening potential in units of the fermi energy
267 :
268 4 : select case (which_Vs_coulomb)
269 : case (None)
270 4 : Vs = 0
271 : case (ThomasFermi)
272 0 : Vs = Vs_coulomb_TF(cc, Z)
273 : case (Itoh2002)
274 6 : Vs = Vs_coulomb_Itoh(cc, Z)
275 : case default
276 10 : stop "STOP (eval_coulomb.f: do_Vs_coulomb): incorrect value of which_Vs_coulomb"
277 : end select
278 :
279 10 : return
280 :
281 10 : end function do_Vs_coulomb
282 :
283 :
284 0 : function Vs_coulomb_TF(cc, Z) result(Vs)
285 10 : use const_def, only : fine, pi
286 : use math_lib
287 : use rates_def, only: Coulomb_Info
288 :
289 :
290 : type(Coulomb_Info), intent(in) :: cc
291 : real(dp), intent(in) :: Z ! nuclear charge
292 :
293 : type(auto_diff_real_2var_order1) :: Vs ! the screening potential in units of the fermi energy
294 :
295 : ! this uses the thomas-fermi approximation, in the relativistic limit
296 0 : Vs = 2d0 * Z * fine * sqrt(fine/pi)
297 :
298 0 : return
299 :
300 0 : end function Vs_coulomb_TF
301 :
302 :
303 6 : function Vs_coulomb_Itoh(cc, Z) result(Vs)
304 0 : use math_lib
305 : use rates_def, only: Coulomb_Info
306 :
307 :
308 : type(Coulomb_Info), intent(in) :: cc
309 : real(dp), intent(in) :: Z ! nuclear charge
310 :
311 : type(auto_diff_real_2var_order1) :: Vs ! the screening potential in units of the fermi energy
312 :
313 : ! code from Itoh, N., Tomizawa, N., Tamamura, M., Wanajo, S., & Nozawa, S. 2002, ApJ, 579, 380
314 :
315 : integer :: i
316 : type(auto_diff_real_2var_order1) :: rs, rs0, s, fj
317 : real(dp), dimension(11), parameter :: c(0:10) = [ &
318 : 0.450861D-01, 0.113078D-02, 0.312104D-02, 0.864302D-03, &
319 : 0.157214D-01, 0.816962D-01, 0.784921D-01,-0.680863D-01, &
320 : -0.979967D-01, 0.204907D-01, 0.366713D-01 ]
321 :
322 6 : rs = cc% rs
323 :
324 6 : rs0=0
325 6 : if(rs<0.00001d0)then
326 0 : rs0=1d0-rs
327 0 : rs=0.00001d0
328 : end if
329 6 : s=(log10(rs)+3d0)/2d0
330 6 : fj=0
331 6 : if(s/=0)then
332 72 : do i=0,10
333 72 : fj=c(i)*pow(s,dble(i))+fj
334 : end do
335 : else
336 0 : fj=c(0)
337 : end if
338 6 : if(rs0/=0)then
339 0 : rs=1d0-rs0
340 : end if
341 :
342 6 : Vs = 1.46d-2 * Z * fj
343 :
344 6 : return
345 :
346 6 : end function Vs_coulomb_Itoh
347 :
348 :
349 :
350 : end module eval_coulomb
|
