sneut5.f90
Go to the documentation of this file.
1  subroutine sneut5_aa(temp,den,abar,zbar, &
2  snu,dsnudt,dsnudd,dsnuda,dsnudz, &
3  spair,splas,sphot,sbrem,sreco)
4  include 'implno.dek'
5  include 'const.dek'
6 
7 ! this routine computes thermal neutrino losses from the analytic fits of
8 ! itoh et al. apjs 102, 411, 1996 and also returns their derivatives.
9 
10 ! input:
11 ! temp = temperature
12 ! den = density
13 ! abar = mean atomic weight
14 ! zbar = mean charge
15 
16 ! output:
17 ! snu = total neutrino loss rate in erg/g/sec
18 ! dsnudt = derivative of snu with temperature
19 ! dsnudd = derivative of snu with density
20 ! dsnuda = derivative of snu with abar
21 ! dsnudz = derivative of snu with zbar
22 
23 
24 ! declare the pass
25  double precision temp,den,abar,zbar, &
26  snu,dsnudt,dsnudd,dsnuda,dsnudz
27 
28 ! local variables
29  double precision spair,spairdt,spairdd,spairda,spairdz, &
30  splas,splasdt,splasdd,splasda,splasdz, &
31  sphot,sphotdt,sphotdd,sphotda,sphotdz, &
32  sbrem,sbremdt,sbremdd,sbremda,sbremdz, &
33  sreco,srecodt,srecodd,srecoda,srecodz
34 
35  double precision t9,xl,xldt,xlp5,xl2,xl3,xl4,xl5,xl6,xl7,xl8,xl9, &
36  xlmp5,xlm1,xlm2,xlm3,xlm4,xlnt,cc,den6,tfermi, &
37  a0,a1,a2,a3,b1,b2,c00,c01,c02,c03,c04,c05,c06, &
38  c10,c11,c12,c13,c14,c15,c16,c20,c21,c22,c23,c24, &
39  c25,c26,dd00,dd01,dd02,dd03,dd04,dd05,dd11,dd12, &
40  dd13,dd14,dd15,dd21,dd22,dd23,dd24,dd25,b,c,d,f0, &
41  f1,deni,tempi,abari,zbari,f2,f3,z,xmue,ye, &
42  dum,dumdt,dumdd,dumda,dumdz, &
43  gum,gumdt,gumdd,gumda,gumdz
44 
45 
46 ! pair production
47  double precision rm,rmdd,rmda,rmdz,rmi,gl,gldt, &
48  zeta,zetadt,zetadd,zetada,zetadz,zeta2,zeta3, &
49  xnum,xnumdt,xnumdd,xnumda,xnumdz, &
50  xden,xdendt,xdendd,xdenda,xdendz, &
51  fpair,fpairdt,fpairdd,fpairda,fpairdz, &
52  qpair,qpairdt,qpairdd,qpairda,qpairdz
53 
54 ! plasma
55  double precision gl2,gl2dt,gl2dd,gl2da,gl2dz,gl12,gl32,gl72,gl6, &
56  ft,ftdt,ftdd,ftda,ftdz,fl,fldt,fldd,flda,fldz, &
57  fxy,fxydt,fxydd,fxyda,fxydz
58 
59 ! photo
60  double precision tau,taudt,cos1,cos2,cos3,cos4,cos5,sin1,sin2, &
61  sin3,sin4,sin5,last,xast, &
62  fphot,fphotdt,fphotdd,fphotda,fphotdz, &
63  qphot,qphotdt,qphotdd,qphotda,qphotdz
64 
65 ! brem
66  double precision t8,t812,t832,t82,t83,t85,t86,t8m1,t8m2,t8m3,t8m5, &
67  t8m6, &
68  eta,etadt,etadd,etada,etadz,etam1,etam2,etam3, &
69  fbrem,fbremdt,fbremdd,fbremda,fbremdz, &
70  gbrem,gbremdt,gbremdd,gbremda,gbremdz, &
71  u,gm1,gm2,gm13,gm23,gm43,gm53,v,w,fb,gt,gb, &
72  fliq,fliqdt,fliqdd,fliqda,fliqdz, &
73  gliq,gliqdt,gliqdd,gliqda,gliqdz
74 
75 ! recomb
76  double precision ifermi12,zfermim12,nu,nudt,nudd,nuda,nudz, &
77  nu2,nu3,bigj,bigjdt,bigjdd,bigjda,bigjdz
78 
79 
80 
81 ! numerical constants
82  double precision fac1,fac2,fac3,oneth,twoth,con1,sixth,iln10
83  parameter (fac1 = 5.0d0 * pi / 3.0d0, &
84  fac2 = 10.0d0 * pi, &
85  fac3 = pi / 5.0d0, &
86  oneth = 1.0d0/3.0d0, &
87  twoth = 2.0d0/3.0d0, &
88  con1 = 1.0d0/5.9302d0, &
89  sixth = 1.0d0/6.0d0, &
90  iln10 = 4.342944819032518d-1)
91 
92 
93 ! theta is sin**2(theta_weinberg) = 0.2223 (2019)
94 ! xnufam is the number of neutrino flavors = 3.02 plus/minus 0.005 (1998)
95 ! change theta and xnufam if need be, and the changes will automatically
96 ! propagate through the routine. cv and ca are the vector and axial currents.
97 
98  double precision theta,xnufam,cv,ca,cvp,cap,tfac1,tfac2,tfac3, &
99  tfac4,tfac5,tfac6
100  parameter(theta = 0.2223d0, &
101  xnufam = 3.0d0, &
102  cv = 0.5d0 + 2.0d0 * theta, &
103  cvp = 1.0d0 - cv, &
104  ca = 0.5d0, &
105  cap = 1.0d0 - ca, &
106  tfac1 = cv*cv + ca*ca + &
107  (xnufam-1.0d0) * (cvp*cvp+cap*cap), &
108  tfac2 = cv*cv - ca*ca + &
109  (xnufam-1.0d0) * (cvp*cvp - cap*cap), &
110  tfac3 = tfac2/tfac1, &
111  tfac4 = 0.5d0 * tfac1, &
112  tfac5 = 0.5d0 * tfac2, &
113  tfac6 = cv*cv + 1.5d0*ca*ca + (xnufam - 1.0d0)* &
114  (cvp*cvp + 1.5d0*cap*cap))
115 
116 
117 
118 ! initialize
119  spair = 0.0d0
120  spairdt = 0.0d0
121  spairdd = 0.0d0
122  spairda = 0.0d0
123  spairdz = 0.0d0
124 
125  splas = 0.0d0
126  splasdt = 0.0d0
127  splasdd = 0.0d0
128  splasda = 0.0d0
129  splasdz = 0.0d0
130 
131  sphot = 0.0d0
132  sphotdt = 0.0d0
133  sphotdd = 0.0d0
134  sphotda = 0.0d0
135  sphotdz = 0.0d0
136 
137  sbrem = 0.0d0
138  sbremdt = 0.0d0
139  sbremdd = 0.0d0
140  sbremda = 0.0d0
141  sbremdz = 0.0d0
142 
143  sreco = 0.0d0
144  srecodt = 0.0d0
145  srecodd = 0.0d0
146  srecoda = 0.0d0
147  srecodz = 0.0d0
148 
149  snu = 0.0d0
150  dsnudt = 0.0d0
151  dsnudd = 0.0d0
152  dsnuda = 0.0d0
153  dsnudz = 0.0d0
154 
155  if (temp .lt. 1.0e7) return
156 
157 
158 ! to avoid lots of divisions
159  deni = 1.0d0/den
160  tempi = 1.0d0/temp
161  abari = 1.0d0/abar
162  zbari = 1.0d0/zbar
163 
164 
165 ! some composition variables
166  ye = zbar*abari
167  xmue = abar*zbari
168 
169 
170 
171 
172 ! some frequent factors
173  t9 = temp * 1.0d-9
174  xl = t9 * con1
175  xldt = 1.0d-9 * con1
176  xlp5 = sqrt(xl)
177  xl2 = xl*xl
178  xl3 = xl2*xl
179  xl4 = xl3*xl
180  xl5 = xl4*xl
181  xl6 = xl5*xl
182  xl7 = xl6*xl
183  xl8 = xl7*xl
184  xl9 = xl8*xl
185  xlmp5 = 1.0d0/xlp5
186  xlm1 = 1.0d0/xl
187  xlm2 = xlm1*xlm1
188  xlm3 = xlm1*xlm2
189  xlm4 = xlm1*xlm3
190 
191  rm = den*ye
192  rmdd = ye
193  rmda = -rm*abari
194  rmdz = den*abari
195  rmi = 1.0d0/rm
196 
197  a0 = rm * 1.0d-9
198  a1 = a0**oneth
199  zeta = a1 * xlm1
200  zetadt = -a1 * xlm2 * xldt
201  a2 = oneth * a1*rmi * xlm1
202  zetadd = a2 * rmdd
203  zetada = a2 * rmda
204  zetadz = a2 * rmdz
205 
206  zeta2 = zeta * zeta
207  zeta3 = zeta2 * zeta
208 
209 
210 
211 
212 ! pair neutrino section
213 ! for reactions like e+ + e- => nu_e + nubar_e
214 
215 ! equation 2.8
216  gl = 1.0d0 - 13.04d0*xl2 +133.5d0*xl4 +1534.0d0*xl6 +918.6d0*xl8
217  gldt = xldt*(-26.08d0*xl +534.0d0*xl3 +9204.0d0*xl5 +7348.8d0*xl7)
218 
219 ! equation 2.7
220 
221  a1 = 6.002d19 + 2.084d20*zeta + 1.872d21*zeta2
222  a2 = 2.084d20 + 2.0d0*1.872d21*zeta
223 
224  if (t9 .lt. 10.0) then
225  b1 = exp(-5.5924d0*zeta)
226  b2 = -b1*5.5924d0
227  else
228  b1 = exp(-4.9924d0*zeta)
229  b2 = -b1*4.9924d0
230  end if
231 
232  xnum = a1 * b1
233  c = a2*b1 + a1*b2
234  xnumdt = c*zetadt
235  xnumdd = c*zetadd
236  xnumda = c*zetada
237  xnumdz = c*zetadz
238 
239  if (t9 .lt. 10.0) then
240  a1 = 9.383d-1*xlm1 - 4.141d-1*xlm2 + 5.829d-2*xlm3
241  a2 = -9.383d-1*xlm2 + 2.0d0*4.141d-1*xlm3 - 3.0d0*5.829d-2*xlm4
242  else
243  a1 = 1.2383d0*xlm1 - 8.141d-1*xlm2
244  a2 = -1.2383d0*xlm2 + 2.0d0*8.141d-1*xlm3
245  end if
246 
247  b1 = 3.0d0*zeta2
248 
249  xden = zeta3 + a1
250  xdendt = b1*zetadt + a2*xldt
251  xdendd = b1*zetadd
252  xdenda = b1*zetada
253  xdendz = b1*zetadz
254 
255  a1 = 1.0d0/xden
256  fpair = xnum*a1
257  fpairdt = (xnumdt - fpair*xdendt)*a1
258  fpairdd = (xnumdd - fpair*xdendd)*a1
259  fpairda = (xnumda - fpair*xdenda)*a1
260  fpairdz = (xnumdz - fpair*xdendz)*a1
261 
262 
263 ! equation 2.6
264  a1 = 10.7480d0*xl2 + 0.3967d0*xlp5 + 1.005d0
265  a2 = xldt*(2.0d0*10.7480d0*xl + 0.5d0*0.3967d0*xlmp5)
266  xnum = 1.0d0/a1
267  xnumdt = -xnum*xnum*a2
268 
269  a1 = 7.692d7*xl3 + 9.715d6*xlp5
270  a2 = xldt*(3.0d0*7.692d7*xl2 + 0.5d0*9.715d6*xlmp5)
271 
272  c = 1.0d0/a1
273  b1 = 1.0d0 + rm*c
274 
275  xden = b1**(-0.3d0)
276 
277  d = -0.3d0*xden/b1
278  xdendt = -d*rm*c*c*a2
279  xdendd = d*rmdd*c
280  xdenda = d*rmda*c
281  xdendz = d*rmdz*c
282 
283  qpair = xnum*xden
284  qpairdt = xnumdt*xden + xnum*xdendt
285  qpairdd = xnum*xdendd
286  qpairda = xnum*xdenda
287  qpairdz = xnum*xdendz
288 
289 
290 
291 ! equation 2.5
292  a1 = exp(-2.0d0*xlm1)
293  a2 = a1*2.0d0*xlm2*xldt
294 
295  spair = a1*fpair
296  spairdt = a2*fpair + a1*fpairdt
297  spairdd = a1*fpairdd
298  spairda = a1*fpairda
299  spairdz = a1*fpairdz
300 
301  a1 = spair
302  spair = gl*a1
303  spairdt = gl*spairdt + gldt*a1
304  spairdd = gl*spairdd
305  spairda = gl*spairda
306  spairdz = gl*spairdz
307 
308  a1 = tfac4*(1.0d0 + tfac3 * qpair)
309  a2 = tfac4*tfac3
310 
311  a3 = spair
312  spair = a1*a3
313  spairdt = a1*spairdt + a2*qpairdt*a3
314  spairdd = a1*spairdd + a2*qpairdd*a3
315  spairda = a1*spairda + a2*qpairda*a3
316  spairdz = a1*spairdz + a2*qpairdz*a3
317 
318 
319 
320 
321 ! plasma neutrino section
322 ! for collective reactions like gamma_plasmon => nu_e + nubar_e
323 ! equation 4.6
324 
325  a1 = 1.019d-6*rm
326  a2 = a1**twoth
327  a3 = twoth*a2/a1
328 
329  b1 = sqrt(1.0d0 + a2)
330  b2 = 1.0d0/b1
331 
332  c00 = 1.0d0/(temp*temp*b1)
333 
334  gl2 = 1.1095d11 * rm * c00
335 
336  gl2dt = -2.0d0*gl2*tempi
337  d = rm*c00*b2*0.5d0*b2*a3*1.019d-6
338  gl2dd = 1.1095d11 * (rmdd*c00 - d*rmdd)
339  gl2da = 1.1095d11 * (rmda*c00 - d*rmda)
340  gl2dz = 1.1095d11 * (rmdz*c00 - d*rmdz)
341 
342 
343  gl = sqrt(gl2)
344  gl12 = sqrt(gl)
345  gl32 = gl * gl12
346  gl72 = gl2 * gl32
347  gl6 = gl2 * gl2 * gl2
348 
349 
350 ! equation 4.7
351  ft = 2.4d0 + 0.6d0*gl12 + 0.51d0*gl + 1.25d0*gl32
352  gum = 1.0d0/gl2
353  a1 =(0.25d0*0.6d0*gl12 +0.5d0*0.51d0*gl +0.75d0*1.25d0*gl32)*gum
354  ftdt = a1*gl2dt
355  ftdd = a1*gl2dd
356  ftda = a1*gl2da
357  ftdz = a1*gl2dz
358 
359 
360 ! equation 4.8
361  a1 = 8.6d0*gl2 + 1.35d0*gl72
362  a2 = 8.6d0 + 1.75d0*1.35d0*gl72*gum
363 
364  b1 = 225.0d0 - 17.0d0*gl + gl2
365  b2 = -0.5d0*17.0d0*gl*gum + 1.0d0
366 
367  c = 1.0d0/b1
368  fl = a1*c
369 
370  d = (a2 - fl*b2)*c
371  fldt = d*gl2dt
372  fldd = d*gl2dd
373  flda = d*gl2da
374  fldz = d*gl2dz
375 
376 
377 ! equation 4.9 and 4.10
378  cc = log10(2.0d0*rm)
379  xlnt = log10(temp)
380 
381  xnum = sixth * (17.5d0 + cc - 3.0d0*xlnt)
382  xnumdt = -iln10*0.5d0*tempi
383  a2 = iln10*sixth*rmi
384  xnumdd = a2*rmdd
385  xnumda = a2*rmda
386  xnumdz = a2*rmdz
387 
388  xden = sixth * (-24.5d0 + cc + 3.0d0*xlnt)
389  xdendt = iln10*0.5d0*tempi
390  xdendd = a2*rmdd
391  xdenda = a2*rmda
392  xdendz = a2*rmdz
393 
394 
395 ! equation 4.11
396  if (abs(xnum) .gt. 0.7d0 .or. xden .lt. 0.0d0) then
397  fxy = 1.0d0
398  fxydt = 0.0d0
399  fxydd = 0.0d0
400  fxydz = 0.0d0
401  fxyda = 0.0d0
402 
403  else
404 
405  a1 = 0.39d0 - 1.25d0*xnum - 0.35d0*sin(4.5d0*xnum)
406  a2 = -1.25d0 - 4.5d0*0.35d0*cos(4.5d0*xnum)
407 
408  b1 = 0.3d0 * exp(-1.0d0*(4.5d0*xnum + 0.9d0)**2)
409  b2 = -b1*2.0d0*(4.5d0*xnum + 0.9d0)*4.5d0
410 
411  c = min(0.0d0, xden - 1.6d0 + 1.25d0*xnum)
412  if (c .eq. 0.0) then
413  dumdt = 0.0d0
414  dumdd = 0.0d0
415  dumda = 0.0d0
416  dumdz = 0.0d0
417  else
418  dumdt = xdendt + 1.25d0*xnumdt
419  dumdd = xdendd + 1.25d0*xnumdd
420  dumda = xdenda + 1.25d0*xnumda
421  dumdz = xdendz + 1.25d0*xnumdz
422  end if
423 
424  d = 0.57d0 - 0.25d0*xnum
425  a3 = c/d
426  c00 = exp(-1.0d0*a3**2)
427 
428  f1 = -c00*2.0d0*a3/d
429  c01 = f1*(dumdt + a3*0.25d0*xnumdt)
430  c02 = f1*(dumdd + a3*0.25d0*xnumdd)
431  c03 = f1*(dumda + a3*0.25d0*xnumda)
432  c04 = f1*(dumdz + a3*0.25d0*xnumdz)
433 
434  fxy = 1.05d0 + (a1 - b1)*c00
435  fxydt = (a2*xnumdt - b2*xnumdt)*c00 + (a1-b1)*c01
436  fxydd = (a2*xnumdd - b2*xnumdd)*c00 + (a1-b1)*c02
437  fxyda = (a2*xnumda - b2*xnumda)*c00 + (a1-b1)*c03
438  fxydz = (a2*xnumdz - b2*xnumdz)*c00 + (a1-b1)*c04
439 
440  end if
441 
442 
443 
444 ! equation 4.1 and 4.5
445  splas = (ft + fl) * fxy
446  splasdt = (ftdt + fldt)*fxy + (ft+fl)*fxydt
447  splasdd = (ftdd + fldd)*fxy + (ft+fl)*fxydd
448  splasda = (ftda + flda)*fxy + (ft+fl)*fxyda
449  splasdz = (ftdz + fldz)*fxy + (ft+fl)*fxydz
450 
451  a2 = exp(-gl)
452  a3 = -0.5d0*a2*gl*gum
453 
454  a1 = splas
455  splas = a2*a1
456  splasdt = a2*splasdt + a3*gl2dt*a1
457  splasdd = a2*splasdd + a3*gl2dd*a1
458  splasda = a2*splasda + a3*gl2da*a1
459  splasdz = a2*splasdz + a3*gl2dz*a1
460 
461  a2 = gl6
462  a3 = 3.0d0*gl6*gum
463 
464  a1 = splas
465  splas = a2*a1
466  splasdt = a2*splasdt + a3*gl2dt*a1
467  splasdd = a2*splasdd + a3*gl2dd*a1
468  splasda = a2*splasda + a3*gl2da*a1
469  splasdz = a2*splasdz + a3*gl2dz*a1
470 
471 
472  a2 = 0.93153d0 * 3.0d21 * xl9
473  a3 = 0.93153d0 * 3.0d21 * 9.0d0*xl8*xldt
474 
475  a1 = splas
476  splas = a2*a1
477  splasdt = a2*splasdt + a3*a1
478  splasdd = a2*splasdd
479  splasda = a2*splasda
480  splasdz = a2*splasdz
481 
482 
483 
484 
485 ! photoneutrino process section
486 ! for reactions like e- + gamma => e- + nu_e + nubar_e
487 ! e+ + gamma => e+ + nu_e + nubar_e
488 ! equation 3.8 for tau, equation 3.6 for cc,
489 ! and table 2 written out for speed
490  if (temp .ge. 1.0d7 .and. temp .lt. 1.0d8) then
491  tau = log10(temp * 1.0d-7)
492  cc = 0.5654d0 + tau
493  c00 = 1.008d11
494  c01 = 0.0d0
495  c02 = 0.0d0
496  c03 = 0.0d0
497  c04 = 0.0d0
498  c05 = 0.0d0
499  c06 = 0.0d0
500  c10 = 8.156d10
501  c11 = 9.728d8
502  c12 = -3.806d9
503  c13 = -4.384d9
504  c14 = -5.774d9
505  c15 = -5.249d9
506  c16 = -5.153d9
507  c20 = 1.067d11
508  c21 = -9.782d9
509  c22 = -7.193d9
510  c23 = -6.936d9
511  c24 = -6.893d9
512  c25 = -7.041d9
513  c26 = -7.193d9
514  dd01 = 0.0d0
515  dd02 = 0.0d0
516  dd03 = 0.0d0
517  dd04 = 0.0d0
518  dd05 = 0.0d0
519  dd11 = -1.879d10
520  dd12 = -9.667d9
521  dd13 = -5.602d9
522  dd14 = -3.370d9
523  dd15 = -1.825d9
524  dd21 = -2.919d10
525  dd22 = -1.185d10
526  dd23 = -7.270d9
527  dd24 = -4.222d9
528  dd25 = -1.560d9
529 
530  else if (temp .ge. 1.0d8 .and. temp .lt. 1.0d9) then
531  tau = log10(temp * 1.0d-8)
532  cc = 1.5654d0
533  c00 = 9.889d10
534  c01 = -4.524d8
535  c02 = -6.088d6
536  c03 = 4.269d7
537  c04 = 5.172d7
538  c05 = 4.910d7
539  c06 = 4.388d7
540  c10 = 1.813d11
541  c11 = -7.556d9
542  c12 = -3.304d9
543  c13 = -1.031d9
544  c14 = -1.764d9
545  c15 = -1.851d9
546  c16 = -1.928d9
547  c20 = 9.750d10
548  c21 = 3.484d10
549  c22 = 5.199d9
550  c23 = -1.695d9
551  c24 = -2.865d9
552  c25 = -3.395d9
553  c26 = -3.418d9
554  dd01 = -1.135d8
555  dd02 = 1.256d8
556  dd03 = 5.149d7
557  dd04 = 3.436d7
558  dd05 = 1.005d7
559  dd11 = 1.652d9
560  dd12 = -3.119d9
561  dd13 = -1.839d9
562  dd14 = -1.458d9
563  dd15 = -8.956d8
564  dd21 = -1.548d10
565  dd22 = -9.338d9
566  dd23 = -5.899d9
567  dd24 = -3.035d9
568  dd25 = -1.598d9
569 
570  else if (temp .ge. 1.0d9) then
571  tau = log10(t9)
572  cc = 1.5654d0
573  c00 = 9.581d10
574  c01 = 4.107d8
575  c02 = 2.305d8
576  c03 = 2.236d8
577  c04 = 1.580d8
578  c05 = 2.165d8
579  c06 = 1.721d8
580  c10 = 1.459d12
581  c11 = 1.314d11
582  c12 = -1.169d11
583  c13 = -1.765d11
584  c14 = -1.867d11
585  c15 = -1.983d11
586  c16 = -1.896d11
587  c20 = 2.424d11
588  c21 = -3.669d9
589  c22 = -8.691d9
590  c23 = -7.967d9
591  c24 = -7.932d9
592  c25 = -7.987d9
593  c26 = -8.333d9
594  dd01 = 4.724d8
595  dd02 = 2.976d8
596  dd03 = 2.242d8
597  dd04 = 7.937d7
598  dd05 = 4.859d7
599  dd11 = -7.094d11
600  dd12 = -3.697d11
601  dd13 = -2.189d11
602  dd14 = -1.273d11
603  dd15 = -5.705d10
604  dd21 = -2.254d10
605  dd22 = -1.551d10
606  dd23 = -7.793d9
607  dd24 = -4.489d9
608  dd25 = -2.185d9
609  end if
610 
611  taudt = iln10*tempi
612 
613 
614 ! equation 3.7, compute the expensive trig functions only one time
615  cos1 = cos(fac1*tau)
616  cos2 = cos(fac1*2.0d0*tau)
617  cos3 = cos(fac1*3.0d0*tau)
618  cos4 = cos(fac1*4.0d0*tau)
619  cos5 = cos(fac1*5.0d0*tau)
620  last = cos(fac2*tau)
621 
622  sin1 = sin(fac1*tau)
623  sin2 = sin(fac1*2.0d0*tau)
624  sin3 = sin(fac1*3.0d0*tau)
625  sin4 = sin(fac1*4.0d0*tau)
626  sin5 = sin(fac1*5.0d0*tau)
627  xast = sin(fac2*tau)
628 
629  a0 = 0.5d0*c00 &
630  + c01*cos1 + dd01*sin1 + c02*cos2 + dd02*sin2 &
631  + c03*cos3 + dd03*sin3 + c04*cos4 + dd04*sin4 &
632  + c05*cos5 + dd05*sin5 + 0.5d0*c06*last
633 
634  f0 = taudt*fac1*(-c01*sin1 + dd01*cos1 - c02*sin2*2.0d0 &
635  + dd02*cos2*2.0d0 - c03*sin3*3.0d0 + dd03*cos3*3.0d0 &
636  - c04*sin4*4.0d0 + dd04*cos4*4.0d0 &
637  - c05*sin5*5.0d0 + dd05*cos5*5.0d0) &
638  - 0.5d0*c06*xast*fac2*taudt
639 
640  a1 = 0.5d0*c10 &
641  + c11*cos1 + dd11*sin1 + c12*cos2 + dd12*sin2 &
642  + c13*cos3 + dd13*sin3 + c14*cos4 + dd14*sin4 &
643  + c15*cos5 + dd15*sin5 + 0.5d0*c16*last
644 
645  f1 = taudt*fac1*(-c11*sin1 + dd11*cos1 - c12*sin2*2.0d0 &
646  + dd12*cos2*2.0d0 - c13*sin3*3.0d0 + dd13*cos3*3.0d0 &
647  - c14*sin4*4.0d0 + dd14*cos4*4.0d0 - c15*sin5*5.0d0 &
648  + dd15*cos5*5.0d0) - 0.5d0*c16*xast*fac2*taudt
649 
650  a2 = 0.5d0*c20 &
651  + c21*cos1 + dd21*sin1 + c22*cos2 + dd22*sin2 &
652  + c23*cos3 + dd23*sin3 + c24*cos4 + dd24*sin4 &
653  + c25*cos5 + dd25*sin5 + 0.5d0*c26*last
654 
655  f2 = taudt*fac1*(-c21*sin1 + dd21*cos1 - c22*sin2*2.0d0 &
656  + dd22*cos2*2.0d0 - c23*sin3*3.0d0 + dd23*cos3*3.0d0 &
657  - c24*sin4*4.0d0 + dd24*cos4*4.0d0 - c25*sin5*5.0d0 &
658  + dd25*cos5*5.0d0) - 0.5d0*c26*xast*fac2*taudt
659 
660 ! equation 3.4
661  dum = a0 + a1*zeta + a2*zeta2
662  dumdt = f0 + f1*zeta + a1*zetadt + f2*zeta2 + 2.0d0*a2*zeta*zetadt
663  dumdd = a1*zetadd + 2.0d0*a2*zeta*zetadd
664  dumda = a1*zetada + 2.0d0*a2*zeta*zetada
665  dumdz = a1*zetadz + 2.0d0*a2*zeta*zetadz
666 
667  z = exp(-cc*zeta)
668 
669  xnum = dum*z
670  xnumdt = dumdt*z - dum*z*cc*zetadt
671  xnumdd = dumdd*z - dum*z*cc*zetadd
672  xnumda = dumda*z - dum*z*cc*zetada
673  xnumdz = dumdz*z - dum*z*cc*zetadz
674 
675  xden = zeta3 + 6.290d-3*xlm1 + 7.483d-3*xlm2 + 3.061d-4*xlm3
676 
677  dum = 3.0d0*zeta2
678  xdendt = dum*zetadt - xldt*(6.290d-3*xlm2 &
679  + 2.0d0*7.483d-3*xlm3 + 3.0d0*3.061d-4*xlm4)
680  xdendd = dum*zetadd
681  xdenda = dum*zetada
682  xdendz = dum*zetadz
683 
684  dum = 1.0d0/xden
685  fphot = xnum*dum
686  fphotdt = (xnumdt - fphot*xdendt)*dum
687  fphotdd = (xnumdd - fphot*xdendd)*dum
688  fphotda = (xnumda - fphot*xdenda)*dum
689  fphotdz = (xnumdz - fphot*xdendz)*dum
690 
691 
692 ! equation 3.3
693  a0 = 1.0d0 + 2.045d0 * xl
694  xnum = 0.666d0*a0**(-2.066d0)
695  xnumdt = -2.066d0*xnum/a0 * 2.045d0*xldt
696 
697  dum = 1.875d8*xl + 1.653d8*xl2 + 8.499d8*xl3 - 1.604d8*xl4
698  dumdt = xldt*(1.875d8 + 2.0d0*1.653d8*xl + 3.0d0*8.449d8*xl2 &
699  - 4.0d0*1.604d8*xl3)
700 
701  z = 1.0d0/dum
702  xden = 1.0d0 + rm*z
703  xdendt = -rm*z*z*dumdt
704  xdendd = rmdd*z
705  xdenda = rmda*z
706  xdendz = rmdz*z
707 
708  z = 1.0d0/xden
709  qphot = xnum*z
710  qphotdt = (xnumdt - qphot*xdendt)*z
711  dum = -qphot*z
712  qphotdd = dum*xdendd
713  qphotda = dum*xdenda
714  qphotdz = dum*xdendz
715 
716 ! equation 3.2
717  sphot = xl5 * fphot
718  sphotdt = 5.0d0*xl4*xldt*fphot + xl5*fphotdt
719  sphotdd = xl5*fphotdd
720  sphotda = xl5*fphotda
721  sphotdz = xl5*fphotdz
722 
723  a1 = sphot
724  sphot = rm*a1
725  sphotdt = rm*sphotdt
726  sphotdd = rm*sphotdd + rmdd*a1
727  sphotda = rm*sphotda + rmda*a1
728  sphotdz = rm*sphotdz + rmdz*a1
729 
730  a1 = tfac4*(1.0d0 - tfac3 * qphot)
731  a2 = -tfac4*tfac3
732 
733  a3 = sphot
734  sphot = a1*a3
735  sphotdt = a1*sphotdt + a2*qphotdt*a3
736  sphotdd = a1*sphotdd + a2*qphotdd*a3
737  sphotda = a1*sphotda + a2*qphotda*a3
738  sphotdz = a1*sphotdz + a2*qphotdz*a3
739 
740  if (sphot .le. 0.0) then
741  sphot = 0.0d0
742  sphotdt = 0.0d0
743  sphotdd = 0.0d0
744  sphotda = 0.0d0
745  sphotdz = 0.0d0
746  end if
747 
748 
749 
750 
751 
752 ! bremsstrahlung neutrino section
753 ! for reactions like e- + (z,a) => e- + (z,a) + nu + nubar
754 ! n + n => n + n + nu + nubar
755 ! n + p => n + p + nu + nubar
756 ! equation 4.3
757 
758  den6 = den * 1.0d-6
759  t8 = temp * 1.0d-8
760  t812 = sqrt(t8)
761  t832 = t8 * t812
762  t82 = t8*t8
763  t83 = t82*t8
764  t85 = t82*t83
765  t86 = t85*t8
766  t8m1 = 1.0d0/t8
767  t8m2 = t8m1*t8m1
768  t8m3 = t8m2*t8m1
769  t8m5 = t8m3*t8m2
770  t8m6 = t8m5*t8m1
771 
772 
773  tfermi = 5.9302d9*(sqrt(1.0d0+1.018d0*(den6*ye)**twoth)-1.0d0)
774 
775 ! "weak" degenerate electrons only
776  if (temp .gt. 0.3d0 * tfermi) then
777 
778 ! equation 5.3
779  dum = 7.05d6 * t832 + 5.12d4 * t83
780  dumdt = (1.5d0*7.05d6*t812 + 3.0d0*5.12d4*t82)*1.0d-8
781 
782  z = 1.0d0/dum
783  eta = rm*z
784  etadt = -rm*z*z*dumdt
785  etadd = rmdd*z
786  etada = rmda*z
787  etadz = rmdz*z
788 
789  etam1 = 1.0d0/eta
790  etam2 = etam1 * etam1
791  etam3 = etam2 * etam1
792 
793 
794 ! equation 5.2
795  a0 = 23.5d0 + 6.83d4*t8m2 + 7.81d8*t8m5
796  f0 = (-2.0d0*6.83d4*t8m3 - 5.0d0*7.81d8*t8m6)*1.0d-8
797  xnum = 1.0d0/a0
798 
799  dum = 1.0d0 + 1.47d0*etam1 + 3.29d-2*etam2
800  z = -1.47d0*etam2 - 2.0d0*3.29d-2*etam3
801  dumdt = z*etadt
802  dumdd = z*etadd
803  dumda = z*etada
804  dumdz = z*etadz
805 
806  c00 = 1.26d0 * (1.0d0+etam1)
807  z = -1.26d0*etam2
808  c01 = z*etadt
809  c02 = z*etadd
810  c03 = z*etada
811  c04 = z*etadz
812 
813  z = 1.0d0/dum
814  xden = c00*z
815  xdendt = (c01 - xden*dumdt)*z
816  xdendd = (c02 - xden*dumdd)*z
817  xdenda = (c03 - xden*dumda)*z
818  xdendz = (c04 - xden*dumdz)*z
819 
820  fbrem = xnum + xden
821  fbremdt = -xnum*xnum*f0 + xdendt
822  fbremdd = xdendd
823  fbremda = xdenda
824  fbremdz = xdendz
825 
826 
827 ! equation 5.9
828  a0 = 230.0d0 + 6.7d5*t8m2 + 7.66d9*t8m5
829  f0 = (-2.0d0*6.7d5*t8m3 - 5.0d0*7.66d9*t8m6)*1.0d-8
830 
831  z = 1.0d0 + rm*1.0d-9
832  dum = a0*z
833  dumdt = f0*z
834  z = a0*1.0d-9
835  dumdd = z*rmdd
836  dumda = z*rmda
837  dumdz = z*rmdz
838 
839  xnum = 1.0d0/dum
840  z = -xnum*xnum
841  xnumdt = z*dumdt
842  xnumdd = z*dumdd
843  xnumda = z*dumda
844  xnumdz = z*dumdz
845 
846  c00 = 7.75d5*t832 + 247.0d0*t8**(3.85d0)
847  dd00 = (1.5d0*7.75d5*t812 + 3.85d0*247.0d0*t8**(2.85d0))*1.0d-8
848 
849  c01 = 4.07d0 + 0.0240d0 * t8**(1.4d0)
850  dd01 = 1.4d0*0.0240d0*t8**(0.4d0)*1.0d-8
851 
852  c02 = 4.59d-5 * t8**(-0.110d0)
853  dd02 = -0.11d0*4.59d-5 * t8**(-1.11d0)*1.0d-8
854 
855  z = den**(0.656d0)
856  dum = c00*rmi + c01 + c02*z
857  dumdt = dd00*rmi + dd01 + dd02*z
858  z = -c00*rmi*rmi
859  dumdd = z*rmdd + 0.656d0*c02*den**(-0.454d0)
860  dumda = z*rmda
861  dumdz = z*rmdz
862 
863  xden = 1.0d0/dum
864  z = -xden*xden
865  xdendt = z*dumdt
866  xdendd = z*dumdd
867  xdenda = z*dumda
868  xdendz = z*dumdz
869 
870  gbrem = xnum + xden
871  gbremdt = xnumdt + xdendt
872  gbremdd = xnumdd + xdendd
873  gbremda = xnumda + xdenda
874  gbremdz = xnumdz + xdendz
875 
876 
877 ! equation 5.1
878  dum = 0.5738d0*zbar*ye*t86*den
879  dumdt = 0.5738d0*zbar*ye*6.0d0*t85*den*1.0d-8
880  dumdd = 0.5738d0*zbar*ye*t86
881  dumda = -dum*abari
882  dumdz = 0.5738d0*2.0d0*ye*t86*den
883 
884  z = tfac4*fbrem - tfac5*gbrem
885  sbrem = dum * z
886  sbremdt = dumdt*z + dum*(tfac4*fbremdt - tfac5*gbremdt)
887  sbremdd = dumdd*z + dum*(tfac4*fbremdd - tfac5*gbremdd)
888  sbremda = dumda*z + dum*(tfac4*fbremda - tfac5*gbremda)
889  sbremdz = dumdz*z + dum*(tfac4*fbremdz - tfac5*gbremdz)
890 
891 
892 
893 
894 ! liquid metal with c12 parameters (not too different for other elements)
895 ! equation 5.18 and 5.16
896 
897  else
898  u = fac3 * (log10(den) - 3.0d0)
899  a0 = iln10*fac3*deni
900 
901 ! compute the expensive trig functions of equation 5.21 only once
902  cos1 = cos(u)
903  cos2 = cos(2.0d0*u)
904  cos3 = cos(3.0d0*u)
905  cos4 = cos(4.0d0*u)
906  cos5 = cos(5.0d0*u)
907 
908  sin1 = sin(u)
909  sin2 = sin(2.0d0*u)
910  sin3 = sin(3.0d0*u)
911  sin4 = sin(4.0d0*u)
912  sin5 = sin(5.0d0*u)
913 
914 ! equation 5.21
915  fb = 0.5d0 * 0.17946d0 + 0.00945d0*u + 0.34529d0 &
916  - 0.05821d0*cos1 - 0.04969d0*sin1 &
917  - 0.01089d0*cos2 - 0.01584d0*sin2 &
918  - 0.01147d0*cos3 - 0.00504d0*sin3 &
919  - 0.00656d0*cos4 - 0.00281d0*sin4 &
920  - 0.00519d0*cos5
921 
922  c00 = a0*(0.00945d0 &
923  + 0.05821d0*sin1 - 0.04969d0*cos1 &
924  + 0.01089d0*sin2*2.0d0 - 0.01584d0*cos2*2.0d0 &
925  + 0.01147d0*sin3*3.0d0 - 0.00504d0*cos3*3.0d0 &
926  + 0.00656d0*sin4*4.0d0 - 0.00281d0*cos4*4.0d0 &
927  + 0.00519d0*sin5*5.0d0)
928 
929 
930 ! equation 5.22
931  ft = 0.5d0 * 0.06781d0 - 0.02342d0*u + 0.24819d0 &
932  - 0.00944d0*cos1 - 0.02213d0*sin1 &
933  - 0.01289d0*cos2 - 0.01136d0*sin2 &
934  - 0.00589d0*cos3 - 0.00467d0*sin3 &
935  - 0.00404d0*cos4 - 0.00131d0*sin4 &
936  - 0.00330d0*cos5
937 
938  c01 = a0*(-0.02342d0 &
939  + 0.00944d0*sin1 - 0.02213d0*cos1 &
940  + 0.01289d0*sin2*2.0d0 - 0.01136d0*cos2*2.0d0 &
941  + 0.00589d0*sin3*3.0d0 - 0.00467d0*cos3*3.0d0 &
942  + 0.00404d0*sin4*4.0d0 - 0.00131d0*cos4*4.0d0 &
943  + 0.00330d0*sin5*5.0d0)
944 
945 
946 ! equation 5.23
947  gb = 0.5d0 * 0.00766d0 - 0.01259d0*u + 0.07917d0 &
948  - 0.00710d0*cos1 + 0.02300d0*sin1 &
949  - 0.00028d0*cos2 - 0.01078d0*sin2 &
950  + 0.00232d0*cos3 + 0.00118d0*sin3 &
951  + 0.00044d0*cos4 - 0.00089d0*sin4 &
952  + 0.00158d0*cos5
953 
954  c02 = a0*(-0.01259d0 &
955  + 0.00710d0*sin1 + 0.02300d0*cos1 &
956  + 0.00028d0*sin2*2.0d0 - 0.01078d0*cos2*2.0d0 &
957  - 0.00232d0*sin3*3.0d0 + 0.00118d0*cos3*3.0d0 &
958  - 0.00044d0*sin4*4.0d0 - 0.00089d0*cos4*4.0d0 &
959  - 0.00158d0*sin5*5.0d0)
960 
961 
962 ! equation 5.24
963  gt = -0.5d0 * 0.00769d0 - 0.00829d0*u + 0.05211d0 &
964  + 0.00356d0*cos1 + 0.01052d0*sin1 &
965  - 0.00184d0*cos2 - 0.00354d0*sin2 &
966  + 0.00146d0*cos3 - 0.00014d0*sin3 &
967  + 0.00031d0*cos4 - 0.00018d0*sin4 &
968  + 0.00069d0*cos5
969 
970  c03 = a0*(-0.00829d0 &
971  - 0.00356d0*sin1 + 0.01052d0*cos1 &
972  + 0.00184d0*sin2*2.0d0 - 0.00354d0*cos2*2.0d0 &
973  - 0.00146d0*sin3*3.0d0 - 0.00014d0*cos3*3.0d0 &
974  - 0.00031d0*sin4*4.0d0 - 0.00018d0*cos4*4.0d0 &
975  - 0.00069d0*sin5*5.0d0)
976 
977 
978  dum = 2.275d-1 * zbar * zbar*t8m1 * (den6*abari)**oneth
979  dumdt = -dum*tempi
980  dumdd = oneth*dum*deni
981  dumda = -oneth*dum*abari
982  dumdz = 2.0d0*dum*zbari
983 
984  gm1 = 1.0d0/dum
985  gm2 = gm1*gm1
986  gm13 = gm1**oneth
987  gm23 = gm13 * gm13
988  gm43 = gm13*gm1
989  gm53 = gm23*gm1
990 
991 
992 ! equation 5.25 and 5.26
993  v = -0.05483d0 - 0.01946d0*gm13 + 1.86310d0*gm23 - 0.78873d0*gm1
994  a0 = oneth*0.01946d0*gm43 - twoth*1.86310d0*gm53 + 0.78873d0*gm2
995 
996  w = -0.06711d0 + 0.06859d0*gm13 + 1.74360d0*gm23 - 0.74498d0*gm1
997  a1 = -oneth*0.06859d0*gm43 - twoth*1.74360d0*gm53 + 0.74498d0*gm2
998 
999 
1000 ! equation 5.19 and 5.20
1001  fliq = v*fb + (1.0d0 - v)*ft
1002  fliqdt = a0*dumdt*(fb - ft)
1003  fliqdd = a0*dumdd*(fb - ft) + v*c00 + (1.0d0 - v)*c01
1004  fliqda = a0*dumda*(fb - ft)
1005  fliqdz = a0*dumdz*(fb - ft)
1006 
1007  gliq = w*gb + (1.0d0 - w)*gt
1008  gliqdt = a1*dumdt*(gb - gt)
1009  gliqdd = a1*dumdd*(gb - gt) + w*c02 + (1.0d0 - w)*c03
1010  gliqda = a1*dumda*(gb - gt)
1011  gliqdz = a1*dumdz*(gb - gt)
1012 
1013 
1014 ! equation 5.17
1015  dum = 0.5738d0*zbar*ye*t86*den
1016  dumdt = 0.5738d0*zbar*ye*6.0d0*t85*den*1.0d-8
1017  dumdd = 0.5738d0*zbar*ye*t86
1018  dumda = -dum*abari
1019  dumdz = 0.5738d0*2.0d0*ye*t86*den
1020 
1021  z = tfac4*fliq - tfac5*gliq
1022  sbrem = dum * z
1023  sbremdt = dumdt*z + dum*(tfac4*fliqdt - tfac5*gliqdt)
1024  sbremdd = dumdd*z + dum*(tfac4*fliqdd - tfac5*gliqdd)
1025  sbremda = dumda*z + dum*(tfac4*fliqda - tfac5*gliqda)
1026  sbremdz = dumdz*z + dum*(tfac4*fliqdz - tfac5*gliqdz)
1027 
1028  end if
1029 
1030 
1031 
1032 
1033 ! recombination neutrino section
1034 ! for reactions like e- (continuum) => e- (bound) + nu_e + nubar_e
1035 ! equation 6.11 solved for nu
1036  xnum = 1.10520d8 * den * ye /(temp*sqrt(temp))
1037  xnumdt = -1.50d0*xnum*tempi
1038  xnumdd = xnum*deni
1039  xnumda = -xnum*abari
1040  xnumdz = xnum*zbari
1041 
1042 ! the chemical potential
1043  nu = ifermi12(xnum)
1044 
1045 ! a0 is d(nu)/d(xnum)
1046  a0 = 1.0d0/(0.5d0*zfermim12(nu))
1047  nudt = a0*xnumdt
1048  nudd = a0*xnumdd
1049  nuda = a0*xnumda
1050  nudz = a0*xnumdz
1051 
1052  nu2 = nu * nu
1053  nu3 = nu2 * nu
1054 
1055 ! table 12
1056  if (nu .ge. -20.0 .and. nu .lt. 0.0) then
1057  a1 = 1.51d-2
1058  a2 = 2.42d-1
1059  a3 = 1.21d0
1060  b = 3.71d-2
1061  c = 9.06e-1
1062  d = 9.28d-1
1063  f1 = 0.0d0
1064  f2 = 0.0d0
1065  f3 = 0.0d0
1066  else if (nu .ge. 0.0 .and. nu .le. 10.0) then
1067  a1 = 1.23d-2
1068  a2 = 2.66d-1
1069  a3 = 1.30d0
1070  b = 1.17d-1
1071  c = 8.97e-1
1072  d = 1.77d-1
1073  f1 = -1.20d-2
1074  f2 = 2.29d-2
1075  f3 = -1.04d-3
1076  end if
1077 
1078 
1079 ! equation 6.7, 6.13 and 6.14
1080  if (nu .ge. -20.0 .and. nu .le. 10.0) then
1081 
1082  zeta = 1.579d5*zbar*zbar*tempi
1083  zetadt = -zeta*tempi
1084  zetadd = 0.0d0
1085  zetada = 0.0d0
1086  zetadz = 2.0d0*zeta*zbari
1087 
1088  c00 = 1.0d0/(1.0d0 + f1*nu + f2*nu2 + f3*nu3)
1089  c01 = f1 + f2*2.0d0*nu + f3*3.0d0*nu2
1090  dum = zeta*c00
1091  dumdt = zetadt*c00 + zeta*c01*nudt
1092  dumdd = zeta*c01*nudd
1093  dumda = zeta*c01*nuda
1094  dumdz = zetadz*c00 + zeta*c01*nudz
1095 
1096 
1097  z = 1.0d0/dum
1098  dd00 = dum**(-2.25)
1099  dd01 = dum**(-4.55)
1100  c00 = a1*z + a2*dd00 + a3*dd01
1101  c01 = -(a1*z + 2.25*a2*dd00 + 4.55*a3*dd01)*z
1102 
1103 
1104  z = exp(c*nu)
1105  dd00 = b*z*(1.0d0 + d*dum)
1106  gum = 1.0d0 + dd00
1107  gumdt = dd00*c*nudt + b*z*d*dumdt
1108  gumdd = dd00*c*nudd + b*z*d*dumdd
1109  gumda = dd00*c*nuda + b*z*d*dumda
1110  gumdz = dd00*c*nudz + b*z*d*dumdz
1111 
1112 
1113  z = exp(nu)
1114  a1 = 1.0d0/gum
1115 
1116  bigj = c00 * z * a1
1117  bigjdt = c01*dumdt*z*a1 + c00*z*nudt*a1 - c00*z*a1*a1 * gumdt
1118  bigjdd = c01*dumdd*z*a1 + c00*z*nudd*a1 - c00*z*a1*a1 * gumdd
1119  bigjda = c01*dumda*z*a1 + c00*z*nuda*a1 - c00*z*a1*a1 * gumda
1120  bigjdz = c01*dumdz*z*a1 + c00*z*nudz*a1 - c00*z*a1*a1 * gumdz
1121 
1122 
1123 ! equation 6.5
1124  z = exp(zeta + nu)
1125  dum = 1.0d0 + z
1126  a1 = 1.0d0/dum
1127  a2 = 1.0d0/bigj
1128 
1129  sreco = tfac6 * 2.649d-18 * ye * zbar**13 * den * bigj*a1
1130  srecodt = sreco*(bigjdt*a2 - z*(zetadt + nudt)*a1)
1131  srecodd = sreco*(1.0d0*deni + bigjdd*a2 - z*(zetadd + nudd)*a1)
1132  srecoda = sreco*(-1.0d0*abari + bigjda*a2 - z*(zetada+nuda)*a1)
1133  srecodz = sreco*(14.0d0*zbari + bigjdz*a2 - z*(zetadz+nudz)*a1)
1134 
1135  end if
1136 
1137 
1138 ! convert from erg/cm^3/s to erg/g/s
1139 ! comment these out to duplicate the itoh et al plots
1140 
1141  spair = spair*deni
1142  spairdt = spairdt*deni
1143  spairdd = spairdd*deni - spair*deni
1144  spairda = spairda*deni
1145  spairdz = spairdz*deni
1146 
1147  splas = splas*deni
1148  splasdt = splasdt*deni
1149  splasdd = splasdd*deni - splas*deni
1150  splasda = splasda*deni
1151  splasdz = splasdz*deni
1152 
1153  sphot = sphot*deni
1154  sphotdt = sphotdt*deni
1155  sphotdd = sphotdd*deni - sphot*deni
1156  sphotda = sphotda*deni
1157  sphotdz = sphotdz*deni
1158 
1159  sbrem = sbrem*deni
1160  sbremdt = sbremdt*deni
1161  sbremdd = sbremdd*deni - sbrem*deni
1162  sbremda = sbremda*deni
1163  sbremdz = sbremdz*deni
1164 
1165  sreco = sreco*deni
1166  srecodt = srecodt*deni
1167  srecodd = srecodd*deni - sreco*deni
1168  srecoda = srecoda*deni
1169  srecodz = srecodz*deni
1170 
1171 
1172 ! the total neutrino loss rate
1173  snu = splas + spair + sphot + sbrem + sreco
1174  dsnudt = splasdt + spairdt + sphotdt + sbremdt + srecodt
1175  dsnudd = splasdd + spairdd + sphotdd + sbremdd + srecodd
1176  dsnuda = splasda + spairda + sphotda + sbremda + srecoda
1177  dsnudz = splasdz + spairdz + sphotdz + sbremdz + srecodz
1178 
1179  return
1180  end
sneut5_aa
subroutine sneut5_aa(temp, den, abar, zbar, snu, dsnudt, dsnudd, dsnuda, dsnudz, spair, splas, sphot, sbrem, sreco)
Definition: sneut5.f90:4
pi
real(r_kind) function pi()
Further information: http://netlib.org/quadpack/index.html https://orion.math.iastate....
Definition: quadpack_module.f90:193