subroutine gnutrino(u2,t2) * * Neutrino rates from Itoh, Adachi, Nakagawa, Kohyama, & Munkata, 1989, * ApJ, 339, 354. (+Errata 1990, ApJ, 360, 741) * (Plasma, Photo, Pair) * Recombination neutrinos from Kohyama, Itoh, Obama, & Mutoh 1993, ApJ, * 415, 267 * Bremsstrahlung from: * 1. Liquid Metal: Itoh & Kohyama, 1983, ApJ, 275, 858. * 2. Crystal Lattice: Itoh etal., 1984, ApJ, 279, 413. * 3. Low T Quantum Corr.: Itoh etal., 1984, ApJ, 280, 787. * 4. Phonon Contr.: Itoh etal., 1984, ApJ, 285, 304. * 5. Partially Degen. e: Munkata, Kohyama & Itoh, 1987, ApJ, 316, 708. * Prepared for WDXD by PAB: Oct. 1990 * Last messed with on 08 17 94 * implicit double precision (a-h,o-z) common/neut/en(5),fn common/temp/s1,r1,s2,r2,b2,p2,tx2,e2,xc2,xo2,fca2,f2,q2,w2,c common/je/je common/comp/amhyhe,amheca,x(4) common/c/az(3),z(3),taux,ts,ps,sm,rs,err,te,gs,rho,ol1x,rl1x,el1x 1,atgx,rtgx,tl,pl,xmass,acu,am,bk,an0,pi common/rhomooe/rhomooe common/phys/flam,gamma,degtemp,tfermi double precision mue data gamc,gamo,gamhe,gamh/3577200.,5778200.,573264.,227500./ * * initialize * t=10.d0**t2 rho=10.d0**u2 do 10 i=1,5 en(i) = 0.0 10 continue fn = 0.0 qbrem = 0.0 qplas = 0.0 qphot = 0.0 qpair = 0.0 qrec = 0.0 * * Escape clause (enter when T < 10^7K or Rho < 1 g/cc) * With this, I don't need an escape clause for flam * if(t2.lt.7.0 .or. u2.lt.0.0)then return endif * * compute gamma and mue for mixture * if(je.eq.1)then * * assume pure H for now. (Note: for H, I shouldn't be here anyway!) * ce = 1. mue = 1. gamma = gamh * 10.**(u2/3. - t2) zbar = 1.0 aavg = 1.0 else if(je.eq.2)then * * pure helium * omue = 1./(x(2)/2.) gamma = gamhe * 10.**(u2/3. - t2) zbar = 2.0 aavg = 4.0 else * * helium, carbon, oxygen, and iron (the latter for tests) * gammahe = gamhe * 10.**(u2/3. - t2) gammac = gamc * 10.**(u2/3. - t2) gammao = gamo * 10.**(u2/3. - t2) gamma = xc2*gammac + xo2*gammao omue = xc2*6./12. + xo2*8./16. zbar = xc2*6.0 + xo2*8.0 aavg = xc2*12.0 + xo2*16.0 endif endif * * zbar is the mean charge on a nucleus, for bremm and pair. * flam was x in neu * flam=t/5.9302d+9 mue = 1.d0 / omue rhomooe = rho/mue rhomue23 = (rhomooe)**(2./3.) sqfac = 1. + 1.018*rhomue23 tfermi = 5.9302d+9*(dsqrt(sqfac) - 1.) degtemp = t/tfermi y = (rhomooe)**(1./3.)/(1000.*flam) s2theta = 0.217 cv=0.5d0+2.0d0*s2theta cvp=1.d0-cv ca=0.5d0 cap=1.d0-ca if(y .le. 100.)then call plasmon(flam,rho,t2,mue,cv,cvp,ca,cap,qplas) endif if(y .le. 100.)then call photon(flam,t2,rho,mue,cv,cvp,ca,cap,qphot) endif if(flam.ge.0.10)then call pairn(flam,rho,mue,cv,cvp,ca,cap,qpair) endif call recombine(aavg,zbar,t,rho,cv,cvp,ca,cap,qrec) if(u2.gt.4.0 .and. degtemp.gt.0.3)then call brempde(u2,t2,gamma,qbrempde) qbrem = qbrempde if(qbrem.lt.0.0)then qbrem = 0.0 endif * * skip over other neutrino rates if we came here. * go to 50 endif if(u2.ge.1.0)then if(gamma.gt.1.0 .and. gamma.lt.171.0)then if(x(2).gt.0.00001)then call bremliqh(u2,t2,gamma,qbremlh) else qbremlh = 0.0 endif if(x(3).gt.0.00001)then call bremliqc(u2,t2,gamma,qbremlc) else qbremlc = 0.0 endif if(x(4).gt.0.00001)then call bremliqo(u2,t2,gamma,qbremlo) else qbremlo = 0.0 endif qbrem = qbremlh*x(2) + qbremlc*x(3) + qbremlo*x(4) elseif(gamma.ge.171.0)then if(x(3).gt.0.00001)then call bremsolc(u2,t2,gamma,qbremsc) else qbremsc = 0.0 endif if(x(4).gt.0.00001)then call bremsolo(u2,t2,gamma,qbremso) else qbremso = 0.0 endif qbrem = qbremsc*x(3) + qbremso*x(4) if(qbrem.lt.0.0)then qbrem = 0.0 endif endif endif * * convert to lumonisities * 50 eplas=qplas/rho en(1) = eplas ephot=qphot/rho en(2) = ephot epair=qpair/rho en(3) = epair erec=qrec/rho en(4) = erec * * output of Bremsstrahlung routines already is divided by rho * ebrem=qbrem en(5) = ebrem etot=eplas+ephot+epair+erec+ebrem fn=etot return end ************************************************************************