Line data Source code
1 : ! ***********************************************************************
2 : !
3 : ! Copyright (C) 2018 Sam Jones, Robert Farmer & 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 : ! Implement screening a la Chugunov, DeWitt & Yakovlev 2007, PhRvD, 76, 025028
21 :
22 : module screening_chugunov
23 : use math_lib
24 : use const_def, only: dp, amu, hbar, kerg, four_thirds, pi, pi2, pi4, qe, rbohr
25 : use rates_def, only: screen_info
26 : use math_def
27 :
28 : implicit none
29 : private
30 :
31 : ! there are various realms of applicability in the Chugunov paper, for
32 : ! the particle-in-cell Monte-Carlo calculations, for the MKB
33 : ! approximation, and for the Potekhin-Cabrier-based fit to the MKB
34 : ! results for the mean field potential H(r).
35 : ! hence, here are parameters for smoothly limiting the domain of application
36 : real(dp), parameter :: gamfitlim = 600.0d0 !< upper limit for applicability of fit for h(gamma)
37 : real(dp), parameter :: g0 = 590.0d0 !< lower limit to start blending from
38 : real(dp), parameter :: deltagam = gamfitlim - g0 ! How much to blend over the gamma boundaries
39 : real(dp), parameter :: tp2 = 0.2d0 !< minimum allowed ratio T/T_p before fading out
40 : real(dp), parameter :: tp1 = 0.1d0 !< floor value of T/T_p (T_p = plasma temperature)
41 : real(dp), parameter :: deltatp = tp2 - tp1 !< Blend over the tp boundary
42 :
43 : ! coefficients from chugunov
44 : real(dp), parameter :: c_a1 = 2.7822d0
45 : real(dp), parameter :: c_a2 = 98.34d0
46 : real(dp), parameter :: c_a3 = sqrt(3.0d0) - c_a1/sqrt(c_a2)
47 : real(dp), parameter :: c_b1 = -1.7476d0
48 : real(dp), parameter :: c_b2 = 66.07d0
49 : real(dp), parameter :: c_b3 = 1.12d0
50 : real(dp), parameter :: c_b4 = 65.d0
51 : real(dp), parameter :: alfa = 0.022d0
52 :
53 : real(dp), parameter :: x13 = 1.0d0/3.0d0
54 : real(dp), parameter :: x12 = 1.0d0/2.0d0
55 : real(dp), parameter :: x23 = 2.0d0/3.0d0
56 : real(dp), parameter :: x32 = 3.0d0/2.0d0
57 : real(dp), parameter :: x43 = 4.0d0/3.0d0
58 : real(dp), parameter :: x53 = 5.0d0/3.0d0
59 :
60 : real(dp), parameter :: h0fitlim = 300d0, h0fit0 = 295d0
61 : real(dp), parameter :: deltah0fit = h0fitlim - h0fit0
62 :
63 : public eval_screen_chugunov
64 :
65 : contains
66 :
67 :
68 2079846 : subroutine eval_screen_chugunov(sc, z1, z2, a1, a2, screen, dscreendt, dscreendd, ierr)
69 :
70 : type (Screen_Info) :: sc
71 : real(dp),intent(in) :: z1, z2 !< charge numbers of reactants
72 : real(dp),intent(in) :: a1, a2 !< mass numbers of reactants
73 : real(dp),intent(out) :: screen !< on return, screening factor for this reaction
74 : real(dp),intent(out) :: dscreendt !< on return, temperature derivative of the screening factor
75 : real(dp),intent(out) :: dscreendd !< on return, density derivative of the screening factor
76 : integer, intent(out) :: ierr
77 :
78 2079846 : real(dp) :: z1z2 !< z1*z2
79 2079846 : real(dp) :: t1, dt1dt, dt1dd
80 2079846 : real(dp) :: t2, dt2dt, dt2dd
81 2079846 : real(dp) :: t3, dt3dt, dt3dd !< the three terms in the fitting formula for h0fit
82 2079846 : real(dp) :: h0fit, dh0fitdt, dh0fitdd !< mean field fit
83 2079846 : real(dp) :: denom, ddenomdt, ddenomdd, denom2
84 2079846 : real(dp) :: tp, dtpdd, wp, dwpdd !< plasma temperature and frequency
85 2079846 : real(dp) :: tn, dtndd, dtndt !< normalised temperature (tk/tp)
86 2079846 : real(dp) :: beta, dbetadt, dbetadd, gama, dgamadt, dgamadd !< coeffs of zeta formula
87 2079846 : real(dp) :: gam, dgamdt, dgamdd
88 2079846 : real(dp) :: gamtild, dgamtilddt, dgamtilddd !< "effective" coulomb coupling parameter
89 2079846 : real(dp) :: gamtild2 !< "effective" coulomb coupling parameter squared
90 2079846 : real(dp) :: zeta, dzetadt, dzetadd !< function of plasma temperature
91 2079846 : real(dp) :: zeta2, zeta3
92 2079846 : real(dp) :: a_1, da_1dd, a_2,da_2dd !< ion sphere radii
93 2079846 : real(dp) :: a_av, da_avdd !< average ion sphere radii
94 2079846 : real(dp) :: a_e !< electron sphere radius
95 2079846 : real(dp) :: a_b !< bohr radius
96 2079846 : real(dp) :: rs !< ion sphere radius normalised to bohr radius
97 : !real(dp) :: m1, m2 !< ion masses
98 : !real(dp) :: n1, n2 !< number density of ions of types 1 and 2
99 2079846 : real(dp) :: ntot, dntotdd !< total number density of these ions
100 : !real(dp) :: s !< evaluated sigmoid function
101 2079846 : real(dp) :: mav !< average ion mass
102 2079846 : real(dp) :: A, B, C, U, dAdt, dAdd, dBdt, dBdd, dCdt, dCdd
103 2079846 : real(dp) :: temp, rho, abar, zbar, rr
104 2079846 : real(dp) :: alpha, dalphadgam, dbetadgam, dalphadtn,dbetadtn
105 2079846 : real(dp) :: tk, dtkdt, dtkdd
106 :
107 : ! check whether both reactants are charged ions
108 :
109 2079846 : screen = 1d0
110 2079846 : dscreendt = 0d0
111 2079846 : dscreendd = 0d0
112 2079846 : ierr = 0
113 :
114 : ! Must be charged ions
115 2079846 : if (z1 <= 0d0 .or. z2 <= 0d0) then
116 0 : return
117 : end if
118 :
119 : ! No free charges to screen things
120 2079846 : if(sc% zbar<=1d0) return
121 :
122 :
123 2079846 : rho = sc% den
124 2079846 : temp = sc% temp
125 2079846 : zbar = sc% zbar
126 2079846 : abar = sc% abar
127 2079846 : rr = sc% rr ! den/abar
128 :
129 : ! ion masses and number densities
130 2079846 : mav = abar * amu
131 : ! ntot = den / (amu*abar)
132 2079846 : ntot = sc% ntot
133 2079846 : dntotdd = 1d0/(amu*abar)
134 :
135 : ! ion and electron sphere radii (itoh 1979 eq 1-3)
136 : !a_e = pow((3.d0 /(pi4 * zbar * ntot)),x13)
137 2079846 : a_e = sc% a_e
138 :
139 2079846 : a_1 = a_e * pre_z(int(z1))%z1_3 !pow(z1,x13)
140 2079846 : a_2 = a_e * pre_z(int(z2))%z1_3 !pow(z2,x13)
141 :
142 2079846 : da_1dd = -x13 * a_1/rho
143 2079846 : da_2dd = -x13 * a_2/rho
144 :
145 2079846 : a_av = 0.5d0 * (a_1 + a_2)
146 2079846 : da_avdd = 0.5d0 * (da_1dd + da_2dd)
147 :
148 2079846 : z1z2 = z1 * z2
149 :
150 : ! bohr radius and normalised ion sphere radius
151 :
152 2079846 : a_b = rbohr/z1z2/ntot
153 2079846 : rs = a_av / a_b
154 :
155 : ! plasma frequency and temperature (chugunov 2007 eq 2)
156 :
157 2079846 : wp = sqrt((pi4 * z1z2 * qe * qe * ntot/mav))
158 2079846 : dwpdd = x12 * wp/rho
159 :
160 2079846 : tp = hbar*wp/kerg
161 2079846 : dtpdd = (hbar/kerg) * dwpdd
162 :
163 2079846 : tn = temp/tp
164 2079846 : dtndt = 1d0/tp
165 2079846 : dtndd = -tn/tp * dtpdd
166 :
167 : ! Revise temperature used if nearing low tp values
168 2079846 : tk = temp
169 2079846 : dtkdt = 1d0
170 2079846 : dtkdd = 0d0
171 2079846 : denom = tp * (tk * dtpdd - tp * dtkdd)/(tk * tk * tk)
172 :
173 : ! Blend out for low Tp values
174 2079846 : if ( tn <= tp1)then
175 0 : tk = tp1 * tp
176 0 : dtkdt = 0d0
177 0 : dtkdd = tp1 * dtpdd
178 0 : denom = 0d0 ! Set as zero otherwise floating point issues cause small error in the subtraction (tk * dtpdd - tp * dtkdd)
179 :
180 2079846 : else if (tn >tp1 .and. tn <= tp2) then
181 :
182 0 : alpha = 0.5d0 * (1d0 - cospi((tn-tp1)/deltatp))
183 0 : dalphadtn = 0.5d0 * (pi/deltatp) * sinpi((tn-tp1)/deltatp)
184 :
185 0 : beta = 1.d0 - alpha
186 0 : dbetadtn = -dalphadtn
187 :
188 0 : dtkdt = (dbetadtn * dtndt * tp1 * tp) + (dalphadtn * dtndt * tp2 * tp)
189 : dtkdd = (beta * tp1 * dtpdd) + (dbetadtn * dtndd * tp1 * tp) + &
190 0 : (alpha * tp2 * dtpdd) + (dalphadtn * dtndd * tp2 * tp)
191 :
192 0 : tk = beta * (tp1 * tp) + alpha * (tp2 * tp)
193 0 : denom = tp * (tk * dtpdd - tp * dtkdd)/(tk * tk * tk)
194 : end if
195 :
196 : ! zeta (chugunov 2007 eq 3)
197 2079846 : U = four_thirds*(tp * tp)/(pi2 * tk * tk)
198 2079846 : zeta = pow(U,x13)
199 2079846 : dzetadt = -x23 * (zeta/tk) * dtkdt
200 2079846 : dzetadd = x13 * (zeta/U) * (four_thirds/pi2) * 2.d0 * denom
201 :
202 : ! coulomb coupling parameter gamma (itoh 1979 eq 4)
203 2079846 : gam = z1z2 * qe * qe/(a_av * tk * kerg)
204 2079846 : dgamdt = -(gam / tk) * dtkdt
205 2079846 : dgamdd = -(gam/(a_av * tk)) * (tk * da_avdd + a_av * dtkdd)
206 :
207 2079846 : if (gam >=gamfitlim) then
208 :
209 : gam = gamfitlim
210 : dgamdt = 0d0
211 : dgamdd = 0d0
212 :
213 2079846 : else if (gam <= gamfitlim .and. gam >= g0 ) then
214 0 : alpha = 0.5d0 * (1d0 - cospi((gam-g0)/deltagam))
215 0 : dalphadgam = 0.5d0 * pi * sinpi((gam-g0)/deltagam) * (1d0/deltagam)
216 0 : beta = 1.d0 - alpha
217 0 : dbetadgam = -dalphadgam
218 :
219 0 : dgamdt = beta * dgamdt + (dbetadgam * dgamdt) * gam + (dalphadgam * dgamdt) * gamfitlim
220 0 : dgamdd = beta * dgamdd + (dbetadgam * dgamdd) * gam + (dalphadgam * dgamdd) * gamfitlim
221 :
222 0 : gam = beta * gam + alpha * gamfitlim
223 :
224 : end if
225 :
226 : ! coefficients of zeta dependent on the ion coulomb coupling
227 : ! parameter (gam) (chugunov 2007, just after eq 21)
228 :
229 2079846 : beta = 0.41d0 - 0.6d0/gam
230 2079846 : dbetadt = 0.6d0/(gam * gam) * dgamdt
231 2079846 : dbetadd = 0.6d0/(gam * gam) * dgamdd
232 :
233 2079846 : gama = 0.06d0 + 2.2d0/gam
234 2079846 : dgamadt = -2.2d0/(gam * gam) * dgamdt
235 2079846 : dgamadd = -2.2d0/(gam * gam) * dgamdd
236 :
237 : ! gamma tilda (chugunov 2007 eq 21)
238 2079846 : zeta2 = zeta * zeta
239 2079846 : zeta3 = zeta2 * zeta
240 :
241 2079846 : U = (1d0 + alfa * zeta + beta * zeta2 + gama * zeta3)
242 2079846 : denom = pow(U,-x13)
243 2079846 : denom2 = denom/U ! pow(U,-x43)
244 : ddenomdt = -x13 * denom2 * &
245 : ((alfa * dzetadt) + (zeta2 * dbetadt + 2d0 * dzetadt * beta * zeta) + &
246 2079846 : (zeta3 * dgamadt + 3d0 * dzetadt * gama * zeta2 ))
247 : ddenomdd = -x13 * denom2 * &
248 : ((alfa * dzetadd) + (zeta2 * dbetadd + 2d0 * dzetadd * beta * zeta) + &
249 2079846 : (zeta3 * dgamadd + 3d0 * dzetadd * gama * zeta2 ))
250 :
251 2079846 : gamtild = gam * denom
252 2079846 : dgamtilddt = ddenomdt * gam + dgamdt * denom
253 2079846 : dgamtilddd = ddenomdd * gam + dgamdd * denom
254 :
255 2079846 : gamtild2 = gamtild * gamtild
256 : ! for for mean field potential H (chugunov 2007 eq 19)
257 :
258 2079846 : A = gamtild * sqrt(gamtild)
259 2079846 : dAdt = x32 * (A/gamtild) * dgamtilddt
260 2079846 : dAdd = x32 * (A/gamtild) * dgamtilddd
261 :
262 2079846 : B = c_a1/sqrt(c_a2 + gamtild)
263 2079846 : dBdt = -1d0/2d0 * B/(c_a2 + gamtild) * dgamtilddt
264 2079846 : dBdd = -1d0/2d0 * B/(c_a2 + gamtild) * dgamtilddd
265 :
266 2079846 : C = c_a3/(1d0 + gamtild)
267 2079846 : dCdt = -C/(1d0 + gamtild) * dgamtilddt
268 2079846 : dCdd = -C/(1d0 + gamtild) * dgamtilddd
269 :
270 2079846 : t1 = A * ( B + C )
271 2079846 : dt1dt = A * ( dBdt + dCdt ) + dAdt * ( B + C )
272 2079846 : dt1dd = A * ( dBdd + dCdd ) + dAdd * ( B + C )
273 :
274 2079846 : denom = c_b2 + gamtild
275 2079846 : ddenomdt = dgamtilddt
276 2079846 : ddenomdd = dgamtilddd
277 2079846 : t2 = c_b1 * gamtild2 / denom
278 2079846 : dt2dt = (( denom * c_b1 * 2d0 * gamtild * dgamtilddt ) - (c_b1 * gamtild2 * ddenomdt)) / (denom * denom)
279 2079846 : dt2dd = (( denom * c_b1 * 2d0 * gamtild * dgamtilddd ) - (c_b1 * gamtild2 * ddenomdd))/ (denom * denom)
280 :
281 2079846 : denom = c_b4 + gamtild2
282 2079846 : ddenomdt = 2d0 * gamtild * dgamtilddt
283 2079846 : ddenomdd = 2d0 * gamtild * dgamtilddd
284 2079846 : t3 = c_b3 * gamtild2 / denom
285 2079846 : dt3dt = (( denom * c_b3 * 2d0 * gamtild * dgamtilddt ) - (c_b3 * gamtild2 * ddenomdt)) / (denom * denom)
286 2079846 : dt3dd = (( denom * c_b3 * 2d0 * gamtild * dgamtilddd ) - (c_b3 * gamtild2 * ddenomdd)) / (denom * denom)
287 :
288 2079846 : h0fit = t1 + t2 + t3
289 2079846 : dh0fitdt = dt1dt + dt2dt + dt3dt
290 2079846 : dh0fitdd = dt1dd + dt2dd + dt3dd
291 :
292 2079846 : dscreendt = dh0fitdt * screen
293 2079846 : dscreendd = dh0fitdd * screen
294 :
295 : ! limit screening factor to h0fitlim
296 2079846 : if (h0fit > h0fitlim) then
297 0 : h0fit = h0fitlim
298 0 : dh0fitdt = 0d0
299 0 : dh0fitdd = 0d0
300 : end if
301 :
302 2079846 : screen = exp(h0fit)
303 2079846 : dscreendt = dh0fitdt * screen
304 2079846 : dscreendd = dh0fitdd * screen
305 :
306 : end subroutine eval_screen_chugunov
307 :
308 : end module screening_chugunov
309 :
|
