PROGRAM NEUT IMPLICIT REAL*8 (A-H,O-Z) RMU=2.0 10 WRITE (6,11) 11 FORMAT (' please input t9 and rho. Y_e is assumed 0.5',/) READ (5,*) T9,RHO IF (T9.LT.0.) GO TO 20 zbar=6.0d0 abar=12.0d0 temp=1.0d9*T9 call sneut (temp,rho,zbar,abar,pair,plas,phot,Q) write (6,13) 13 format (' neutrino energy loss rates in erg/gm/s',/) WRITE (6,12) T9,RHO,Q 12 FORMAT (' T9 = ',F6.2,3X,' RHO =',1PE12.2,' TOTAL = ',E10.2) WRITE (6,15) PLAS,PAIR,PHOT 15 FORMAT (' PLAS = ',1PE10.2,' PAIR = ',E10.2,' PHOT = ', 1 E10.2,/) GO TO 10 20 CALL EXIT END CXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX subroutine sneut(temp,den,zbar,abar,pair,plas,phot,total) implicit none save c.. c..this routine computes neutrino losses from the analytic fits of c..itoh et al. apjs 102, 411, 1996 c... this version for 220c ignores bremms and recomb neutrinos c.. c..input is the temperature temp, density den, c..mean number of nucleons abar, and mean charge zbar. c.. c..output is the c..pair neutrino contributions pair, c..plasma neutrino contributions plas, c..photoneutrino contributions phot, c..and their sum total, c..all in erg/g/s. c.. c..declare REAL*8 temp,den,zbar,abar, 1 pair,plas,phot,brem,recomb,total INTEGER iold REAL*8 xmue,t9,xl,xlp5,xl2,xl3,xl4,xl5,xl6,xl8,xl9, 1 xlm1,xlm2,xlm3,rm,xnum,xden,fpair,fphoto,gl, 2 zeta,zeta2,zeta3,qpair,gl2,gl12,gl32,gl72,gl6, 3 ft,fl,cc,xlnt,fxy,qv,tau,qphoto,den6,tfermi,t8, 4 t832,t83,t86,t8m2,t8m5,eta,etam1,etam2,fbrem, 5 gbrem,a0,a1,a2,a3,b1,b2,b3,c00, 6 c01,c02,c03,c04,c05,c06,c10,c11,c12,c13,c14,c15, 7 c16,c20,c21,c22,c23,c24,c25,c26,dd01,dd02,dd03, 8 dd04,dd05,dd11,dd12,dd13,dd14,dd15,dd21,dd22, 9 dd23,dd24,dd25,u,gamma,gm1,gm13,gm23,v,w,fb, 1 gb,gt,fliq,gliq,ifermi12,nu,nu2,nu3,b,c,d,f1, 2 f2,f3,z,bigj,cos1,cos2,cos3,cos4,cos5,sin1,sin2, 3 sin3,sin4,sin5,ye,last,deni REAL*8 pi,fac1,fac2,fac3,third,twoth,con1,sixth parameter (pi = 3.1415926535897932384d0, 1 fac1 = 5.0d0 * pi / 3.0d0, 2 fac2 = 10.0d0 * pi, 3 fac3 = pi / 5.0d0, 4 third = 1.0d0/3.0d0, 5 twoth = 2.0d0/3.0d0, 6 con1 = 1.0d0/5.9302d0, 7 sixth = 1.0d0/6.0d0) c.. theta is sin**2(theta_weinberg) = 0.2319 plus/minus 0.00005 (1996) c.. xnufam is the number of neutrino flavors = 3.02 plus/minus 0.005 (1998) c.. change theta and xnufam if need be, and the changes will automatically c.. propagate through the routine. cv and ca are the vector and axial currents. REAL*8 theta,xnufam,cv,ca,cvp,cap,tfac1,tfac2,tfac3, 1 tfac4,tfac5,tfac6 parameter (theta = 0.2319d0, 1 xnufam = 3.0d0, 2 cv = 0.5d0 + 2.0d0 * theta, 3 cvp = 1.0d0 - cv, 4 ca = 0.5d0, 5 cap = 1.0d0 - ca, 6 tfac1 = cv*cv + ca*ca + 7 (xnufam-1.0d0) * (cvp*cvp+cap*cap), 8 tfac2 = cv*cv - ca*ca + 9 (xnufam-1.0d0) * (cvp*cvp - cap*cap), & tfac3 = tfac2/tfac1, 1 tfac4 = 0.5d0 * tfac1, 2 tfac5 = 0.5d0 * tfac2, 3 tfac6 = cv*cv + 1.5d0*ca*ca + (xnufam - 1.0d0)* 4 (cvp*cvp + 1.5d0*cap*cap)) c.. variable for old neutrino loss REAL*8 t9o,rho,rmu,q c.. initialize and bail if its too cold pair = 0.0d0 plas = 0.0d0 phot = 0.0d0 brem = 0.0d0 recomb = 0.0d0 total = 0.0d0 if (temp .lt. 1.0e7) return c.. some frequent factors xmue = abar/zbar ye = zbar/abar t9 = temp * 1.0d-9 xl = t9 * con1 xlp5 = sqrt(xl) xl2 = xl*xl xl3 = xl2*xl xl4 = xl2*xl2 xl5 = xl*xl4 xl6 = xl2*xl4 xl8 = xl4*xl4 xl9 = xl8 * xl xlm1 = 1.0d0/xl xlm2 = xlm1*xlm1 xlm3 = xlm1*xlm2 rm = den/xmue zeta = (rm * 1.0d-9)**third * xlm1 zeta2 = zeta * zeta zeta3 = zeta2 * zeta c..pair neutrino section c..for reactions like e+ + e- => nu_e + nubar_e c..equation 2.8 gl = 1.0d0-13.04d0*xl 2 + 133.5d0*xl4 +1534.0d0*xl6 + 918.6d0*xl8 c..equation 2.7 if (t9 .lt. 10.0) then xnum = (6.002d19 + 2.084d20*zeta + 1.872d21*zeta2) 1 * exp(-5.5924d0*zeta) xden = zeta3 + 9.383d-1*xlm1 - 4.141d-1*xlm2 + 5.829d-2*xlm3 else xnum = (6.002d19 + 2.084d20*zeta + 1.872d21*zeta2) 1 * exp(-4.9924d0*zeta) xden = zeta3 + 1.2383d0*xlm1 - 8.141d-1*xlm2 end if fpair = xnum/xden c..equation 2.6 xnum = 1.0d0 / (10.7480d0*xl2 + 0.3967d0*xlp5 + 1.005d0) xden = (1.0d0 + rm/(7.692d7*xl3 + 9.715d6*xlp5) )**(-0.3d0) qpair = xnum*xden c..equation 2.5 pair = tfac4 * (1.0d0 + tfac3 * qpair) 1 * gl * exp(-2.0d0*xlm1) * fpair c..plasma neutrino section c..for collective reactions like gamma_plasmon => nu_e + nubar_e c..equation 4.6 gl2 = 1.1095d11*rm/(temp*temp * sqrt(1.0d0+(1.019d-6*rm)**twoth)) gl = sqrt(gl2) gl12 = sqrt(gl) gl32 = gl * gl12 gl72 = gl2 * gl32 gl6 = gl2 * gl2 * gl2 c.. equation 4.7 ft = 2.4d0 + 0.6d0*gl12 + 0.51d0*gl + 1.25d0*gl32 c.. equation 4.8 fl = (8.6d0*gl2 + 1.35d0*gl72)/(225.0d0 - 17.0d0*gl + gl2) c.. equation 4.9 and 4.10 cc = log10(2.0d0*rm) xlnt = log10(temp) xnum = sixth * ( 17.5d0 + cc - 3.0d0 * xlnt) xden = sixth * (-24.5d0 + cc + 3.0d0 * xlnt) c..equation 4.11 if (abs(xnum) .gt. 0.7d0 .or. xden .lt. 0.0d0) then fxy = 1.0d0 else fxy = 1.05d0 + (0.39d0 - 1.25d0*xnum - 0.35d0*sin(4.5d0*xnum) 1 - 0.3d0 * exp(-1.0d0*(4.5d0*xnum + 0.9d0)**2)) 2 * exp(-1.0d0* ( min(0.0d0, xden - 1.6d0 + 1.25d0*xnum) 3 / (0.57d0 - 0.25d0*xnum) )**2 ) end if c.. equation 4.5 qv = 3.0d21 * xl9 * gl6 * exp(-gl) * (ft + fl) * fxy c.. equation 4.1 plas = 0.93153d0 * qv c.. photoneutrino process section c.. for reactions like e- + gamma => e- + nu_e + nubar_e c.. e+ + gamma => e+ + nu_e + nubar_e c.. equation 3.8 for tau, equation 3.6 for cc, c.. and table 2 written out for speed if (temp .ge. 1.0d7 .and. temp .lt. 1.0d8) then tau = log10(temp * 1.0d-7) cc = 0.5654d0 + tau c00 = 1.008d11 c01 = 0.0d0 c02 = 0.0d0 c03 = 0.0d0 c04 = 0.0d0 c05 = 0.0d0 c06 = 0.0d0 c10 = 8.156d10 c11 = 9.728d8 c12 = -3.806d9 c13 = -4.384d9 c14 = -5.774d9 c15 = -5.249d9 c16 = -5.153d9 c20 = 1.067d11 c21 = -9.782d9 c22 = -7.193d9 c23 = -6.936d9 c24 = -6.893d9 c25 = -7.041d9 c26 = -7.193d9 dd01 = 0.0d0 dd02 = 0.0d0 dd03 = 0.0d0 dd04 = 0.0d0 dd05 = 0.0d0 dd11 = -1.879d10 dd12 = -9.667d9 dd13 = -5.602d9 dd14 = -3.370d9 dd15 = -1.825d9 dd21 = -2.919d10 dd22 = -1.185d10 dd23 = -7.270d9 dd24 = -4.222d9 dd25 = -1.560d9 else if (temp .ge. 1.0d8 .and. temp .lt. 1.0d9) then tau = log10(temp * 1.0d-8) cc = 1.5654d0 c00 = 9.889d10 c01 = -4.524d8 c02 = -6.088d6 c03 = 4.269d7 c04 = 5.172d7 c05 = 4.910d7 c06 = 4.388d7 c10 = 1.813d11 c11 = -7.556d9 c12 = -3.304d9 c13 = -1.031d9 c14 = -1.764d9 c15 = -1.851d9 c16 = -1.928d9 c20 = 9.750d10 c21 = 3.484d10 c22 = 5.199d9 c23 = -1.695d9 c24 = -2.865d9 c25 = -3.395d9 c26 = -3.418d9 dd01 = -1.135d8 dd02 = 1.256d8 dd03 = 5.149d7 dd04 = 3.436d7 dd05 = 1.005d7 dd11 = 1.652d9 dd12 = -3.119d9 dd13 = -1.839d9 dd14 = -1.458d9 dd15 = -8.956d8 dd21 = -1.549d10 dd22 = -9.338d9 dd23 = -5.899d9 dd24 = -3.035d9 dd25 = -1.598d9 else if (temp .ge. 1.0d9) then tau = log10(t9) cc = 1.5654d0 c00 = 9.581d10 c01 = 4.107d8 c02 = 2.305d8 c03 = 2.236d8 c04 = 1.580d8 c05 = 2.165d8 c06 = 1.721d8 c10 = 1.459d12 c11 = 1.314d11 c12 = -1.169d11 c13 = -1.765d11 c14 = -1.867d11 c15 = -1.983d11 c16 = -1.896d11 c20 = 2.424d11 c21 = -3.669d9 c22 = -8.691d9 c23 = -7.967d9 c24 = -7.932d9 c25 = -7.987d9 c26 = -8.333d9 dd01 = 4.724d8 dd02 = 2.976d8 dd03 = 2.242d8 dd04 = 7.937d7 dd05 = 4.859d7 dd11 = -7.094d11 dd12 = -3.697d11 dd13 = -2.189d11 dd14 = -1.273d11 dd15 = -5.705d10 dd21 = -2.254d10 dd22 = -1.551d10 dd23 = -7.793d9 dd24 = -4.489d9 dd25 = -2.185d9 end if c.. equation 3.7, compute the expensive trig functions only one time cos1 = cos(fac1*tau) cos2 = cos(fac1*2.0d0*tau) cos3 = cos(fac1*3.0d0*tau) cos4 = cos(fac1*4.0d0*tau) cos5 = cos(fac1*5.0d0*tau) last = cos(fac2*tau) sin1 = sin(fac1*tau) sin2 = sin(fac1*2.0d0*tau) sin3 = sin(fac1*3.0d0*tau) sin4 = sin(fac1*4.0d0*tau) sin5 = sin(fac1*5.0d0*tau) a0 = 0.5d0*c00 1 + c01*cos1 + dd01*sin1 + c02*cos2 + dd02*sin2 2 + c03*cos3 + dd03*sin3 + c04*cos4 + dd04*sin4 3 + c05*cos5 + dd05*sin5 + 0.5d0*c06*last a1 = 0.5d0*c10 1 + c11*cos1 + dd11*sin1 + c12*cos2 + dd12*sin2 2 + c13*cos3 + dd13*sin3 + c14*cos4 + dd14*sin4 3 + c15*cos5 + dd15*sin5 + 0.5d0*c16*last a2 = 0.5d0*c20 1 + c21*cos1 + dd21*sin1 + c22*cos2 + dd22*sin2 2 + c23*cos3 + dd23*sin3 + c24*cos4 + dd24*sin4 3 + c25*cos5 + dd25*sin5 + 0.5d0*c26*last c..equation 3.4 xnum = (a0 + a1*zeta + a2*zeta2)*exp(-cc*zeta) xden = zeta3 + 6.290d-3*xlm1 + 7.483d-3*xlm2 + 3.061d-4*xlm3 fphoto = xnum/xden c..equation 3.3 xnum = 0.666d0*((1.0d0 + 2.045d0 * xl)**(-2.066d0)) xden = 1.0d0 + rm / (1.875d8*xl + 1.653d8*xl2 1 + 8.449d8*xl3 - 1.604d8*xl4) qphoto = xnum/xden c..equation 3.2 phot = tfac4 * (1.0d0 - tfac3 * qphoto) * rm * xl5 * fphoto c.. convert from erg/cm^3/s to erg/g/s c.. comment these out to duplicate the itoh et al plots deni = 1.0d0/den pair = pair*deni plas = plas*deni phot = phot*deni c..the total neutrino loss rate total = plas + pair + phot end