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 439 : 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 439 : ierr = 0
50 :
51 439 : 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 439 : zs13 = pow(z1 + z2, x13)
62 439 : zs13inv = 1.0d0/zs13
63 439 : zhat = pow(z1 + z2, x53) - pow(z1,x53) - pow(z2,x53)
64 439 : zhat2 = pow(z1 + z2, x512) - pow(z1,x512) - pow(z2,x512)
65 439 : lzav = x53 * log(max(1d-99,z1*z2/(z1 + z2)))
66 439 : aznut = pow(z1*z1*z2*z2*a1*a2/(a1 + a2), x13)
67 :
68 : end subroutine screen5_init_AZ_info
69 :
70 173149 : 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 : real(dp) :: aa, daadt, daadd, bb, cc, dccdt, dccdd, &
112 : qq, dqqdt, dqqdd, rr, drrdt, drrdd, &
113 : ss, dssdt, dssdd, tt, dttdt, dttdd, uu, duudt, duudd, &
114 : vv, dvvdt, dvvdd, a3, da3, &
115 : qlam0z, qlam0zdt, qlam0zdd, &
116 : h12w, dh12wdt, dh12wdd, h12, dh12dt, dh12dd, &
117 : taufac, taufacdt, gamp, gampdt, gampdd, &
118 : gamef, gamefdt, gamefdd, &
119 : tau12, tau12dt, alph12, alph12dt, alph12dd, &
120 : xlgfac, dxlgfacdt, dxlgfacdd, &
121 : gamp14, gamp14dt, gamp14dd, h12x, dh12xdt, dh12xdd, &
122 : gamefx, gamefs, temp, den, zbar, abar, z2bar,&
123 : 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 173149 : ierr = 0
152 :
153 173149 : 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 173149 : zbar = sc% zbar
161 173149 : abar = sc% abar
162 173149 : z2bar = sc% z2bar
163 173149 : temp = sc% temp
164 173149 : den = sc% den
165 :
166 : !..unload the screen info
167 : !ytot = sc% ytot
168 173149 : 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 173149 : qlam0z = sc% qlam0z
178 173149 : qlam0zdt = sc% qlam0zdt
179 173149 : qlam0zdd = sc% qlam0zdd
180 :
181 173149 : taufac = sc% taufac
182 173149 : taufacdt = sc% taufacdt
183 :
184 : !xni = sc% xni
185 : !dxnidd = sc% dxnidd
186 :
187 173149 : aa = sc% aa
188 173149 : daadt = sc% daadt
189 173149 : daadd = sc% daadd
190 :
191 :
192 : !..calculate individual screening factors
193 173149 : bb = z1 * z2
194 173149 : gamp = aa
195 173149 : gampdt = daadt
196 173149 : gampdd = daadd
197 :
198 173149 : qq = fact * bb * zs13inv
199 173149 : gamef = qq * gamp
200 173149 : gamefdt = qq * gampdt
201 173149 : gamefdd = qq * gampdd
202 :
203 173149 : tau12 = taufac * aznut
204 173149 : tau12dt = taufacdt * aznut
205 :
206 173149 : qq = 1.0d0/tau12
207 173149 : alph12 = gamef * qq
208 173149 : alph12dt = (gamefdt - alph12*tau12dt) * qq
209 173149 : 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 173149 : if (alph12 > alph12_lim) then
214 :
215 38 : alph12 = alph12_lim
216 38 : alph12dt = 0.0d0
217 38 : alph12dd = 0.0d0
218 :
219 38 : gamef = alph12 * tau12
220 38 : gamefdt = alph12 * tau12dt
221 38 : gamefdd = 0.0d0
222 :
223 38 : qq = zs13/(fact * bb)
224 38 : gamp = gamef * qq
225 38 : gampdt = gamefdt * qq
226 38 : gampdd = 0.0d0
227 :
228 : end if
229 :
230 :
231 : !..weak screening regime
232 173149 : h12w = bb * qlam0z
233 173149 : dh12wdt = bb * qlam0zdt
234 173149 : dh12wdd = bb * qlam0zdd
235 :
236 173149 : h12 = h12w
237 173149 : dh12dt = dh12wdt
238 173149 : 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 173149 : gamefx = 0.3d0
252 173149 : if (gamef > gamefx) then
253 : if (debug) write(*,1) 'intermediate and strong sceening regime'
254 :
255 59119 : gamp14 = pow(gamp,x14)
256 59119 : rr = 1.0d0/gamp
257 59119 : qq = 0.25d0*gamp14*rr
258 59119 : gamp14dt = qq * gampdt
259 59119 : gamp14dd = qq * gampdd
260 :
261 : cc = 0.896434d0 * gamp * zhat &
262 : - 3.44740d0 * gamp14 * zhat2 &
263 : - 0.5551d0 * (log(gamp) + lzav) &
264 59119 : - 2.996d0
265 :
266 : dccdt = 0.896434d0 * gampdt * zhat &
267 : - 3.44740d0 * gamp14dt * zhat2 &
268 59119 : - 0.5551d0*rr*gampdt
269 :
270 : dccdd = 0.896434d0 * gampdd * zhat &
271 : - 3.44740d0 * gamp14dd * zhat2 &
272 59119 : - 0.5551d0*rr*gampdd
273 :
274 59119 : a3 = alph12 * alph12 * alph12
275 59119 : da3 = 3.0d0 * alph12 * alph12
276 :
277 59119 : qq = 0.014d0 + 0.0128d0*alph12
278 59119 : dqqdt = 0.0128d0*alph12dt
279 59119 : dqqdd = 0.0128d0*alph12dd
280 :
281 59119 : rr = x532 - alph12*qq
282 59119 : drrdt = -(alph12dt*qq + alph12*dqqdt)
283 59119 : drrdd = -(alph12dd*qq + alph12*dqqdd)
284 :
285 59119 : ss = tau12*rr
286 59119 : dssdt = tau12dt*rr + tau12*drrdt
287 59119 : dssdd = tau12*drrdd
288 :
289 59119 : tt = -0.0098d0 + 0.0048d0*alph12
290 59119 : dttdt = 0.0048d0*alph12dt
291 59119 : dttdd = 0.0048d0*alph12dd
292 :
293 59119 : uu = 0.0055d0 + alph12*tt
294 59119 : duudt = alph12dt*tt + alph12*dttdt
295 59119 : duudd = alph12dd*tt + alph12*dttdd
296 :
297 59119 : vv = gamef * alph12 * uu
298 59119 : dvvdt= gamefdt*alph12*uu + gamef*alph12dt*uu + gamef*alph12*duudt
299 59119 : dvvdd= gamefdd*alph12*uu + gamef*alph12dd*uu + gamef*alph12*duudd
300 :
301 59119 : h12 = cc - a3 * (ss + vv)
302 59119 : rr = da3 * (ss + vv)
303 59119 : dh12dt = dccdt - rr*alph12dt - a3*(dssdt + dvvdt)
304 59119 : dh12dd = dccdd - rr*alph12dd - a3*(dssdd + dvvdd)
305 :
306 59119 : rr = 1.0d0 - 0.0562d0*a3
307 59119 : ss = -0.0562d0*da3
308 59119 : drrdt = ss*alph12dt
309 59119 : drrdd = ss*alph12dd
310 :
311 :
312 : if (debug) write(*, 1) 'rr', rr
313 :
314 59119 : if (rr >= 0.77d0) then
315 59081 : xlgfac = rr
316 59081 : dxlgfacdt = drrdt
317 59081 : dxlgfacdd = drrdd
318 : else
319 38 : xlgfac = 0.77d0
320 38 : dxlgfacdt = 0.0d0
321 38 : dxlgfacdd = 0.0d0
322 : end if
323 :
324 59119 : h12 = log(xlgfac) + h12
325 59119 : rr = 1.0d0/xlgfac
326 59119 : dh12dt = rr*dxlgfacdt + dh12dt
327 59119 : dh12dd = rr*dxlgfacdd + dh12dd
328 :
329 :
330 :
331 : if (debug) write(*, 1) 'gamef', gamef
332 :
333 59119 : gamefs = 0.8d0
334 59119 : if (gamef <= gamefs) then
335 5690 : dgamma = 1.0d0/(gamefs - gamefx)
336 :
337 5690 : rr = dgamma*(gamefs - gamef)
338 5690 : drrdt = -dgamma*gamefdt
339 5690 : drrdd = -dgamma*gamefdd
340 :
341 5690 : ss = dgamma*(gamef - gamefx)
342 5690 : dssdt = dgamma*gamefdt
343 5690 : dssdd = dgamma*gamefdd
344 :
345 5690 : vv = h12
346 :
347 5690 : h12x = h12
348 5690 : dh12xdt = dh12dt
349 5690 : dh12xdd = dh12dd
350 :
351 5690 : h12 = h12w*rr + vv*ss
352 5690 : dh12dt = dh12wdt*rr + h12w*drrdt + dh12dt*ss + vv*dssdt
353 5690 : 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 173149 : h12 = max(min(h12, h12_max), 0.0d0)
378 173149 : scor = exp(h12)
379 173149 : if (h12 == h12_max) then
380 25 : scordt = 0.0d0
381 25 : scordd = 0.0d0
382 : else
383 173124 : scordt = scor * dh12dt
384 173124 : 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
|
