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 mod_neu
21 : use const_def, only: dp, pi, ln10, weinberg_theta, num_neu_fam, iln10, one_third, two_thirds, one_sixth, arg_not_provided
22 : use neu_def
23 : use math_lib
24 : use utils_lib, only: mesa_error
25 :
26 : implicit none
27 :
28 : !..various numerical constants
29 :
30 : real(dp), parameter :: con1 = 1.0d0/5.9302d0
31 :
32 : !..cv and ca are the vector and axial currents.
33 :
34 : real(dp), parameter :: cv = 0.5d0 + 2.0d0 * weinberg_theta
35 : real(dp), parameter :: cvp = 1.0d0 - cv
36 : real(dp), parameter :: ca = 0.5d0
37 : real(dp), parameter :: cap = 1.0d0 - ca
38 : real(dp), parameter :: tfac1 = cv*cv + ca*ca + (num_neu_fam-1.0d0) * (cvp*cvp+cap*cap)
39 : real(dp), parameter :: tfac2 = cv*cv - ca*ca + (num_neu_fam-1.0d0) * (cvp*cvp - cap*cap)
40 : real(dp), parameter :: tfac3 = tfac2/tfac1
41 : real(dp), parameter :: tfac4 = 0.5d0 * tfac1
42 : real(dp), parameter :: tfac5 = 0.5d0 * tfac2
43 : real(dp), parameter :: tfac6 = cv*cv + 1.5d0*ca*ca + (num_neu_fam - 1.0d0)*(cvp*cvp + 1.5d0*cap*cap)
44 :
45 : real(dp), parameter :: fac1 = 5.0d0 * pi / 3.0d0
46 : real(dp), parameter :: fac2 = 10.0d0 * pi
47 : real(dp), parameter :: fac3 = pi / 5.0d0
48 :
49 :
50 : type t8s
51 : real(dp) :: t8,t812,t832,t82,t83,t85,t86,t8m1,t8m2,t8m3,t8m5,t8m6
52 : end type t8s
53 :
54 :
55 : type inputs
56 : real(dp) :: temp,logtemp,den,logden,den6
57 : real(dp) :: ye,deni,tempi,abari,zbari,abar,zbar
58 : real(dp) :: t9,xl,xldt,xlp5,xl2,xl3,xl4,xl5,xl6,xl7,xl8,xl9,xlmp5,xlm1,xlm2,xlm3,xlm4
59 : real(dp) :: rm,rmdd,rmda,rmdz,rmi
60 : real(dp) :: zeta,zetadt,zetadd,zetada,zetadz,zeta2,zeta3
61 : end type inputs
62 :
63 :
64 : contains
65 :
66 33684 : real(dp) function ifermi12(f)
67 :
68 : !..this routine applies a rational function expansion to get the inverse
69 : !..fermi-dirac integral of order 1/2 when it is equal to f.
70 : !..maximum error is 4.19d-9. reference: antia apjs 84,101 1993
71 :
72 : !..declare
73 : real(dp), intent(in) :: f
74 : integer :: i,m1,k1,m2,k2
75 33684 : real(dp) :: an,a1(12),b1(12),a2(12),b2(12),rn,den,ff
76 :
77 :
78 : !..load the coefficients of the expansion
79 : data an,m1,k1,m2,k2 /0.5d0, 4, 3, 6, 5/
80 : data (a1(i),i=1,5)/ 1.999266880833d4, 5.702479099336d3, &
81 : 6.610132843877d2, 3.818838129486d1, &
82 : 1.0d0/
83 : data (b1(i),i=1,4)/ 1.771804140488d4, -2.014785161019d3, &
84 : 9.130355392717d1, -1.670718177489d0/
85 : data (a2(i),i=1,7)/-1.277060388085d-2, 7.187946804945d-2, &
86 : -4.262314235106d-1, 4.997559426872d-1, &
87 : -1.285579118012d0, -3.930805454272d-1, &
88 : 1.0d0/
89 : data (b2(i),i=1,6)/-9.745794806288d-3, 5.485432756838d-2, &
90 : -3.299466243260d-1, 4.077841975923d-1, &
91 : -1.145531476975d0, -6.067091689181d-2/
92 :
93 :
94 33684 : if (f < 4.0d0) then
95 27774 : rn = f + a1(m1)
96 111096 : do i=m1-1,1,-1
97 111096 : rn = rn*f + a1(i)
98 : end do
99 27774 : den = b1(k1+1)
100 111096 : do i=k1,1,-1
101 111096 : den = den*f + b1(i)
102 : end do
103 27774 : ifermi12 = log(f * rn/den)
104 :
105 : else
106 5910 : ff = 1.0d0/pow(f,1.0d0/(1.0d0 + an))
107 5910 : rn = ff + a2(m2)
108 35460 : do i=m2-1,1,-1
109 35460 : rn = rn*ff + a2(i)
110 : end do
111 5910 : den = b2(k2+1)
112 35460 : do i=k2,1,-1
113 35460 : den = den*ff + b2(i)
114 : end do
115 5910 : ifermi12 = rn/(den*ff)
116 : end if
117 :
118 33684 : end function ifermi12
119 :
120 :
121 33684 : real(dp) function zfermim12(x)
122 :
123 : !..this routine applies a rational function expansion to get the fermi-dirac
124 : !..integral of order -1/2 evaluated at x. maximum error is 1.23d-12.
125 : !..reference: antia apjs 84,101 1993
126 :
127 : !..declare
128 : real(dp), intent(in) :: x
129 : integer :: i,m1,k1,m2,k2
130 33684 : real(dp) :: an,a1(12),b1(12),a2(12),b2(12),rn,den,xx
131 :
132 : !..load the coefficients of the expansion
133 : data an,m1,k1,m2,k2 /-0.5d0, 7, 7, 11, 11/
134 : data (a1(i),i=1,8)/ 1.71446374704454d7, 3.88148302324068d7, &
135 : 3.16743385304962d7, 1.14587609192151d7, &
136 : 1.83696370756153d6, 1.14980998186874d5, &
137 : 1.98276889924768d3, 1.0d0/
138 : data (b1(i),i=1,8)/ 9.67282587452899d6, 2.87386436731785d7, &
139 : 3.26070130734158d7, 1.77657027846367d7, &
140 : 4.81648022267831d6, 6.13709569333207d5, &
141 : 3.13595854332114d4, 4.35061725080755d2/
142 : data (a2(i),i=1,12)/-4.46620341924942d-15, -1.58654991146236d-12, &
143 : -4.44467627042232d-10, -6.84738791621745d-8, &
144 : -6.64932238528105d-6, -3.69976170193942d-4, &
145 : -1.12295393687006d-2, -1.60926102124442d-1, &
146 : -8.52408612877447d-1, -7.45519953763928d-1, &
147 : 2.98435207466372d0, 1.0d0/
148 : data (b2(i),i=1,12)/-2.23310170962369d-15, -7.94193282071464d-13, &
149 : -2.22564376956228d-10, -3.43299431079845d-8, &
150 : -3.33919612678907d-6, -1.86432212187088d-4, &
151 : -5.69764436880529d-3, -8.34904593067194d-2, &
152 : -4.78770844009440d-1, -4.99759250374148d-1, &
153 : 1.86795964993052d0, 4.16485970495288d-1/
154 :
155 :
156 33684 : if (x < 2.0d0) then
157 27610 : xx = exp(x)
158 27610 : rn = xx + a1(m1)
159 193270 : do i=m1-1,1,-1
160 193270 : rn = rn*xx + a1(i)
161 : end do
162 27610 : den = b1(k1+1)
163 220880 : do i=k1,1,-1
164 220880 : den = den*xx + b1(i)
165 : end do
166 27610 : zfermim12 = xx * rn/den
167 : else
168 6074 : xx = 1.0d0/(x*x)
169 6074 : rn = xx + a2(m2)
170 66814 : do i=m2-1,1,-1
171 66814 : rn = rn*xx + a2(i)
172 : end do
173 6074 : den = b2(k2+1)
174 72888 : do i=k2,1,-1
175 72888 : den = den*xx + b2(i)
176 : end do
177 6074 : zfermim12 = sqrt(x)*rn/den
178 : end if
179 :
180 33684 : end function zfermim12
181 :
182 :
183 70368 : subroutine neutrinos(T, logT, Rho, logRho, abar, zbar, log10_Tlim, &
184 : flags, loss, sources, info)
185 : use utils_lib, only: is_bad
186 :
187 : !..this routine computes neutrino losses from the analytic fits of
188 : !..itoh et al. apjs 102, 411, 1996, and also returns their derivatives.
189 :
190 : ! provide T or logT or both (the code needs both, so pass 'em if you've got 'em!)
191 : ! same for Rho and logRho
192 :
193 : real(dp), intent(in) :: T ! temperature
194 : real(dp), intent(in) :: logT ! log10 of temperature
195 : real(dp), intent(in) :: Rho ! density
196 : real(dp), intent(in) :: logRho ! log10 of density
197 : real(dp), intent(in) :: abar ! mean atomic weight
198 : real(dp), intent(in) :: zbar ! mean charge
199 : real(dp), intent(in) :: log10_Tlim ! start to cutoff at this temperature
200 : logical, intent(in) :: flags(num_neu_types) ! true if should include the type
201 : ! in most cases for stellar evolution, you may want to include brem, plas, pair, and phot
202 : ! but skip reco. it is fairly expensive to compute and typically makes only a small contribution.
203 :
204 : real(dp), intent(inout) :: loss(num_neu_rvs) ! total from all sources
205 : real(dp), intent(inout) :: sources(num_neu_types, num_neu_rvs)
206 : integer, intent(out) :: info ! 0 means AOK.
207 :
208 : !..local variables
209 : type(inputs) :: input
210 :
211 36684 : real(dp) :: temp,logtemp,den,logden,tcutoff_factor
212 :
213 36684 : real(dp) :: snu, dsnudt, dsnudd, dsnuda, dsnudz, dtcutoff_factordt, dtlim
214 36684 : real(dp) :: spair,spairdt,spairdd,spairda,spairdz, &
215 36684 : splas,splasdt,splasdd,splasda,splasdz, &
216 36684 : sphot,sphotdt,sphotdd,sphotda,sphotdz, &
217 36684 : sbrem,sbremdt,sbremdd,sbremda,sbremdz, &
218 36684 : sreco,srecodt,srecodd,srecoda,srecodz
219 :
220 :
221 36684 : info = 0
222 :
223 36684 : if ((T /= arg_not_provided .and. T <= Tmin_neu) .or. &
224 : (logT /= arg_not_provided .and. logT <= log10Tmin_neu)) then
225 3000 : loss = 0d0
226 3000 : sources(1:num_neu_types, 1:num_neu_rvs) = 0d0
227 3000 : return
228 : end if
229 :
230 33684 : if (T == arg_not_provided) then
231 0 : if (logT == arg_not_provided) then
232 0 : info = -2
233 0 : return
234 : end if
235 0 : temp = exp10(logT)
236 : else
237 33684 : temp = T
238 : end if
239 :
240 33684 : if (T <= 0) then
241 0 : info = -1
242 0 : return
243 : end if
244 :
245 33684 : if (logT == arg_not_provided) then
246 0 : logtemp = log10(T)
247 : else
248 33684 : logtemp = logT
249 : end if
250 :
251 33684 : if (logtemp > 20) then
252 0 : info = -1
253 0 : return
254 : end if
255 :
256 33684 : if (Rho == arg_not_provided) then
257 0 : if (logRho == arg_not_provided) then
258 0 : info = -3
259 0 : return
260 : end if
261 0 : den = exp10(logRho)
262 : else
263 33684 : den = Rho
264 : end if
265 :
266 33684 : if (Rho <= 0) then
267 0 : info = -1
268 0 : return
269 : end if
270 :
271 33684 : if (logRho == arg_not_provided) then
272 0 : logden = log10(Rho)
273 : else
274 33684 : logden = logRho
275 : end if
276 :
277 33684 : if (logden > 20) then
278 0 : info = -1
279 0 : return
280 : end if
281 :
282 : !..initialize
283 33684 : spair = 0.0d0
284 33684 : spairdt = 0.0d0
285 33684 : spairdd = 0.0d0
286 33684 : spairda = 0.0d0
287 33684 : spairdz = 0.0d0
288 :
289 33684 : splas = 0.0d0
290 33684 : splasdt = 0.0d0
291 33684 : splasdd = 0.0d0
292 33684 : splasda = 0.0d0
293 33684 : splasdz = 0.0d0
294 :
295 33684 : sphot = 0.0d0
296 33684 : sphotdt = 0.0d0
297 33684 : sphotdd = 0.0d0
298 33684 : sphotda = 0.0d0
299 33684 : sphotdz = 0.0d0
300 :
301 33684 : sbrem = 0.0d0
302 33684 : sbremdt = 0.0d0
303 33684 : sbremdd = 0.0d0
304 33684 : sbremda = 0.0d0
305 33684 : sbremdz = 0.0d0
306 :
307 33684 : sreco = 0.0d0
308 33684 : srecodt = 0.0d0
309 33684 : srecodd = 0.0d0
310 33684 : srecoda = 0.0d0
311 33684 : srecodz = 0.0d0
312 :
313 33684 : snu = 0.0d0
314 33684 : dsnudt = 0.0d0
315 33684 : dsnudd = 0.0d0
316 33684 : dsnuda = 0.0d0
317 33684 : dsnudz = 0.0d0
318 :
319 :
320 33684 : call set_inputs(input,temp,logtemp,den,logden,abar,zbar)
321 :
322 : !..do the requested types
323 :
324 33684 : if (flags(pair_neu_type)) call pair_neu(spair,spairdt,spairdd,spairda,spairdz, input)
325 33684 : if (flags(plas_neu_type)) call plas_neu(splas,splasdt,splasdd,splasda,splasdz, input)
326 33684 : if (flags(phot_neu_type)) call phot_neu(sphot,sphotdt,sphotdd,sphotda,sphotdz, input)
327 33684 : if (flags(brem_neu_type)) call brem_neu(sbrem,sbremdt,sbremdd,sbremda,sbremdz, input)
328 33684 : if (flags(reco_neu_type)) call reco_neu(sreco,srecodt,srecodd,srecoda,srecodz, input)
329 :
330 : !..convert from erg/cm^3/s to erg/g/s
331 : !..comment these out to duplicate the itoh et al plots
332 :
333 33684 : spair = spair * input% deni
334 33684 : spairdt = spairdt * input% deni
335 33684 : spairdd = spairdd * input% deni - spair * input% deni
336 33684 : spairda = spairda * input% deni
337 33684 : spairdz = spairdz * input% deni
338 :
339 33684 : splas = splas * input% deni
340 33684 : splasdt = splasdt * input% deni
341 33684 : splasdd = splasdd * input% deni - splas * input% deni
342 33684 : splasda = splasda * input% deni
343 33684 : splasdz = splasdz * input% deni
344 :
345 33684 : sphot = sphot * input% deni
346 33684 : sphotdt = sphotdt * input% deni
347 33684 : sphotdd = sphotdd * input% deni - sphot * input% deni
348 33684 : sphotda = sphotda * input% deni
349 33684 : sphotdz = sphotdz * input% deni
350 :
351 33684 : sbrem = sbrem * input% deni
352 33684 : sbremdt = sbremdt * input% deni
353 33684 : sbremdd = sbremdd * input% deni - sbrem * input% deni
354 33684 : sbremda = sbremda * input% deni
355 33684 : sbremdz = sbremdz * input% deni
356 :
357 33684 : sreco = sreco * input% deni
358 33684 : srecodt = srecodt * input% deni
359 33684 : srecodd = srecodd * input% deni - sreco * input% deni
360 33684 : srecoda = srecoda * input% deni
361 33684 : srecodz = srecodz * input% deni
362 :
363 : !..calculate temperature cutoff factor
364 :
365 33684 : if (input% logtemp <= log10_Tlim .and. log10_Tlim > log10Tmin_neu) then
366 27684 : dtlim = log10_Tlim - log10Tmin_neu
367 : tcutoff_factor = 0.5d0* &
368 27684 : (1 - cospi((input% logtemp - log10Tmin_neu)/(log10_Tlim - log10Tmin_neu)))
369 :
370 :
371 : dtcutoff_factordt = 0.5d0 * pi * sinpi((input% logtemp - log10Tmin_neu)/dtlim) * &
372 27684 : 1.d0/(dtlim * temp * ln10)
373 :
374 27684 : splasdt = tcutoff_factor * splasdt + dtcutoff_factordt * splas
375 27684 : splasdd = tcutoff_factor * splasdd
376 27684 : splasda = tcutoff_factor * splasda
377 27684 : splasdz = tcutoff_factor * splasdz
378 27684 : splas = tcutoff_factor * splas
379 :
380 27684 : spairdt = tcutoff_factor * spairdt + dtcutoff_factordt * spair
381 27684 : spairdd = tcutoff_factor * spairdd
382 27684 : spairda = tcutoff_factor * spairda
383 27684 : spairdz = tcutoff_factor * spairdz
384 27684 : spair = tcutoff_factor * spair
385 :
386 27684 : sphotdt = tcutoff_factor * sphotdt + dtcutoff_factordt * sphot
387 27684 : sphotdd = tcutoff_factor * sphotdd
388 27684 : sphotda = tcutoff_factor * sphotda
389 27684 : sphotdz = tcutoff_factor * sphotdz
390 27684 : sphot = tcutoff_factor * sphot
391 :
392 27684 : sbremdt = tcutoff_factor * sbremdt + dtcutoff_factordt * sbrem
393 27684 : sbremdd = tcutoff_factor * sbremdd
394 27684 : sbremda = tcutoff_factor * sbremda
395 27684 : sbremdz = tcutoff_factor * sbremdz
396 27684 : sbrem = tcutoff_factor * sbrem
397 :
398 27684 : srecodt = tcutoff_factor * srecodt + dtcutoff_factordt * sreco
399 27684 : srecodd = tcutoff_factor * srecodd
400 27684 : srecoda = tcutoff_factor * srecoda
401 27684 : srecodz = tcutoff_factor * srecodz
402 27684 : sreco = tcutoff_factor * sreco
403 :
404 :
405 : end if
406 :
407 :
408 : !..the total neutrino loss rate
409 33684 : snu = splas + spair + sphot + sbrem + sreco
410 33684 : dsnudt = splasdt + spairdt + sphotdt + sbremdt + srecodt
411 33684 : dsnudd = splasdd + spairdd + sphotdd + sbremdd + srecodd
412 33684 : dsnuda = splasda + spairda + sphotda + sbremda + srecoda
413 33684 : dsnudz = splasdz + spairdz + sphotdz + sbremdz + srecodz
414 :
415 33684 : if (is_bad(snu)) then
416 0 : info = -1
417 0 : return
418 : end if
419 :
420 : !..packup the results
421 :
422 33684 : loss(1) = snu
423 33684 : loss(2) = dsnudt
424 33684 : loss(3) = dsnudd
425 33684 : loss(4) = dsnuda
426 33684 : loss(5) = dsnudz
427 :
428 33684 : call store(pair_neu_type, spair, spairdt, spairdd, spairda, spairdz)
429 33684 : call store(plas_neu_type, splas, splasdt, splasdd, splasda, splasdz)
430 33684 : call store(phot_neu_type, sphot, sphotdt, sphotdd, sphotda, sphotdz)
431 33684 : call store(brem_neu_type, sbrem, sbremdt, sbremdd, sbremda, sbremdz)
432 33684 : call store(reco_neu_type, sreco, srecodt, srecodd, srecoda, srecodz)
433 :
434 :
435 : contains
436 :
437 :
438 168420 : subroutine store(neu_type, s, sdt, sdd, sda, sdz)
439 : integer, intent(in) :: neu_type
440 : real(dp), intent(in) :: s, sdt, sdd, sda, sdz
441 168420 : sources(neu_type,1) = s
442 168420 : sources(neu_type,2) = sdt
443 168420 : sources(neu_type,3) = sdd
444 168420 : sources(neu_type,4) = sda
445 168420 : sources(neu_type,5) = sdz
446 168420 : end subroutine store
447 :
448 : end subroutine neutrinos
449 :
450 :
451 33684 : subroutine set_inputs(input,temp,logtemp,den,logden,abar,zbar)
452 : type(inputs),intent(out) :: input
453 :
454 : real(dp), intent(in) :: temp,logtemp,den,logden,abar, zbar
455 33684 : real(dp) ::a0,a1,a2
456 :
457 33684 : input% temp = temp
458 33684 : input% logtemp = logtemp
459 33684 : input% den = den
460 33684 : input% logden = logden
461 33684 : input% abar = abar
462 33684 : input% zbar = zbar
463 33684 : input% den6 = input% den * 1.0d-6
464 :
465 : !..to avoid lots of divisions
466 33684 : input% deni = 1.0d0 / input% den
467 33684 : input% tempi = 1.0d0 / input% temp
468 33684 : input% abari = 1.0d0 / input% abar
469 33684 : input% zbari = 1.0d0 / input% zbar
470 :
471 :
472 : !..some composition variables
473 33684 : input% ye = input% zbar * input% abari
474 : !xmue = abar * zbari
475 :
476 :
477 : !..some frequent factors
478 33684 : input% t9 = input% temp * 1.0d-9
479 33684 : input% xl = input% t9 * con1
480 33684 : input% xldt = 1.0d-9 * con1
481 33684 : input% xlp5 = sqrt(input% xl)
482 33684 : input% xl2 = input% xl * input% xl
483 33684 : input% xl3 = input% xl2 * input% xl
484 33684 : input% xl4 = input% xl3 * input% xl
485 33684 : input% xl5 = input% xl4 * input% xl
486 33684 : input% xl6 = input% xl5 * input% xl
487 33684 : input% xl7 = input% xl6 * input% xl
488 33684 : input% xl8 = input% xl7 * input% xl
489 33684 : input% xl9 = input% xl8 * input% xl
490 33684 : input% xlmp5 = 1.0d0 / input% xlp5
491 33684 : input% xlm1 = 1.0d0 / input% xl
492 33684 : input% xlm2 = input% xlm1*input% xlm1
493 33684 : input% xlm3 = input% xlm1*input% xlm2
494 33684 : input% xlm4 = input% xlm1*input% xlm3
495 :
496 :
497 33684 : input% rm = input% den*input% ye
498 33684 : input% rmdd = input% ye
499 33684 : input% rmda = -input% rm*input% abari
500 33684 : input% rmdz = input% den *input% abari
501 33684 : input% rmi = 1.0d0/input% rm
502 :
503 :
504 33684 : a0 = input% rm * 1.0d-9
505 33684 : a1 = pow(a0,one_third)
506 33684 : input% zeta = a1 * input% xlm1
507 33684 : input% zetadt = -a1 * input% xlm2 * input% xldt
508 33684 : a2 = one_third * a1*input%rmi * input% xlm1
509 33684 : input% zetadd = a2 * input%rmdd
510 33684 : input% zetada = a2 * input%rmda
511 33684 : input% zetadz = a2 * input%rmdz
512 :
513 33684 : input% zeta2 = input%zeta * input%zeta
514 33684 : input% zeta3 = input%zeta2 * input%zeta
515 :
516 :
517 33684 : end subroutine set_inputs
518 :
519 :
520 33684 : subroutine phot_neu(sphot,sphotdt,sphotdd,sphotda,sphotdz, input)
521 : real(dp), intent(out) :: sphot,sphotdt,sphotdd,sphotda,sphotdz
522 : type(inputs), intent(in) :: input
523 :
524 33684 : real(dp) :: tau,taudt,cos1,cos2,cos3,cos4,cos5,sin1,sin2, &
525 33684 : sin3,sin4,sin5,last,xast, &
526 33684 : fphot,fphotdt,fphotdd,fphotda,fphotdz, &
527 33684 : qphot,qphotdt,qphotdd,qphotda,qphotdz
528 :
529 33684 : real(dp) :: a0,a1,a2,a3,c00,c01,c02,c03,c04,c05,c06,cc, &
530 33684 : c10,c11,c12,c13,c14,c15,c16,c20,c21,c22,c23,c24, &
531 33684 : c25,c26,dd01,dd02,dd03,dd04,dd05,dd11,dd12, &
532 33684 : dd13,dd14,dd15,dd21,dd22,dd23,dd24,dd25,f0,f1,f2,z, &
533 33684 : dum,dumdt,dumdd,dumda,dumdz
534 :
535 33684 : real(dp) :: xnum,xnumdt,xnumdd,xnumda,xnumdz, &
536 33684 : xden,xdendt,xdendd,xdenda,xdendz
537 :
538 33684 : real(dp) :: dccdt
539 :
540 : !..photoneutrino process section
541 : !..for reactions like e- + gamma => e- + nu_e + nubar_e
542 : !.. e+ + gamma => e+ + nu_e + nubar_e
543 : !..equation 3.8 for tau, equation 3.6 for cc,
544 : !..and table 2 written out for speed
545 33684 : dccdt = 0d0
546 :
547 33684 : sphot = 0.0d0
548 33684 : sphotdt = 0.0d0
549 33684 : sphotdd = 0.0d0
550 33684 : sphotda = 0.0d0
551 33684 : sphotdz = 0.0d0
552 :
553 33684 : if(input% temp <=1.0d7) then
554 0 : return
555 33684 : else if (input% temp >= 1.0d7 .and. input% temp < 1.0d8) then
556 28684 : tau = input% logtemp - 7d0
557 28684 : cc = 0.5654d0 + tau
558 28684 : dccdt = 1d0/(ln10*input% temp)
559 28684 : c00 = 1.008d11
560 28684 : c01 = 0.0d0
561 28684 : c02 = 0.0d0
562 28684 : c03 = 0.0d0
563 28684 : c04 = 0.0d0
564 28684 : c05 = 0.0d0
565 28684 : c06 = 0.0d0
566 28684 : c10 = 8.156d10
567 28684 : c11 = 9.728d8
568 28684 : c12 = -3.806d9
569 28684 : c13 = -4.384d9
570 28684 : c14 = -5.774d9
571 28684 : c15 = -5.249d9
572 28684 : c16 = -5.153d9
573 28684 : c20 = 1.067d11
574 28684 : c21 = -9.782d9
575 28684 : c22 = -7.193d9
576 28684 : c23 = -6.936d9
577 28684 : c24 = -6.893d9
578 28684 : c25 = -7.041d9
579 28684 : c26 = -7.193d9
580 28684 : dd01 = 0.0d0
581 28684 : dd02 = 0.0d0
582 28684 : dd03 = 0.0d0
583 28684 : dd04 = 0.0d0
584 28684 : dd05 = 0.0d0
585 28684 : dd11 = -1.879d10
586 28684 : dd12 = -9.667d9
587 28684 : dd13 = -5.602d9
588 28684 : dd14 = -3.370d9
589 28684 : dd15 = -1.825d9
590 28684 : dd21 = -2.919d10
591 28684 : dd22 = -1.185d10
592 28684 : dd23 = -7.270d9
593 28684 : dd24 = -4.222d9
594 28684 : dd25 = -1.560d9
595 :
596 5000 : else if (input% temp >= 1.0d8 .and. input% temp < 1.0d9) then
597 2000 : tau = input% logtemp - 8d0
598 2000 : cc = 1.5654d0
599 2000 : dccdt = 0d0
600 2000 : c00 = 9.889d10
601 2000 : c01 = -4.524d8
602 2000 : c02 = -6.088d6
603 2000 : c03 = 4.269d7
604 2000 : c04 = 5.172d7
605 2000 : c05 = 4.910d7
606 2000 : c06 = 4.388d7
607 2000 : c10 = 1.813d11
608 2000 : c11 = -7.556d9
609 2000 : c12 = -3.304d9
610 2000 : c13 = -1.031d9
611 2000 : c14 = -1.764d9
612 2000 : c15 = -1.851d9
613 2000 : c16 = -1.928d9
614 2000 : c20 = 9.750d10
615 2000 : c21 = 3.484d10
616 2000 : c22 = 5.199d9
617 2000 : c23 = -1.695d9
618 2000 : c24 = -2.865d9
619 2000 : c25 = -3.395d9
620 2000 : c26 = -3.418d9
621 2000 : dd01 = -1.135d8
622 2000 : dd02 = 1.256d8
623 2000 : dd03 = 5.149d7
624 2000 : dd04 = 3.436d7
625 2000 : dd05 = 1.005d7
626 2000 : dd11 = 1.652d9
627 2000 : dd12 = -3.119d9
628 2000 : dd13 = -1.839d9
629 2000 : dd14 = -1.458d9
630 2000 : dd15 = -8.956d8
631 2000 : dd21 = -1.548d10
632 2000 : dd22 = -9.338d9
633 2000 : dd23 = -5.899d9
634 2000 : dd24 = -3.035d9
635 2000 : dd25 = -1.598d9
636 :
637 3000 : else if (input% temp >= 1.0d9) then
638 3000 : tau = input% logtemp - 9d0
639 3000 : cc = 1.5654d0
640 3000 : dccdt = 0d0
641 3000 : c00 = 9.581d10
642 3000 : c01 = 4.107d8
643 3000 : c02 = 2.305d8
644 3000 : c03 = 2.236d8
645 3000 : c04 = 1.580d8
646 3000 : c05 = 2.165d8
647 3000 : c06 = 1.721d8
648 3000 : c10 = 1.459d12
649 3000 : c11 = 1.314d11
650 3000 : c12 = -1.169d11
651 3000 : c13 = -1.765d11
652 3000 : c14 = -1.867d11
653 3000 : c15 = -1.983d11
654 3000 : c16 = -1.896d11
655 3000 : c20 = 2.424d11
656 3000 : c21 = -3.669d9
657 3000 : c22 = -8.691d9
658 3000 : c23 = -7.967d9
659 3000 : c24 = -7.932d9
660 3000 : c25 = -7.987d9
661 3000 : c26 = -8.333d9
662 3000 : dd01 = 4.724d8
663 3000 : dd02 = 2.976d8
664 3000 : dd03 = 2.242d8
665 3000 : dd04 = 7.937d7
666 3000 : dd05 = 4.859d7
667 3000 : dd11 = -7.094d11
668 3000 : dd12 = -3.697d11
669 3000 : dd13 = -2.189d11
670 3000 : dd14 = -1.273d11
671 3000 : dd15 = -5.705d10
672 3000 : dd21 = -2.254d10
673 3000 : dd22 = -1.551d10
674 3000 : dd23 = -7.793d9
675 3000 : dd24 = -4.489d9
676 3000 : dd25 = -2.185d9
677 : end if
678 :
679 33684 : taudt = iln10*input% tempi
680 :
681 :
682 : !..equation 3.7, compute the expensive trig functions only one time
683 33684 : cos1 = cos(fac1*tau)
684 33684 : sin1 = sin(fac1*tau)
685 :
686 : ! double, triple, etc. angle formulas
687 : ! sin/cos (2 fac1 tau)
688 33684 : sin2 = 2.0d0 * sin1 * cos1
689 33684 : cos2 = 2.0d0 * cos1 * cos1 - 1.0d0
690 :
691 : ! sin/cos (3 fac1 tau)
692 33684 : sin3 = sin1 * (3.0d0 - 4.0d0 * sin1 * sin1)
693 33684 : cos3 = cos1 * (4.0d0 * cos1 * cos1 - 3.0d0)
694 :
695 : ! sin/cos (4 fac1 tau) -- use double angle on sin2/cos2
696 33684 : sin4 = 2.0d0 * sin2 * cos2
697 33684 : cos4 = 2.0d0 * cos2 * cos2 - 1.0d0
698 :
699 : ! sin/cos (5 fac1 tau)
700 33684 : sin5 = sin1 * (5.0d0 - sin1 * sin1 * (20.0d0 - 16.0d0 * sin1 * sin1))
701 33684 : cos5 = cos1 * (cos1 * cos1 * (16.0d0 * cos1 * cos1 - 20.0d0) + 5.0d0)
702 :
703 33684 : last = cos(fac2*tau)
704 33684 : xast = sin(fac2*tau)
705 :
706 : a0 = 0.5d0*c00 &
707 : + c01*cos1 + dd01*sin1 + c02*cos2 + dd02*sin2 &
708 : + c03*cos3 + dd03*sin3 + c04*cos4 + dd04*sin4 &
709 33684 : + c05*cos5 + dd05*sin5 + 0.5d0*c06*last
710 :
711 : f0 = taudt*fac1*(-c01*sin1 + dd01*cos1 - c02*sin2*2.0d0 &
712 : + dd02*cos2*2.0d0 - c03*sin3*3.0d0 + dd03*cos3*3.0d0 &
713 : - c04*sin4*4.0d0 + dd04*cos4*4.0d0 &
714 : - c05*sin5*5.0d0 + dd05*cos5*5.0d0) &
715 33684 : - 0.5d0*c06*xast*fac2*taudt
716 :
717 : a1 = 0.5d0*c10 &
718 : + c11*cos1 + dd11*sin1 + c12*cos2 + dd12*sin2 &
719 : + c13*cos3 + dd13*sin3 + c14*cos4 + dd14*sin4 &
720 33684 : + c15*cos5 + dd15*sin5 + 0.5d0*c16*last
721 :
722 : f1 = taudt*fac1*(-c11*sin1 + dd11*cos1 - c12*sin2*2.0d0 &
723 : + dd12*cos2*2.0d0 - c13*sin3*3.0d0 + dd13*cos3*3.0d0 &
724 : - c14*sin4*4.0d0 + dd14*cos4*4.0d0 - c15*sin5*5.0d0 &
725 33684 : + dd15*cos5*5.0d0) - 0.5d0*c16*xast*fac2*taudt
726 :
727 : a2 = 0.5d0*c20 &
728 : + c21*cos1 + dd21*sin1 + c22*cos2 + dd22*sin2 &
729 : + c23*cos3 + dd23*sin3 + c24*cos4 + dd24*sin4 &
730 33684 : + c25*cos5 + dd25*sin5 + 0.5d0*c26*last
731 :
732 : f2 = taudt*fac1*(-c21*sin1 + dd21*cos1 - c22*sin2*2.0d0 &
733 : + dd22*cos2*2.0d0 - c23*sin3*3.0d0 + dd23*cos3*3.0d0 &
734 : - c24*sin4*4.0d0 + dd24*cos4*4.0d0 - c25*sin5*5.0d0 &
735 33684 : + dd25*cos5*5.0d0) - 0.5d0*c26*xast*fac2*taudt
736 :
737 : !..equation 3.4
738 33684 : dum = a0 + a1*input% zeta + a2*input% zeta2
739 33684 : dumdt = f0 + f1*input% zeta + a1*input% zetadt + f2*input% zeta2 + 2.0d0*a2*input% zeta*input% zetadt
740 33684 : dumdd = a1*input% zetadd + 2.0d0*a2*input% zeta*input% zetadd
741 33684 : dumda = a1*input% zetada + 2.0d0*a2*input% zeta*input% zetada
742 33684 : dumdz = a1*input% zetadz + 2.0d0*a2*input% zeta*input% zetadz
743 :
744 33684 : z = exp(-cc*input% zeta)
745 :
746 33684 : xnum = dum*z
747 33684 : xnumdt = dumdt*z - dum*z*(dccdt*input% zeta + input% zetadt*cc)
748 33684 : xnumdd = dumdd*z - dum*z*cc*input% zetadd
749 33684 : xnumda = dumda*z - dum*z*cc*input% zetada
750 33684 : xnumdz = dumdz*z - dum*z*cc*input% zetadz
751 :
752 33684 : xden = input% zeta3 + 6.290d-3*input% xlm1 + 7.483d-3*input% xlm2 + 3.061d-4*input% xlm3
753 :
754 33684 : dum = 3.0d0*input% zeta2
755 : xdendt = dum*input% zetadt - input% xldt*(6.290d-3*input% xlm2 &
756 33684 : + 2.0d0*7.483d-3*input% xlm3 + 3.0d0*3.061d-4*input% xlm4)
757 33684 : xdendd = dum*input% zetadd
758 33684 : xdenda = dum*input% zetada
759 33684 : xdendz = dum*input% zetadz
760 :
761 33684 : dum = 1.0d0/xden
762 33684 : fphot = xnum*dum
763 33684 : fphotdt = (xnumdt - fphot*xdendt)*dum
764 33684 : fphotdd = (xnumdd - fphot*xdendd)*dum
765 33684 : fphotda = (xnumda - fphot*xdenda)*dum
766 33684 : fphotdz = (xnumdz - fphot*xdendz)*dum
767 :
768 :
769 : !..equation 3.3
770 33684 : a0 = 1.0d0 + 2.045d0 * input% xl
771 33684 : xnum = 0.666d0*pow(a0,-2.066d0)
772 33684 : xnumdt = -2.066d0*xnum/a0 * 2.045d0*input% xldt
773 :
774 33684 : dum = 1.875d8*input% xl + 1.653d8*input% xl2 + 8.499d8*input% xl3 - 1.604d8*input% xl4
775 : dumdt = input% xldt*(1.875d8 + 2.0d0*1.653d8*input% xl + 3.0d0*8.499d8*input% xl2 &
776 33684 : - 4.0d0*1.604d8*input% xl3)
777 :
778 33684 : z = 1.0d0/dum
779 33684 : xden = 1.0d0 + input% rm*z
780 33684 : xdendt = -input% rm*z*z*dumdt
781 33684 : xdendd = input% rmdd*z
782 33684 : xdenda = input% rmda*z
783 33684 : xdendz = input% rmdz*z
784 :
785 33684 : z = 1.0d0/xden
786 33684 : qphot = xnum*z
787 33684 : qphotdt = (xnumdt - qphot*xdendt)*z
788 33684 : dum = -qphot*z
789 33684 : qphotdd = dum*xdendd
790 33684 : qphotda = dum*xdenda
791 33684 : qphotdz = dum*xdendz
792 :
793 : !..equation 3.2
794 33684 : sphot = input% xl5 * fphot
795 33684 : sphotdt = 5.0d0*input% xl4*input% xldt*fphot + input% xl5*fphotdt
796 33684 : sphotdd = input% xl5*fphotdd
797 33684 : sphotda = input% xl5*fphotda
798 33684 : sphotdz = input% xl5*fphotdz
799 :
800 33684 : a1 = sphot
801 33684 : sphot = input% rm*a1
802 33684 : sphotdt = input% rm*sphotdt
803 33684 : sphotdd = input% rm*sphotdd + input% rmdd*a1
804 33684 : sphotda = input% rm*sphotda + input% rmda*a1
805 33684 : sphotdz = input% rm*sphotdz + input% rmdz*a1
806 :
807 33684 : a1 = tfac4*(1.0d0 - tfac3 * qphot)
808 33684 : a2 = -tfac4*tfac3
809 :
810 33684 : a3 = sphot
811 33684 : sphot = a1*a3
812 33684 : sphotdt = a1*sphotdt + a2*qphotdt*a3
813 33684 : sphotdd = a1*sphotdd + a2*qphotdd*a3
814 33684 : sphotda = a1*sphotda + a2*qphotda*a3
815 33684 : sphotdz = a1*sphotdz + a2*qphotdz*a3
816 :
817 :
818 : if (.false.) then
819 : write(*,*) 'logT', input% logtemp
820 : write(*,*) 'logRho', input% logden
821 : write(*,*) 'sphot', sphot
822 : write(*,*)
823 : end if
824 :
825 33684 : if (sphot <= 0.0d0) then
826 24 : sphot = 0.0d0
827 24 : sphotdt = 0.0d0
828 24 : sphotdd = 0.0d0
829 24 : sphotda = 0.0d0
830 24 : sphotdz = 0.0d0
831 : end if
832 :
833 : end subroutine phot_neu
834 :
835 :
836 30182 : subroutine brem_neu_weak_degen(sbrem,sbremdt,sbremdd,sbremda,sbremdz,t8, input)
837 : type(t8s), intent(in) :: t8
838 : real(dp), intent(out) :: sbrem,sbremdt,sbremdd,sbremda,sbremdz
839 : type(inputs), intent(in) :: input
840 :
841 30182 : real(dp) :: a0,c00,c01,c02,c03,c04, &
842 30182 : dd00,dd01,dd02,&
843 30182 : f0,z, &
844 30182 : dum,dumdt,dumdd,dumda,dumdz
845 :
846 30182 : real(dp) :: xnum,xnumdt,xnumdd,xnumda,xnumdz, &
847 30182 : xden,xdendt,xdendd,xdenda,xdendz
848 :
849 : ! brem
850 :
851 30182 : real(dp) :: eta,etadt,etadd,etada,etadz,etam1,etam2,etam3, &
852 30182 : fbrem,fbremdt,fbremdd,fbremda,fbremdz, &
853 30182 : gbrem,gbremdt,gbremdd,gbremda,gbremdz
854 :
855 30182 : real(dp) :: p
856 :
857 30182 : sbrem=0d0
858 30182 : sbremdt=0d0
859 30182 : sbremdd=0d0
860 30182 : sbremda=0d0
861 30182 : sbremdz=0d0
862 :
863 : !..equation 5.3
864 30182 : dum = 7.05d6 * t8% t832 + 5.12d4 * t8% t83
865 30182 : dumdt = (1.5d0*7.05d6*t8% t812 + 3.0d0*5.12d4*t8% t82)*1.0d-8
866 :
867 30182 : z = 1.0d0/dum
868 30182 : eta = input% rm*z
869 30182 : etadt = -input% rm*z*z*dumdt
870 30182 : etadd = input% rmdd*z
871 30182 : etada = input% rmda*z
872 30182 : etadz = input% rmdz*z
873 :
874 30182 : etam1 = 1.0d0/eta
875 30182 : etam2 = etam1 * etam1
876 30182 : etam3 = etam2 * etam1
877 :
878 :
879 : !..equation 5.2
880 30182 : a0 = 23.5d0 + 6.83d4*t8% t8m2 + 7.81d8*t8% t8m5
881 30182 : f0 = (-2.0d0*6.83d4*t8% t8m3 - 5.0d0*7.81d8*t8% t8m6)*1.0d-8
882 30182 : xnum = 1.0d0/a0
883 :
884 30182 : dum = 1.0d0 + 1.47d0*etam1 + 3.29d-2*etam2
885 30182 : z = -1.47d0*etam2 - 2.0d0*3.29d-2*etam3
886 30182 : dumdt = z*etadt
887 30182 : dumdd = z*etadd
888 30182 : dumda = z*etada
889 30182 : dumdz = z*etadz
890 :
891 30182 : c00 = 1.26d0 * (1.0d0+etam1)
892 30182 : z = -1.26d0*etam2
893 30182 : c01 = z*etadt
894 30182 : c02 = z*etadd
895 30182 : c03 = z*etada
896 30182 : c04 = z*etadz
897 :
898 30182 : z = 1.0d0/dum
899 30182 : xden = c00*z
900 30182 : xdendt = (c01 - xden*dumdt)*z
901 30182 : xdendd = (c02 - xden*dumdd)*z
902 30182 : xdenda = (c03 - xden*dumda)*z
903 30182 : xdendz = (c04 - xden*dumdz)*z
904 :
905 30182 : fbrem = xnum + xden
906 30182 : fbremdt = -xnum*xnum*f0 + xdendt
907 30182 : fbremdd = xdendd
908 30182 : fbremda = xdenda
909 30182 : fbremdz = xdendz
910 :
911 :
912 : !..equation 5.9
913 30182 : a0 = 230.0d0 + 6.7d5*t8% t8m2 + 7.66d9*t8% t8m5
914 30182 : f0 = (-2.0d0*6.7d5*t8% t8m3 - 5.0d0*7.66d9*t8% t8m6)*1.0d-8
915 :
916 30182 : z = 1.0d0 + input% rm*1.0d-9
917 30182 : dum = a0*z
918 30182 : dumdt = f0*z
919 30182 : z = a0*1.0d-9
920 30182 : dumdd = z*input% rmdd
921 30182 : dumda = z*input% rmda
922 30182 : dumdz = z*input% rmdz
923 :
924 30182 : xnum = 1.0d0/dum
925 30182 : z = -xnum*xnum
926 30182 : xnumdt = z*dumdt
927 30182 : xnumdd = z*dumdd
928 30182 : xnumda = z*dumda
929 30182 : xnumdz = z*dumdz
930 :
931 30182 : p = pow(t8% t8,3.85d0)
932 30182 : c00 = 7.75d5*t8% t832 + 247.0d0*p
933 30182 : dd00 = (1.5d0*7.75d5*t8% t812 + 3.85d0*247.0d0*p/t8% T8)*1.0d-8
934 :
935 30182 : p = pow(t8% t8,1.4d0)
936 30182 : c01 = 4.07d0 + 0.0240d0 * p
937 30182 : dd01 = 1.4d0*0.0240d0*(p/t8% T8)*1.0d-8
938 :
939 30182 : p = pow(t8% t8,-0.110d0)
940 30182 : c02 = 4.59d-5 * p
941 30182 : dd02 = -0.11d0*4.59d-5*(p/t8% T8)*1.0d-8
942 :
943 30182 : z = pow(input% den,0.656d0)
944 30182 : dum = c00*input% rmi + c01 + c02*z
945 30182 : dumdt = dd00*input% rmi + dd01 + dd02*z
946 30182 : z = -c00*input% rmi*input% rmi
947 30182 : dumdd = z*input% rmdd + 0.656d0*c02*pow(input% den,-0.344d0)
948 30182 : dumda = z*input% rmda
949 30182 : dumdz = z*input% rmdz
950 :
951 30182 : xden = 1.0d0/dum
952 30182 : z = -xden*xden
953 30182 : xdendt = z*dumdt
954 30182 : xdendd = z*dumdd
955 30182 : xdenda = z*dumda
956 30182 : xdendz = z*dumdz
957 :
958 :
959 30182 : gbrem = xnum + xden
960 30182 : gbremdt = xnumdt + xdendt
961 30182 : gbremdd = xnumdd + xdendd
962 30182 : gbremda = xnumda + xdenda
963 30182 : gbremdz = xnumdz + xdendz
964 :
965 :
966 : !..equation 5.1
967 30182 : dum = 0.5738d0*input% zbar*input% ye*t8% t86*input% den
968 30182 : dumdt = 0.5738d0*input% zbar*input% ye*6.0d0*t8% t85*input% den*1.0d-8
969 30182 : dumdd = 0.5738d0*input% zbar*input% ye*t8% t86
970 30182 : dumda = -dum*input% abari
971 30182 : dumdz = 0.5738d0*2.0d0*input% ye*t8% t86*input% den
972 :
973 30182 : z = tfac4*fbrem - tfac5*gbrem
974 30182 : sbrem = dum * z
975 30182 : sbremdt = dumdt*z + dum*(tfac4*fbremdt - tfac5*gbremdt)
976 30182 : sbremdd = dumdd*z + dum*(tfac4*fbremdd - tfac5*gbremdd)
977 30182 : sbremda = dumda*z + dum*(tfac4*fbremda - tfac5*gbremda)
978 30182 : sbremdz = dumdz*z + dum*(tfac4*fbremdz - tfac5*gbremdz)
979 :
980 :
981 30182 : end subroutine brem_neu_weak_degen
982 :
983 :
984 4905 : subroutine brem_neu_liquid_metal(sbrem,sbremdt,sbremdd,sbremda,sbremdz,t8, input)
985 : type(t8s), intent(in) :: t8
986 : real(dp), intent(out) :: sbrem,sbremdt,sbremdd,sbremda,sbremdz
987 : type(inputs), intent(in) :: input
988 :
989 4905 : real(dp) :: a0,a1,c00,c01,c02,c03, &
990 4905 : z, &
991 4905 : dum,dumdt,dumdd,dumda,dumdz
992 :
993 : ! brem
994 :
995 4905 : real(dp) :: u,gm1,gm2,gm13,gm23,gm43,gm53,v,w,fb,gt,gb, &
996 4905 : fliq,fliqdt,fliqdd,fliqda,fliqdz, &
997 4905 : gliq,gliqdt,gliqdd,gliqda,gliqdz
998 :
999 4905 : real(dp) :: cos1,cos2,cos3,cos4,cos5,sin1,sin2, &
1000 4905 : sin3,sin4,sin5
1001 :
1002 4905 : real(dp) :: ft
1003 :
1004 4905 : sbrem=0d0
1005 4905 : sbremdt=0d0
1006 4905 : sbremdd=0d0
1007 4905 : sbremda=0d0
1008 4905 : sbremdz=0d0
1009 :
1010 :
1011 : !..liquid metal with c12 parameters (not too different for other elements)
1012 : !..equation 5.18 and 5.16
1013 4905 : u = fac3 * (input% logden - 3.0d0)
1014 4905 : a0 = iln10*fac3 * input% deni
1015 :
1016 : !..compute the expensive trig functions of equation 5.21 only once
1017 4905 : cos1 = cos(u)
1018 4905 : sin1 = sin(u)
1019 :
1020 : ! double, triple, etc. angle formulas
1021 : ! sin/cos (2 u)
1022 4905 : sin2 = 2.0d0 * sin1 * cos1
1023 4905 : cos2 = 2.0d0 * cos1 * cos1 - 1.0d0
1024 :
1025 : ! sin/cos (3 u)
1026 4905 : sin3 = sin1 * (3.0d0 - 4.0d0 * sin1 * sin1)
1027 4905 : cos3 = cos1 * (4.0d0 * cos1 * cos1 - 3.0d0)
1028 :
1029 : ! sin/cos (4 u) -- use double angle on sin2/cos2
1030 4905 : sin4 = 2.0d0 * sin2 * cos2
1031 4905 : cos4 = 2.0d0 * cos2 * cos2 - 1.0d0
1032 :
1033 : ! sin/cos (5 u)
1034 4905 : sin5 = sin1 * (5.0d0 - sin1 * sin1 * (20.0d0 - 16.0d0 * sin1 * sin1))
1035 4905 : cos5 = cos1 * (cos1 * cos1 * (16.0d0 * cos1 * cos1 - 20.0d0) + 5.0d0)
1036 :
1037 :
1038 : !..equation 5.21
1039 : fb = 0.5d0 * 0.17946d0 + 0.00945d0*u + 0.34529d0 &
1040 : - 0.05821d0*cos1 - 0.04969d0*sin1 &
1041 : - 0.01089d0*cos2 - 0.01584d0*sin2 &
1042 : - 0.01147d0*cos3 - 0.00504d0*sin3 &
1043 : - 0.00656d0*cos4 - 0.00281d0*sin4 &
1044 4905 : - 0.00519d0*cos5
1045 :
1046 : c00 = a0*(0.00945d0 &
1047 : + 0.05821d0*sin1 - 0.04969d0*cos1 &
1048 : + 0.01089d0*sin2*2.0d0 - 0.01584d0*cos2*2.0d0 &
1049 : + 0.01147d0*sin3*3.0d0 - 0.00504d0*cos3*3.0d0 &
1050 : + 0.00656d0*sin4*4.0d0 - 0.00281d0*cos4*4.0d0 &
1051 4905 : + 0.00519d0*sin5*5.0d0)
1052 :
1053 :
1054 : !..equation 5.22
1055 : ft = 0.5d0 * 0.06781d0 - 0.02342d0*u + 0.24819d0 &
1056 : - 0.00944d0*cos1 - 0.02213d0*sin1 &
1057 : - 0.01289d0*cos2 - 0.01136d0*sin2 &
1058 : - 0.00589d0*cos3 - 0.00467d0*sin3 &
1059 : - 0.00404d0*cos4 - 0.00131d0*sin4 &
1060 4905 : - 0.00330d0*cos5
1061 :
1062 : c01 = a0*(-0.02342d0 &
1063 : + 0.00944d0*sin1 - 0.02213d0*cos1 &
1064 : + 0.01289d0*sin2*2.0d0 - 0.01136d0*cos2*2.0d0 &
1065 : + 0.00589d0*sin3*3.0d0 - 0.00467d0*cos3*3.0d0 &
1066 : + 0.00404d0*sin4*4.0d0 - 0.00131d0*cos4*4.0d0 &
1067 4905 : + 0.00330d0*sin5*5.0d0)
1068 :
1069 :
1070 : !..equation 5.23
1071 : gb = 0.5d0 * 0.00766d0 - 0.01259d0*u + 0.07917d0 &
1072 : - 0.00710d0*cos1 + 0.02300d0*sin1 &
1073 : - 0.00028d0*cos2 - 0.01078d0*sin2 &
1074 : + 0.00232d0*cos3 + 0.00118d0*sin3 &
1075 : + 0.00044d0*cos4 - 0.00089d0*sin4 &
1076 4905 : + 0.00158d0*cos5
1077 :
1078 : c02 = a0*(-0.01259d0 &
1079 : + 0.00710d0*sin1 + 0.02300d0*cos1 &
1080 : + 0.00028d0*sin2*2.0d0 - 0.01078d0*cos2*2.0d0 &
1081 : - 0.00232d0*sin3*3.0d0 + 0.00118d0*cos3*3.0d0 &
1082 : - 0.00044d0*sin4*4.0d0 - 0.00089d0*cos4*4.0d0 &
1083 4905 : - 0.00158d0*sin5*5.0d0)
1084 :
1085 :
1086 : !..equation 5.24
1087 : gt = -0.5d0 * 0.00769d0 - 0.00829d0*u + 0.05211d0 &
1088 : + 0.00356d0*cos1 + 0.01052d0*sin1 &
1089 : - 0.00184d0*cos2 - 0.00354d0*sin2 &
1090 : + 0.00146d0*cos3 - 0.00014d0*sin3 &
1091 : + 0.00031d0*cos4 - 0.00018d0*sin4 &
1092 4905 : + 0.00069d0*cos5
1093 :
1094 : c03 = a0*(-0.00829d0 &
1095 : - 0.00356d0*sin1 + 0.01052d0*cos1 &
1096 : + 0.00184d0*sin2*2.0d0 - 0.00354d0*cos2*2.0d0 &
1097 : - 0.00146d0*sin3*3.0d0 - 0.00014d0*cos3*3.0d0 &
1098 : - 0.00031d0*sin4*4.0d0 - 0.00018d0*cos4*4.0d0 &
1099 4905 : - 0.00069d0*sin5*5.0d0)
1100 :
1101 :
1102 4905 : dum = 2.275d-1 * input% zbar * input% zbar*t8% t8m1 * pow(input% den6*input% abari, one_third)
1103 4905 : dumdt = -dum*input% tempi
1104 4905 : dumdd = one_third*dum * input% deni
1105 4905 : dumda = -one_third*dum*input% abari
1106 4905 : dumdz = 2.0d0*dum*input% zbari
1107 :
1108 4905 : gm1 = 1.0d0/dum
1109 4905 : gm2 = gm1*gm1
1110 4905 : gm13 = pow(gm1,one_third)
1111 4905 : gm23 = gm13 * gm13
1112 4905 : gm43 = gm13*gm1
1113 4905 : gm53 = gm23*gm1
1114 :
1115 :
1116 : !..equation 5.25 and 5.26
1117 4905 : v = -0.05483d0 - 0.01946d0*gm13 + 1.86310d0*gm23 - 0.78873d0*gm1
1118 4905 : a0 = one_third*0.01946d0*gm43 - two_thirds*1.86310d0*gm53 + 0.78873d0*gm2
1119 :
1120 4905 : w = -0.06711d0 + 0.06859d0*gm13 + 1.74360d0*gm23 - 0.74498d0*gm1
1121 4905 : a1 = -one_third*0.06859d0*gm43 - two_thirds*1.74360d0*gm53 + 0.74498d0*gm2
1122 :
1123 :
1124 : !..equation 5.19 and 5.20
1125 4905 : fliq = v*fb + (1.0d0 - v)*ft
1126 4905 : fliqdt = a0*dumdt*(fb - ft)
1127 4905 : fliqdd = a0*dumdd*(fb - ft) + v*c00 + (1.0d0 - v)*c01
1128 4905 : fliqda = a0*dumda*(fb - ft)
1129 4905 : fliqdz = a0*dumdz*(fb - ft)
1130 :
1131 4905 : gliq = w*gb + (1.0d0 - w)*gt
1132 4905 : gliqdt = a1*dumdt*(gb - gt)
1133 4905 : gliqdd = a1*dumdd*(gb - gt) + w*c02 + (1.0d0 - w)*c03
1134 4905 : gliqda = a1*dumda*(gb - gt)
1135 4905 : gliqdz = a1*dumdz*(gb - gt)
1136 :
1137 :
1138 : !..equation 5.17
1139 4905 : dum = 0.5738d0*input% zbar*input% ye*t8% t86*input% den
1140 4905 : dumdt = 0.5738d0*input% zbar*input% ye*6.0d0*t8% t85*input% den*1.0d-8
1141 4905 : dumdd = 0.5738d0*input% zbar*input% ye*t8% t86
1142 4905 : dumda = -dum*input% abari
1143 4905 : dumdz = 0.5738d0*2.0d0*input% ye*t8% t86*input% den
1144 :
1145 4905 : z = tfac4*fliq - tfac5*gliq
1146 4905 : sbrem = dum * z
1147 4905 : sbremdt = dumdt*z + dum*(tfac4*fliqdt - tfac5*gliqdt)
1148 4905 : sbremdd = dumdd*z + dum*(tfac4*fliqdd - tfac5*gliqdd)
1149 4905 : sbremda = dumda*z + dum*(tfac4*fliqda - tfac5*gliqda)
1150 4905 : sbremdz = dumdz*z + dum*(tfac4*fliqdz - tfac5*gliqdz)
1151 :
1152 4905 : end subroutine brem_neu_liquid_metal
1153 :
1154 :
1155 33684 : subroutine brem_neu(sbrem,sbremdt,sbremdd,sbremda,sbremdz, input)
1156 : real(dp), intent(out) :: sbrem,sbremdt,sbremdd,sbremda,sbremdz
1157 : type(inputs), intent(in) :: input
1158 :
1159 33684 : real(dp) ::tfermi
1160 :
1161 : ! brem
1162 :
1163 : type(t8s) :: t8
1164 :
1165 :
1166 : !..bremsstrahlung neutrino section
1167 : !..for reactions like e- + (z,a) => e- + (z,a) + nu + nubar
1168 : !.. n + n => n + n + nu + nubar
1169 : !.. n + p => n + p + nu + nubar
1170 :
1171 33684 : real(dp) :: alfa, beta, sb, sbdt, sbdd, sbda, sbdz
1172 33684 : real(dp) :: sb2, sbdt2, sbdd2, sbda2, sbdz2
1173 : real(dp), parameter :: tflo = 0.05d0, tfhi = 0.3d0
1174 33684 : real(dp) :: dtfracdt, dtfracdd, dtfracda, dtfracdz
1175 33684 : real(dp) :: dalfadt, dalfadd, dalfada, dalfadz
1176 33684 : real(dp) :: dbetadt, dbetadd, dbetada, dbetadz
1177 33684 : real(dp) :: dtfermidd, dtfermida, dtfermidz
1178 33684 : real(dp) :: A, B, C, D, U, dtfermidu, dtfracdtfermi, dtf, tfrac
1179 33684 : real(dp) :: dudv, dudd, duda, dudz
1180 :
1181 33684 : sbrem=0d0
1182 33684 : sbremdt=0d0
1183 33684 : sbremdd=0d0
1184 33684 : sbremda=0d0
1185 33684 : sbremdz=0d0
1186 :
1187 : !..equation 4.3
1188 :
1189 33684 : t8% t8 = input% temp * 1.0d-8
1190 33684 : t8% t812 = sqrt(t8% t8)
1191 33684 : t8% t832 = t8% t8 * t8% t812
1192 33684 : t8% t82 = t8% t8*t8% t8
1193 33684 : t8% t83 = t8% t82*t8% t8
1194 33684 : t8% t85 = t8% t82*t8% t83
1195 33684 : t8% t86 = t8% t85*t8% t8
1196 33684 : t8% t8m1 = 1.0d0/t8% t8
1197 33684 : t8% t8m2 = t8% t8m1*t8% t8m1
1198 33684 : t8% t8m3 = t8% t8m2*t8% t8m1
1199 33684 : t8% t8m5 = t8% t8m3*t8% t8m2
1200 33684 : t8% t8m6 = t8% t8m5*t8% t8m1
1201 :
1202 :
1203 33684 : A = 5.9302d9
1204 33684 : B = 1.d0
1205 33684 : C = 1.018d0
1206 33684 : D = 1.0d0
1207 :
1208 33684 : U = pow(input% den6*input% ye,two_thirds)
1209 33684 : tfermi = A * (sqrt(U) - D)
1210 :
1211 33684 : if (input% temp >= tfhi * tfermi) then
1212 :
1213 28779 : call brem_neu_weak_degen(sbrem,sbremdt,sbremdd,sbremda,sbremdz,t8, input)
1214 :
1215 4905 : else if (input% temp <= tflo * tfermi) then
1216 :
1217 3502 : call brem_neu_liquid_metal(sbrem,sbremdt,sbremdd,sbremda,sbremdz,t8, input)
1218 :
1219 : else ! blend
1220 :
1221 1403 : call brem_neu_weak_degen(sbrem,sbremdt,sbremdd,sbremda,sbremdz,t8, input)
1222 1403 : sb = sbrem
1223 1403 : sbdt = sbremdt
1224 1403 : sbdd = sbremdd
1225 1403 : sbda = sbremda
1226 1403 : sbdz = sbremdz
1227 :
1228 1403 : call brem_neu_liquid_metal(sbrem,sbremdt,sbremdd,sbremda,sbremdz,t8, input)
1229 1403 : sb2 = sbrem
1230 1403 : sbdt2 = sbremdt
1231 1403 : sbdd2 = sbremdd
1232 1403 : sbda2 = sbremda
1233 1403 : sbdz2 = sbremdz
1234 :
1235 1403 : dtf = tfhi - tflo
1236 1403 : tfrac = (input% temp / tfermi - tflo) / dtf
1237 1403 : alfa = 0.5d0 * (1d0 - cospi(tfrac))
1238 1403 : beta = 1d0 - alfa
1239 :
1240 :
1241 1403 : dtfermidu = (1d0/2d0) * A * pow(U,-1d0/2d0)
1242 1403 : dtfracdtfermi = -input% temp/(tfermi * tfermi * dtf )
1243 :
1244 : ! v = den6* ye = den *10**-6 * ye
1245 1403 : dudv = two_thirds * pow(input% den6 * input% ye,-1d0/3d0)
1246 :
1247 1403 : dudd = dudv * 1d-6 * input% ye
1248 1403 : duda = dudv * input% den6 * input% ye * input% abari * (-1d0)
1249 1403 : dudz = dudv * input% den6 * input% abari
1250 :
1251 :
1252 1403 : dtfermidd = dtfermidu * dudd
1253 1403 : dtfermida = dtfermidu * duda
1254 1403 : dtfermidz = dtfermidu * dudz
1255 :
1256 1403 : dtfracdt = 1.0d0/(tfermi * dtf)
1257 :
1258 1403 : dtfracdd = dtfermidd * dtfracdtfermi
1259 1403 : dtfracda = dtfermida * dtfracdtfermi
1260 1403 : dtfracdz = dtfermidz * dtfracdtfermi
1261 :
1262 1403 : dalfadt = dtfracdt * 0.5d0 * pi * sinpi(tfrac)
1263 1403 : dalfadd = dtfracdd * 0.5d0 * pi * sinpi(tfrac)
1264 1403 : dalfada = dtfracda * 0.5d0 * pi * sinpi(tfrac)
1265 1403 : dalfadz = dtfracdz * 0.5d0 * pi * sinpi(tfrac)
1266 :
1267 1403 : dbetadt = -dalfadt
1268 1403 : dbetadd = -dalfadd
1269 1403 : dbetada = -dalfada
1270 1403 : dbetadz = -dalfadz
1271 :
1272 1403 : sbrem = alfa * sb + beta * sb2
1273 1403 : sbremdt = alfa * sbdt + beta * sbdt2 + dalfadt * sb + dbetadt * sb2
1274 1403 : sbremdd = alfa * sbdd + beta * sbdd2 + dalfadd * sb + dbetadd * sb2
1275 1403 : sbremda = alfa * sbda + beta * sbda2 + dalfada * sb + dbetada * sb2
1276 1403 : sbremdz = alfa * sbdz + beta * sbdz2 + dalfadz * sb + dbetadz * sb2
1277 :
1278 :
1279 : end if
1280 :
1281 :
1282 33684 : end subroutine brem_neu
1283 :
1284 :
1285 33684 : subroutine reco_neu(sreco,srecodt,srecodd,srecoda,srecodz, input)
1286 : real(dp), intent(out) :: sreco,srecodt,srecodd,srecoda,srecodz
1287 : type(inputs), intent(in) :: input
1288 :
1289 33684 : real(dp) :: a0,a1,a2,a3,c00,c01, &
1290 33684 : dd00,dd01, &
1291 33684 : b,c,d,f1,f2,f3,z, &
1292 33684 : dum,dumdt,dumdd,dumda,dumdz, &
1293 33684 : gum,gumdt,gumdd,gumda,gumdz
1294 :
1295 33684 : real(dp) :: zeta,zetadt,zetadd,zetada,zetadz
1296 33684 : real(dp) :: xnum,xnumdt,xnumdd,xnumda,xnumdz
1297 :
1298 33684 : real(dp) :: nu,nudt,nudd,nuda,nudz, &
1299 33684 : nu2,nu3,bigj,bigjdt,bigjdd,bigjda,bigjdz
1300 :
1301 33684 : sreco=0d0
1302 33684 : srecodt=0d0
1303 33684 : srecodd=0d0
1304 33684 : srecoda=0d0
1305 33684 : srecodz=0d0
1306 :
1307 : !..recombination neutrino section
1308 : !..for reactions like e- (continuum) => e- (bound) + nu_e + nubar_e
1309 :
1310 : !..equation 6.11 solved for nu
1311 33684 : xnum = 1.10520d8 * input% den * input% ye /(input% temp*sqrt(input% temp))
1312 33684 : xnumdt = -1.50d0*xnum*input% tempi
1313 33684 : xnumdd = xnum * input% deni
1314 33684 : xnumda = -xnum*input% abari
1315 33684 : xnumdz = xnum*input% zbari
1316 :
1317 : !..the chemical potential
1318 33684 : nu = ifermi12(xnum)
1319 :
1320 : !..a0 is d(nu)/d(xnum)
1321 33684 : a0 = 1.0d0/(0.5d0*zfermim12(nu))
1322 33684 : nudt = a0*xnumdt
1323 33684 : nudd = a0*xnumdd
1324 33684 : nuda = a0*xnumda
1325 33684 : nudz = a0*xnumdz
1326 :
1327 33684 : nu2 = nu * nu
1328 33684 : nu3 = nu2 * nu
1329 :
1330 : !..table 12
1331 33684 : if (nu >= -20.0d0 .and. nu < 0.0d0) then
1332 : a1 = 1.51d-2
1333 : a2 = 2.42d-1
1334 : a3 = 1.21d0
1335 : b = 3.71d-2
1336 : c = 9.06d-1
1337 : d = 9.28d-1
1338 : f1 = 0.0d0
1339 : f2 = 0.0d0
1340 : f3 = 0.0d0
1341 6446 : else if (nu >= 0.0d0 .and. nu <= 10.0d0) then
1342 33684 : a1 = 1.23d-2
1343 33684 : a2 = 2.66d-1
1344 33684 : a3 = 1.30d0
1345 33684 : b = 1.17d-1
1346 33684 : c = 8.97d-1
1347 33684 : d = 1.77d-1
1348 33684 : f1 = -1.20d-2
1349 33684 : f2 = 2.29d-2
1350 33684 : f3 = -1.04d-3
1351 : end if
1352 :
1353 :
1354 : !..equation 6.7, 6.13 and 6.14
1355 33684 : if (nu >= -20.0d0 .and. nu <= 10.0d0) then
1356 :
1357 28471 : zeta = 1.579d5*input% zbar*input% zbar*input% tempi
1358 28471 : zetadt = -zeta*input% tempi
1359 28471 : zetadd = 0.0d0
1360 28471 : zetada = 0.0d0
1361 28471 : zetadz = 2.0d0*zeta*input% zbari
1362 :
1363 28471 : c00 = 1.0d0/(1.0d0 + f1*nu + f2*nu2 + f3*nu3)
1364 28471 : c01 = f1 + f2*2.0d0*nu + f3*3.0d0*nu2
1365 28471 : dum = zeta*c00
1366 28471 : dumdt = zetadt*c00 -1d0 * c00 * c00 * zeta*c01*nudt
1367 28471 : dumdd = -1d0 * c00 * c00 * zeta*c01*nudd
1368 28471 : dumda = -1d0 * c00 * c00 * zeta*c01*nuda
1369 28471 : dumdz = zetadz*c00 -1d0 * c00 *c00 * zeta*c01*nudz
1370 :
1371 28471 : z = 1.0d0/dum
1372 28471 : dd00 = pow(dum,-2.25d0)
1373 28471 : dd01 = pow(dum,-4.55d0)
1374 28471 : c00 = a1*z + a2*dd00 + a3*dd01
1375 28471 : c01 = -(a1*z + 2.25d0*a2*dd00 + 4.55d0*a3*dd01)*z
1376 :
1377 28471 : z = exp(c*nu)
1378 28471 : dd00 = b*z*(1.0d0 + d*dum)
1379 28471 : gum = 1.0d0 + dd00
1380 28471 : gumdt = dd00*c*nudt + b*z*d*dumdt
1381 28471 : gumdd = dd00*c*nudd + b*z*d*dumdd
1382 28471 : gumda = dd00*c*nuda + b*z*d*dumda
1383 28471 : gumdz = dd00*c*nudz + b*z*d*dumdz
1384 :
1385 28471 : z = exp(nu)
1386 28471 : a1 = 1.0d0/gum
1387 :
1388 28471 : bigj = c00 * z * a1
1389 28471 : bigjdt = c01*dumdt*z*a1 + c00*z*nudt*a1 - c00*z*a1*a1 * gumdt
1390 28471 : bigjdd = c01*dumdd*z*a1 + c00*z*nudd*a1 - c00*z*a1*a1 * gumdd
1391 28471 : bigjda = c01*dumda*z*a1 + c00*z*nuda*a1 - c00*z*a1*a1 * gumda
1392 28471 : bigjdz = c01*dumdz*z*a1 + c00*z*nudz*a1 - c00*z*a1*a1 * gumdz
1393 :
1394 : !..equation 6.5
1395 28471 : z = exp(zeta + nu)
1396 28471 : dum = 1.0d0 + z
1397 28471 : a1 = 1.0d0/dum
1398 28471 : a2 = 1.0d0/bigj
1399 :
1400 28471 : sreco = tfac6 * 2.649d-18 * input% ye * pow(input% zbar,13) * input% den * bigj*a1
1401 28471 : srecodt = sreco*(bigjdt*a2 - z*(zetadt + nudt)*a1)
1402 28471 : srecodd = sreco*(1.0d0 * input% deni + bigjdd*a2 - z*(zetadd + nudd)*a1)
1403 28471 : srecoda = sreco*(-1.0d0*input% abari + bigjda*a2 - z*(zetada+nuda)*a1)
1404 28471 : srecodz = sreco*(14.0d0*input% zbari + bigjdz*a2 - z*(zetadz+nudz)*a1)
1405 :
1406 : end if
1407 :
1408 :
1409 33684 : end subroutine reco_neu
1410 :
1411 33684 : subroutine plas_neu(splas,splasdt,splasdd,splasda,splasdz, input)
1412 : real(dp), intent(out) :: splas,splasdt,splasdd,splasda,splasdz
1413 : type(inputs), intent(in) :: input
1414 :
1415 33684 : real(dp) :: a1,a2,a3,b1,b2,c00,c01,c02,c03,c04,xlnt,cc, &
1416 33684 : c,d,f1,&
1417 33684 : dumdt,dumdd,dumda,dumdz, gum
1418 :
1419 33684 : real(dp) :: xnum,xnumdt,xnumdd,xnumda,xnumdz, &
1420 33684 : xden,xdendt,xdendd,xdenda,xdendz
1421 :
1422 33684 : real(dp) :: gl,gl2,gl2dt,gl2dd,gl2da,gl2dz,gl12,gl32,gl72,gl6, &
1423 33684 : ft,ftdt,ftdd,ftda,ftdz,fl,fldt,fldd,flda,fldz, &
1424 33684 : fxy,fxydt,fxydd,fxyda,fxydz
1425 :
1426 33684 : splas=0d0
1427 33684 : splasdt=0d0
1428 33684 : splasdd=0d0
1429 33684 : splasda=0d0
1430 33684 : splasdz=0d0
1431 :
1432 : !..plasma neutrino section
1433 : !..for collective reactions like gamma_plasmon => nu_e + nubar_e
1434 : !..equation 4.6
1435 :
1436 33684 : a1 = 1.019d-6*input% rm
1437 33684 : a2 = pow(a1,two_thirds)
1438 33684 : a3 = two_thirds*a2/a1
1439 :
1440 33684 : b1 = sqrt(1.0d0 + a2)
1441 33684 : b2 = 1.0d0/b1
1442 :
1443 33684 : c00 = 1.0d0/(input% temp*input% temp*b1)
1444 :
1445 33684 : gl2 = 1.1095d11 * input% rm * c00
1446 :
1447 33684 : gl2dt = -2.0d0*gl2*input% tempi
1448 33684 : d = input% rm*c00*b2*0.5d0*b2*a3*1.019d-6
1449 33684 : gl2dd = 1.1095d11 * (input% rmdd*c00 - d*input% rmdd)
1450 33684 : gl2da = 1.1095d11 * (input% rmda*c00 - d*input% rmda)
1451 33684 : gl2dz = 1.1095d11 * (input% rmdz*c00 - d*input% rmdz)
1452 :
1453 :
1454 33684 : gl = sqrt(gl2)
1455 33684 : gl12 = sqrt(gl)
1456 33684 : gl32 = gl * gl12
1457 33684 : gl72 = gl2 * gl32
1458 33684 : gl6 = gl2 * gl2 * gl2
1459 :
1460 :
1461 : !..equation 4.7
1462 33684 : ft = 2.4d0 + 0.6d0*gl12 + 0.51d0*gl + 1.25d0*gl32
1463 33684 : gum = 1.0d0/gl2
1464 33684 : a1 =(0.25d0*0.6d0*gl12 +0.5d0*0.51d0*gl +0.75d0*1.25d0*gl32)*gum
1465 33684 : ftdt = a1*gl2dt
1466 33684 : ftdd = a1*gl2dd
1467 33684 : ftda = a1*gl2da
1468 33684 : ftdz = a1*gl2dz
1469 :
1470 :
1471 : !..equation 4.8
1472 33684 : a1 = 8.6d0*gl2 + 1.35d0*gl72
1473 33684 : a2 = 8.6d0 + 1.75d0*1.35d0*gl72*gum
1474 :
1475 33684 : b1 = 225.0d0 - 17.0d0*gl + gl2
1476 33684 : b2 = -0.5d0*17.0d0*gl*gum + 1.0d0
1477 :
1478 33684 : c = 1.0d0/b1
1479 33684 : fl = a1*c
1480 :
1481 33684 : d = (a2 - fl*b2)*c
1482 33684 : fldt = d*gl2dt
1483 33684 : fldd = d*gl2dd
1484 33684 : flda = d*gl2da
1485 33684 : fldz = d*gl2dz
1486 :
1487 :
1488 : !..equation 4.9 and 4.10
1489 33684 : cc = log10(2.0d0*input% rm)
1490 33684 : xlnt = input% logtemp
1491 :
1492 33684 : xnum = one_sixth * (17.5d0 + cc - 3.0d0*xlnt)
1493 33684 : xnumdt = -iln10*0.5d0*input% tempi
1494 33684 : a2 = iln10*one_sixth*input% rmi
1495 33684 : xnumdd = a2*input% rmdd
1496 33684 : xnumda = a2*input% rmda
1497 33684 : xnumdz = a2*input% rmdz
1498 :
1499 33684 : xden = one_sixth * (-24.5d0 + cc + 3.0d0*xlnt)
1500 33684 : xdendt = iln10*0.5d0*input% tempi
1501 33684 : xdendd = a2*input% rmdd
1502 33684 : xdenda = a2*input% rmda
1503 33684 : xdendz = a2*input% rmdz
1504 :
1505 :
1506 : !..equation 4.11
1507 33684 : if (abs(xnum) > 0.7d0 .or. xden < 0.0d0) then
1508 : fxy = 1.0d0
1509 : fxydt = 0.0d0
1510 : fxydd = 0.0d0
1511 : fxydz = 0.0d0
1512 : fxyda = 0.0d0
1513 :
1514 : else
1515 :
1516 5642 : a1 = 0.39d0 - 1.25d0*xnum - 0.35d0*sin(4.5d0*xnum)
1517 5642 : a2 = -1.25d0 - 4.5d0*0.35d0*cos(4.5d0*xnum)
1518 :
1519 5642 : b1 = 0.3d0 * exp(-1.0d0*pow((4.5d0*xnum + 0.9d0),2))
1520 5642 : b2 = -b1*2.0d0*(4.5d0*xnum + 0.9d0)*4.5d0
1521 :
1522 5642 : c = min(0.0d0, xden - 1.6d0 + 1.25d0*xnum)
1523 5642 : if (c == 0.0d0) then
1524 : dumdt = 0.0d0
1525 : dumdd = 0.0d0
1526 : dumda = 0.0d0
1527 : dumdz = 0.0d0
1528 : else
1529 3045 : dumdt = xdendt + 1.25d0*xnumdt
1530 3045 : dumdd = xdendd + 1.25d0*xnumdd
1531 3045 : dumda = xdenda + 1.25d0*xnumda
1532 3045 : dumdz = xdendz + 1.25d0*xnumdz
1533 : end if
1534 :
1535 5642 : d = 0.57d0 - 0.25d0*xnum
1536 5642 : a3 = c/d
1537 5642 : c00 = exp(-1.0d0*a3*a3)
1538 :
1539 5642 : f1 = -c00*2.0d0*a3/d
1540 5642 : c01 = f1*(dumdt + a3*0.25d0*xnumdt)
1541 5642 : c02 = f1*(dumdd + a3*0.25d0*xnumdd)
1542 5642 : c03 = f1*(dumda + a3*0.25d0*xnumda)
1543 5642 : c04 = f1*(dumdz + a3*0.25d0*xnumdz)
1544 :
1545 5642 : fxy = 1.05d0 + (a1 - b1)*c00
1546 5642 : fxydt = (a2*xnumdt - b2*xnumdt)*c00 + (a1-b1)*c01
1547 5642 : fxydd = (a2*xnumdd - b2*xnumdd)*c00 + (a1-b1)*c02
1548 5642 : fxyda = (a2*xnumda - b2*xnumda)*c00 + (a1-b1)*c03
1549 5642 : fxydz = (a2*xnumdz - b2*xnumdz)*c00 + (a1-b1)*c04
1550 :
1551 : end if
1552 :
1553 :
1554 : !..equation 4.1 and 4.5
1555 33684 : splas = (ft + fl) * fxy
1556 33684 : splasdt = (ftdt + fldt)*fxy + (ft+fl)*fxydt
1557 33684 : splasdd = (ftdd + fldd)*fxy + (ft+fl)*fxydd
1558 33684 : splasda = (ftda + flda)*fxy + (ft+fl)*fxyda
1559 33684 : splasdz = (ftdz + fldz)*fxy + (ft+fl)*fxydz
1560 :
1561 33684 : a2 = exp(-gl)
1562 33684 : a3 = -0.5d0*a2*gl*gum
1563 :
1564 : a1 = splas
1565 33684 : splas = a2*a1
1566 33684 : splasdt = a2*splasdt + a3*gl2dt*a1
1567 33684 : splasdd = a2*splasdd + a3*gl2dd*a1
1568 33684 : splasda = a2*splasda + a3*gl2da*a1
1569 33684 : splasdz = a2*splasdz + a3*gl2dz*a1
1570 :
1571 33684 : a2 = gl6
1572 33684 : a3 = 3.0d0*gl6*gum
1573 :
1574 : a1 = splas
1575 33684 : splas = a2*a1
1576 33684 : splasdt = a2*splasdt + a3*gl2dt*a1
1577 33684 : splasdd = a2*splasdd + a3*gl2dd*a1
1578 33684 : splasda = a2*splasda + a3*gl2da*a1
1579 33684 : splasdz = a2*splasdz + a3*gl2dz*a1
1580 :
1581 :
1582 33684 : a2 = 0.93153d0 * 3.0d21 * input% xl9
1583 33684 : a3 = 0.93153d0 * 3.0d21 * 9.0d0*input% xl8*input% xldt
1584 :
1585 : a1 = splas
1586 33684 : splas = a2*a1
1587 33684 : splasdt = a2*splasdt + a3*a1
1588 33684 : splasdd = a2*splasdd
1589 33684 : splasda = a2*splasda
1590 33684 : splasdz = a2*splasdz
1591 :
1592 :
1593 33684 : end subroutine plas_neu
1594 :
1595 :
1596 33684 : subroutine pair_neu(spair,spairdt,spairdd,spairda,spairdz, input)
1597 : real(dp), intent(out) :: spair,spairdt,spairdd,spairda,spairdz
1598 : type(inputs), intent(in) :: input
1599 :
1600 33684 : real(dp) :: a1,a2,a3,b1,b2,c,d, gl,gldt
1601 :
1602 33684 : real(dp) :: xnum,xnumdt,xnumdd,xnumda,xnumdz, &
1603 33684 : xden,xdendt,xdendd,xdenda,xdendz
1604 :
1605 33684 : real(dp) :: fpair,fpairdt,fpairdd,fpairda,fpairdz, &
1606 33684 : qpair,qpairdt,qpairdd,qpairda,qpairdz
1607 :
1608 : !..pair neutrino section
1609 : !..for reactions like e+ + e- => nu_e + nubar_e
1610 :
1611 : !..equation 2.8
1612 33684 : gl = 1.0d0 - 13.04d0*input% xl2 +133.5d0*input% xl4 +1534.0d0*input% xl6 +918.6d0*input% xl8
1613 33684 : gldt = input% xldt*(-26.08d0*input% xl +534.0d0*input% xl3 +9204.0d0*input% xl5 +7348.8d0*input% xl7)
1614 :
1615 : !..equation 2.7
1616 :
1617 33684 : a1 = 6.002d19 + 2.084d20*input% zeta + 1.872d21*input% zeta2
1618 33684 : a2 = 2.084d20 + 2.0d0*1.872d21*input% zeta
1619 :
1620 33684 : if (input% t9 < 10.0d0) then
1621 32684 : b1 = exp(-5.5924d0*input% zeta)
1622 32684 : b2 = -b1*5.5924d0
1623 : else
1624 1000 : b1 = exp(-4.9924d0*input% zeta)
1625 1000 : b2 = -b1*4.9924d0
1626 : end if
1627 :
1628 33684 : xnum = a1 * b1
1629 33684 : c = a2*b1 + a1*b2
1630 33684 : xnumdt = c*input% zetadt
1631 33684 : xnumdd = c*input% zetadd
1632 33684 : xnumda = c*input% zetada
1633 33684 : xnumdz = c*input% zetadz
1634 :
1635 33684 : if (input% t9 < 10.0d0) then
1636 32684 : a1 = 9.383d-1*input% xlm1 - 4.141d-1*input% xlm2 + 5.829d-2*input% xlm3
1637 32684 : a2 = -9.383d-1*input% xlm2 + 2.0d0*4.141d-1*input% xlm3 - 3.0d0*5.829d-2*input% xlm4
1638 : else
1639 1000 : a1 = 1.2383d0*input% xlm1 - 8.141d-1*input% xlm2
1640 1000 : a2 = -1.2383d0*input% xlm2 + 2.0d0*8.141d-1*input% xlm3
1641 : end if
1642 :
1643 33684 : b1 = 3.0d0*input% zeta2
1644 :
1645 33684 : xden = input% zeta3 + a1
1646 33684 : xdendt = b1*input% zetadt + a2*input% xldt
1647 33684 : xdendd = b1*input% zetadd
1648 33684 : xdenda = b1*input% zetada
1649 33684 : xdendz = b1*input% zetadz
1650 :
1651 33684 : a1 = 1.0d0/xden
1652 33684 : fpair = xnum*a1
1653 33684 : fpairdt = (xnumdt - fpair*xdendt)*a1
1654 33684 : fpairdd = (xnumdd - fpair*xdendd)*a1
1655 33684 : fpairda = (xnumda - fpair*xdenda)*a1
1656 33684 : fpairdz = (xnumdz - fpair*xdendz)*a1
1657 :
1658 :
1659 : !..equation 2.6
1660 33684 : a1 = 10.7480d0*input% xl2 + 0.3967d0*input% xlp5 + 1.005d0
1661 33684 : a2 = input% xldt*(2.0d0*10.7480d0*input% xl + 0.5d0*0.3967d0*input% xlmp5)
1662 33684 : xnum = 1.0d0/a1
1663 33684 : xnumdt = -xnum*xnum*a2
1664 :
1665 33684 : a1 = 7.692d7*input% xl3 + 9.715d6*input% xlp5
1666 33684 : a2 = input% xldt*(3.0d0*7.692d7*input% xl2 + 0.5d0*9.715d6*input% xlmp5)
1667 :
1668 33684 : c = 1.0d0/a1
1669 33684 : b1 = 1.0d0 + input% rm*c
1670 :
1671 33684 : xden = pow(b1,-0.3d0)
1672 :
1673 33684 : d = -0.3d0*xden/b1
1674 33684 : xdendt = -d*input% rm*c*c*a2
1675 33684 : xdendd = d*input% rmdd*c
1676 33684 : xdenda = d*input% rmda*c
1677 33684 : xdendz = d*input% rmdz*c
1678 :
1679 33684 : qpair = xnum*xden
1680 33684 : qpairdt = xnumdt*xden + xnum*xdendt
1681 33684 : qpairdd = xnum*xdendd
1682 33684 : qpairda = xnum*xdenda
1683 33684 : qpairdz = xnum*xdendz
1684 :
1685 : !..equation 2.5
1686 33684 : a1 = exp(-2.0d0*input% xlm1)
1687 33684 : a2 = a1*2.0d0*input% xlm2*input% xldt
1688 :
1689 33684 : spair = a1*fpair
1690 33684 : spairdt = a2*fpair + a1*fpairdt
1691 33684 : spairdd = a1*fpairdd
1692 33684 : spairda = a1*fpairda
1693 33684 : spairdz = a1*fpairdz
1694 :
1695 33684 : a1 = spair
1696 33684 : spair = gl*a1
1697 33684 : spairdt = gl*spairdt + gldt*a1
1698 33684 : spairdd = gl*spairdd
1699 33684 : spairda = gl*spairda
1700 33684 : spairdz = gl*spairdz
1701 :
1702 33684 : a1 = tfac4*(1.0d0 + tfac3 * qpair)
1703 33684 : a2 = tfac4*tfac3
1704 :
1705 33684 : a3 = spair
1706 33684 : spair = a1*a3
1707 33684 : spairdt = a1*spairdt + a2*qpairdt*a3
1708 33684 : spairdd = a1*spairdd + a2*qpairdd*a3
1709 33684 : spairda = a1*spairda + a2*qpairda*a3
1710 33684 : spairdz = a1*spairdz + a2*qpairdz*a3
1711 :
1712 33684 : end subroutine pair_neu
1713 :
1714 :
1715 0 : end module mod_neu
|
