Line data Source code
1 : ! ***********************************************************************
2 : !
3 : ! Copyright (C) 2010-2019 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 screen5
21 :
22 : use rates_def
23 : use const_def, only: ln10, pi, two_13
24 : use math_lib
25 :
26 : implicit none
27 :
28 : contains
29 :
30 852 : subroutine screen5_init_AZ_info( &
31 : zs13, zhat, zhat2, lzav, aznut, zs13inv, a1, z1, a2, z2, ierr)
32 : !..compute and store things that only depend on reaction A's and Z's
33 : real(dp), intent(out) :: zs13, zhat, zhat2, lzav, aznut, zs13inv
34 : ! zs13 = (z1+z2)**(1./3.)
35 : ! zhat = combination of z1 and z2 raised to the 5/3 power
36 : ! zhat2 = combination of z1 and z2 raised to the 5/12 power
37 : ! lzav = log of effective charge
38 : ! aznut = combination of a1, z1, a2, z2 raised to 1/3 power
39 : ! zs13inv = 1 / zs13
40 : real(dp), intent(in) :: a1, z1, a2, z2
41 : integer, intent(out) :: ierr
42 :
43 : real(dp), parameter :: x13 = 1.0d0/3.0d0
44 : real(dp), parameter :: x14 = 1.0d0/4.0d0
45 : real(dp), parameter :: x53 = 5.0d0/3.0d0
46 : real(dp), parameter :: x532 = 5.0d0/32.0d0
47 : real(dp), parameter :: x512 = 5.0d0/12.0d0
48 :
49 852 : ierr = 0
50 :
51 852 : if (z1 <= 0 .or. z2 <= 0) then
52 0 : zs13 = 0
53 0 : zs13inv = 0
54 0 : zhat = 0
55 0 : zhat2 = 0
56 0 : lzav = 0
57 0 : aznut = 0
58 0 : return
59 : end if
60 :
61 852 : zs13 = pow(z1 + z2, x13)
62 852 : zs13inv = 1.0d0/zs13
63 852 : zhat = pow(z1 + z2, x53) - pow(z1,x53) - pow(z2,x53)
64 852 : zhat2 = pow(z1 + z2, x512) - pow(z1,x512) - pow(z2,x512)
65 852 : lzav = x53 * log(max(1d-99,z1*z2/(z1 + z2)))
66 852 : aznut = pow(z1*z1*z2*z2*a1*a2/(a1 + a2), x13)
67 :
68 : end subroutine screen5_init_AZ_info
69 :
70 346298 : subroutine fxt_screen5(sc, zs13, zhat, zhat2, lzav, aznut, zs13inv, &
71 : a1, z1, a2, z2, scor, scordt, scordd, ierr)
72 : !..this subroutine calculates screening factors and their derivatives
73 : !..for nuclear reaction rates in the weak, intermediate and strong regimes.
74 :
75 : !..based on graboske, dewit, grossman and cooper apj 181 457 1973 for weak screening.
76 :
77 : !..based on alastuey and jancovici apj 226 1034 1978,
78 : !..with plasma parameters from itoh et al apj 234 1079 1979, for strong screening.
79 :
80 : !..input:
81 : !..temp = temperature
82 : !..den = density
83 : !..zbar = mean charge per nucleus
84 : !..abar = mean number of nucleons per nucleus
85 : !..z2bar = mean square charge per nucleus
86 : !..z1 a1 = charge and number in the entrance channel
87 : !..z2 a2 = charge and number in the exit channel
88 : !..jscreen = counter of which reaction is being calculated
89 : !..init = flag to compute the more expensive functions just once
90 :
91 : !..output:
92 : !..scor = screening correction
93 : !..scordt = derivative of screening correction with temperature
94 : !..scordd = derivative of screening correction with density
95 :
96 :
97 : !..declare the pass
98 : type (Screen_Info) :: sc
99 : real(dp), intent(in) :: zs13, zhat, zhat2, lzav, aznut, zs13inv
100 : ! zs13 = (z1+z2)**(1./3.)
101 : ! zhat = combination of z1 and z2 raised to the 5/3 power
102 : ! zhat2 = combination of z1 and z2 raised to the 5/12 power
103 : ! lzav = log of effective charge
104 : ! aznut = combination of a1, z1, a2, z2 raised to 1/3 power
105 : ! zs13inv = 1 / zs13
106 : real(dp), intent(in) :: a1, z1, a2, z2
107 : real(dp), intent(out) :: scor, scordt, scordd
108 : integer, intent(out) :: ierr
109 :
110 : !..local variables
111 346298 : real(dp) :: aa, daadt, daadd, bb, cc, dccdt, dccdd, &
112 346298 : qq, dqqdt, dqqdd, rr, drrdt, drrdd, &
113 346298 : ss, dssdt, dssdd, tt, dttdt, dttdd, uu, duudt, duudd, &
114 346298 : vv, dvvdt, dvvdd, a3, da3, &
115 346298 : qlam0z, qlam0zdt, qlam0zdd, &
116 346298 : h12w, dh12wdt, dh12wdd, h12, dh12dt, dh12dd, &
117 346298 : taufac, taufacdt, gamp, gampdt, gampdd, &
118 346298 : gamef, gamefdt, gamefdd, &
119 346298 : tau12, tau12dt, alph12, alph12dt, alph12dd, &
120 346298 : xlgfac, dxlgfacdt, dxlgfacdd, &
121 346298 : gamp14, gamp14dt, gamp14dd, h12x, dh12xdt, dh12xdd, &
122 346298 : gamefx, gamefs, temp, den, zbar, abar, z2bar,&
123 346298 : dgamma
124 :
125 : !..screening variables
126 : !..zs13 = (z1+z2)**(1./3.)
127 : !..zhat = combination of z1 and z2 raised to the 5/3 power
128 : !..zhat2 = combination of z1 and z2 raised to the 5/12 power
129 : !..lzav = log of effective charge
130 : !..aznut = combination of a1, z1, a2, z2 raised to 1/3 power
131 :
132 :
133 : real(dp), parameter :: alph12_lim = 1.6d0 ! ln(10)
134 : real(dp), parameter :: h12_max = 300d0
135 : real(dp), parameter :: x13 = 1.0d0/3.0d0
136 : real(dp), parameter :: x14 = 1.0d0/4.0d0
137 : real(dp), parameter :: x53 = 5.0d0/3.0d0
138 : real(dp), parameter :: x532 = 5.0d0/32.0d0
139 : real(dp), parameter :: x512 = 5.0d0/12.0d0
140 : real(dp), parameter :: fact = two_13 ! the cube root of 2
141 : real(dp), parameter :: co2 = x13 * 4.248719d3
142 :
143 : logical, parameter :: debug = .false.
144 : !logical :: debug
145 :
146 : include 'formats'
147 :
148 :
149 : !debug = (abs(a1 - 4d0) < 1d-14 .and. abs(a2 - 12d0) < 1d-14)
150 :
151 346298 : ierr = 0
152 :
153 346298 : if (z1 <= 0d0 .or. z2 <= 0d0) then
154 0 : scor = 1d0
155 0 : scordt = 0d0
156 0 : scordd = 0d0
157 0 : return
158 : end if
159 :
160 346298 : zbar = sc% zbar
161 346298 : abar = sc% abar
162 346298 : z2bar = sc% z2bar
163 346298 : temp = sc% temp
164 346298 : den = sc% den
165 :
166 : !..unload the screen info
167 : !ytot = sc% ytot
168 346298 : rr = sc% rr
169 : !tempi = sc% tempi
170 : !dtempi = sc% dtempi
171 : !deni = sc% deni
172 :
173 : !pp = sc% pp
174 : !dppdt = sc% dppdt
175 : !dppdd = sc% dppdd
176 :
177 346298 : qlam0z = sc% qlam0z
178 346298 : qlam0zdt = sc% qlam0zdt
179 346298 : qlam0zdd = sc% qlam0zdd
180 :
181 346298 : taufac = sc% taufac
182 346298 : taufacdt = sc% taufacdt
183 :
184 : !xni = sc% xni
185 : !dxnidd = sc% dxnidd
186 :
187 346298 : aa = sc% aa
188 346298 : daadt = sc% daadt
189 346298 : daadd = sc% daadd
190 :
191 :
192 : !..calculate individual screening factors
193 346298 : bb = z1 * z2
194 346298 : gamp = aa
195 346298 : gampdt = daadt
196 346298 : gampdd = daadd
197 :
198 346298 : qq = fact * bb * zs13inv
199 346298 : gamef = qq * gamp
200 346298 : gamefdt = qq * gampdt
201 346298 : gamefdd = qq * gampdd
202 :
203 346298 : tau12 = taufac * aznut
204 346298 : tau12dt = taufacdt * aznut
205 :
206 346298 : qq = 1.0d0/tau12
207 346298 : alph12 = gamef * qq
208 346298 : alph12dt = (gamefdt - alph12*tau12dt) * qq
209 346298 : alph12dd = gamefdd * qq
210 :
211 : !..limit alph12 to alph12_lim to prevent unphysical behavior.
212 : !..this should really be replaced by a pycnonuclear reaction rate formula
213 346298 : if (alph12 > alph12_lim) then
214 :
215 76 : alph12 = alph12_lim
216 76 : alph12dt = 0.0d0
217 76 : alph12dd = 0.0d0
218 :
219 76 : gamef = alph12 * tau12
220 76 : gamefdt = alph12 * tau12dt
221 76 : gamefdd = 0.0d0
222 :
223 76 : qq = zs13/(fact * bb)
224 76 : gamp = gamef * qq
225 76 : gampdt = gamefdt * qq
226 76 : gampdd = 0.0d0
227 :
228 : end if
229 :
230 :
231 : !..weak screening regime
232 346298 : h12w = bb * qlam0z
233 346298 : dh12wdt = bb * qlam0zdt
234 346298 : dh12wdd = bb * qlam0zdd
235 :
236 346298 : h12 = h12w
237 346298 : dh12dt = dh12wdt
238 346298 : dh12dd = dh12wdd
239 :
240 : !..intermediate and strong sceening regime
241 :
242 : if (debug) write(*, 1) 'gamef', gamef
243 : if (debug) write(*, 1) 'fact', fact
244 : if (debug) write(*, 1) 'z1', z1
245 : if (debug) write(*, 1) 'z2', z2
246 : if (debug) write(*, 1) 'gamp', gamp
247 : if (debug) write(*, 1) 'sc% aa', sc% aa
248 : if (debug) write(*, 1) 'sc% tempi', sc% tempi
249 : if (debug) write(*, 1) 'sc% xni', sc% xni
250 :
251 346298 : gamefx = 0.3d0
252 346298 : if (gamef > gamefx) then
253 : if (debug) write(*,1) 'intermediate and strong sceening regime'
254 :
255 118238 : gamp14 = pow(gamp,x14)
256 118238 : rr = 1.0d0/gamp
257 118238 : qq = 0.25d0*gamp14*rr
258 118238 : gamp14dt = qq * gampdt
259 118238 : gamp14dd = qq * gampdd
260 :
261 : cc = 0.896434d0 * gamp * zhat &
262 : - 3.44740d0 * gamp14 * zhat2 &
263 : - 0.5551d0 * (log(gamp) + lzav) &
264 118238 : - 2.996d0
265 :
266 : dccdt = 0.896434d0 * gampdt * zhat &
267 : - 3.44740d0 * gamp14dt * zhat2 &
268 118238 : - 0.5551d0*rr*gampdt
269 :
270 : dccdd = 0.896434d0 * gampdd * zhat &
271 : - 3.44740d0 * gamp14dd * zhat2 &
272 118238 : - 0.5551d0*rr*gampdd
273 :
274 118238 : a3 = alph12 * alph12 * alph12
275 118238 : da3 = 3.0d0 * alph12 * alph12
276 :
277 118238 : qq = 0.014d0 + 0.0128d0*alph12
278 118238 : dqqdt = 0.0128d0*alph12dt
279 118238 : dqqdd = 0.0128d0*alph12dd
280 :
281 118238 : rr = x532 - alph12*qq
282 118238 : drrdt = -(alph12dt*qq + alph12*dqqdt)
283 118238 : drrdd = -(alph12dd*qq + alph12*dqqdd)
284 :
285 118238 : ss = tau12*rr
286 118238 : dssdt = tau12dt*rr + tau12*drrdt
287 118238 : dssdd = tau12*drrdd
288 :
289 118238 : tt = -0.0098d0 + 0.0048d0*alph12
290 118238 : dttdt = 0.0048d0*alph12dt
291 118238 : dttdd = 0.0048d0*alph12dd
292 :
293 118238 : uu = 0.0055d0 + alph12*tt
294 118238 : duudt = alph12dt*tt + alph12*dttdt
295 118238 : duudd = alph12dd*tt + alph12*dttdd
296 :
297 118238 : vv = gamef * alph12 * uu
298 118238 : dvvdt= gamefdt*alph12*uu + gamef*alph12dt*uu + gamef*alph12*duudt
299 118238 : dvvdd= gamefdd*alph12*uu + gamef*alph12dd*uu + gamef*alph12*duudd
300 :
301 118238 : h12 = cc - a3 * (ss + vv)
302 118238 : rr = da3 * (ss + vv)
303 118238 : dh12dt = dccdt - rr*alph12dt - a3*(dssdt + dvvdt)
304 118238 : dh12dd = dccdd - rr*alph12dd - a3*(dssdd + dvvdd)
305 :
306 118238 : rr = 1.0d0 - 0.0562d0*a3
307 118238 : ss = -0.0562d0*da3
308 118238 : drrdt = ss*alph12dt
309 118238 : drrdd = ss*alph12dd
310 :
311 :
312 : if (debug) write(*, 1) 'rr', rr
313 :
314 118238 : if (rr >= 0.77d0) then
315 118162 : xlgfac = rr
316 118162 : dxlgfacdt = drrdt
317 118162 : dxlgfacdd = drrdd
318 : else
319 76 : xlgfac = 0.77d0
320 76 : dxlgfacdt = 0.0d0
321 76 : dxlgfacdd = 0.0d0
322 : end if
323 :
324 118238 : h12 = log(xlgfac) + h12
325 118238 : rr = 1.0d0/xlgfac
326 118238 : dh12dt = rr*dxlgfacdt + dh12dt
327 118238 : dh12dd = rr*dxlgfacdd + dh12dd
328 :
329 :
330 :
331 : if (debug) write(*, 1) 'gamef', gamef
332 :
333 118238 : gamefs = 0.8d0
334 118238 : if (gamef <= gamefs) then
335 11380 : dgamma = 1.0d0/(gamefs - gamefx)
336 :
337 11380 : rr = dgamma*(gamefs - gamef)
338 11380 : drrdt = -dgamma*gamefdt
339 11380 : drrdd = -dgamma*gamefdd
340 :
341 11380 : ss = dgamma*(gamef - gamefx)
342 11380 : dssdt = dgamma*gamefdt
343 11380 : dssdd = dgamma*gamefdd
344 :
345 11380 : vv = h12
346 :
347 11380 : h12x = h12
348 11380 : dh12xdt = dh12dt
349 11380 : dh12xdd = dh12dd
350 :
351 11380 : h12 = h12w*rr + vv*ss
352 11380 : dh12dt = dh12wdt*rr + h12w*drrdt + dh12dt*ss + vv*dssdt
353 11380 : dh12dd = dh12wdd*rr + h12w*drrdd + dh12dd*ss + vv*dssdd
354 :
355 : if (debug) write(*, 1) 'gamef', gamef
356 : if (debug) write(*, 1) 'gamefx', gamefx
357 :
358 : if (debug) write(*,*) 'intermediate screening'
359 :
360 : else
361 : if (debug) write(*,*) 'strong screening'
362 :
363 : end if
364 :
365 : !..end of intermediate and strong screening if
366 :
367 : else if (debug) then
368 : write(*,*) 'weak screening'
369 :
370 : end if
371 :
372 : if (debug) write(*, 1) 'h12', h12
373 : if (debug) write(*, 1) 'h12/ln10', h12/ln10
374 : if (debug) write(*, *)
375 :
376 : !..machine limit the output
377 346298 : h12 = max(min(h12, h12_max), 0.0d0)
378 346298 : scor = exp(h12)
379 346298 : if (h12 == h12_max) then
380 50 : scordt = 0.0d0
381 50 : scordd = 0.0d0
382 : else
383 346248 : scordt = scor * dh12dt
384 346248 : scordd = scor * dh12dd
385 : end if
386 :
387 : if (debug) write(*, 1) 'scor', scor
388 : if (debug) write(*, 1) 'scordt', scordt
389 : if (debug) write(*, 1) 'scordd', scordd
390 : if (debug) write(*, *)
391 : if (debug) write(*, *)
392 : if (debug) write(*, *)
393 :
394 : end subroutine fxt_screen5
395 :
396 : end module screen5
|
