program approx1 c c Various drivers can be substituted for this routine. c The present version follows a path of rho prop T**3 c through carbon, neon, oxygen, and silicon burning. c It also seeks to maintain balanced power with c nuclear energy generation about three times the neutrino c energy loss rate. See notes for a justification. c The results are a little noisy. One could iterate on the c temperature each step and obtain smooher behavior but c for present purposes it did not seem worth the effort. c c If you change the initial conditions or any time step criteria c you will need to reset timer implicit real*8 (a-h,o-z) common /yycom/ yp,yhe3,yalp,yn14,yc12,yo16,yne,ymg,ysi,ys32, 1 ya36,yca,yti,ycr,yfe,yni,y54,yprot,yneut dimension x(19),y(19),dy(19),ai(19) equivalence (y(1),yp) data ai / 1.D0,3.D0,4.D0,14.D0,12.D0,16.D0,20.D0,24.D0,28.D0, 1 32.D0,36.D0,40.D0,44.D0,48.D0,52.D0,56.D0,54.D0,1.D0,1.D0 / call init jj = 1 do 10 i = 1,19 x(i) = 1.0d-20 y(i) = x(i)/ai(i) 10 continue c set temperature, density and initial composition c c 1 = h 7 = ne0 13 = ti44 19 = n-phot c 2 = he3 8 = mg24 14 = cr48 c 3 = he4 9 = si28 15 = fe52 c 4 = n14 10 = s32 16 = ni56 c 5 = c12 11 = ar36 17 = fe54 c 6 = o16 12 = ca40 18 = p-phot c t9 = 0.6 rho = 2.0d+5 x(5) = 0.20d0 x(6) = 0.80d0 y(5) = x(5)/ai(5) y(6) = x(6)/ai(6) c set initial time step and controllers h = 1.0d-10 iburn = 1 iold = 0 lburn = 0 irprox = 0 dtcp = 0.02d0 chimin = 1.d-4 jj = 1 dt9max = 0.01d0 write (6,20) (ai(n),n=5,12),ai(17) 20 format(' step',' time ',' dtime ',' t9 ',' rho', 1 ' enuc ','neutrino',9f10.0) snutmt1 = 1.0d0 snuwmt1 = 1.0d0 snucmt1 = 1.0d0 lburn = 0 iold = 0 irprox = 0 c this is the final time where the calculation stops c we want to run the clock backwards c you will need to run the program once to get c this value if you change the starting conditions time = 0.0d0 timer = 3.3457201432d+10 do 100 i = 1,2000 call burn(jj,t9,rho,h,dy,snuw,snubps,sneut,enuc,sdot,iburn, & snutmt1,snuwmt1,snucmt1,abar1,zbar1,etanx,iold,lburn,irprox) time = time+h timer = timer -h c try to follow path of enuc = 3 eneut c defac = (3.0d0*sneut-enuc)/abs(sneut+snuc) defac = (3.0d0*sneut-enuc)/abs(sneut+enuc) dt9 = 0.05d0*defac*t9 if (dt9.lt.(-dt9max)) dt9 = -dt9max if (dt9.gt.dt9max) dt9 = dt9max t9old=t9 t9 = t9+dt9 rho = rho*(t9/t9old)**3 c calculate new time step dp = 1.0d-10 do 31 n = 1,19 if (y(n).lt.chimin) go to 31 dp = max(dp,dy(n)/y(n)) 31 continue dtmult = min(2.0,dtcp/dp) h = h*dtmult c update abundances do 32 n = 1,19 y(n) = y(n)+dy(n) 32 x(n) = y(n)*ai(n) write (6,30) i,timer,h,t9,rho,enuc,sneut,(x(n),n=5,12),x(17) 30 format (i5,1pe10.2,14e10.2) if ((x(6)+x(9)).lt.1.0d-2) go to 101 100 continue c write final time 101 write (6,102) time 102 format ('final time = ', 1pe20.10) call exit(1) end subroutine init implicit real*8 (a-h,o-z) save common/ratc/rpp,r33,r34,rcpg,rnpg,ropg,r3a,rg3a,rnag,rcag, 1 r1212,r1216,r1616,roga,roag,rnega,rneag,rmgga,rmgag, 2 rsiga,rmgap,ralpa,ralpg,rsigp,rsiag,rsga,rsiap,rppa,rppg, 3 rsgp,rsag,rarga,rsap,rclpa,rclpg,rargp,rarag,rcaga,rarap, 4 rkpa,rkpg,rcagp,rcaag,rtiga,rcaap,rscpa,rscpg,rtigp, 5 rtiag,rcrga,rtiap,rvpa,rvpg,rcrgp,rcrag,rfega,rcrap, 6 rmnpa,rmnpg,rfegp,rfeag,rniga,rfeap,rcopa,rcopg,rnigp, 7 rfepg,rcogp,r52ng,r53gn,r53ng,r54gn,rheng,rhegn,rhng, 8 rdgn,rdpg,rhegp,rpen,rnep,xlam,sumb,fa,fg,r1,s1,t1,u1, 9 v1,w1,x1,c12agmlt,e1mltc12,e2mltc12,r3amult r3amult = 1.0d0 c12agmlt = 1.2d0 e1mltc12 = 1.7d0 e2mltc12 = 1.7d0 return end subroutine burn(jj,t9,rho,h,dy,snuw,snubps,sneut,enuc,sdot,iburn, & snutmt1,snuwmt1,snucmt1,abar1,zbar1,etanx,iold,lburn,irprox) implicit real*8 (a-h,o-z) save c.... for additional weak rates (N14(EC)) INTEGER nwwdim PARAMETER (nwwdim=1) REAL*8 recww(nwwdim),redww(nwwdim),rpdww(nwwdim), & eecww(nwwdim),eedww(nwwdim) INTEGER nww,nwwnot INTEGER idxww(nwwdim,2) COMMON /wwcom/ nww,nwwnot, & recww,redww,rpdww,eecww,eedww, & idxww INTEGER nwwc14 PARAMETER (nwwc14=1) c c this subroutine computes nucleosynthesis and energy generation c.... required input c jj = zone number (needed for bburn calls) c t9 = temperature in billions of degrees at new time point c rho =density in gm/cm**3 at new time point c h = time step c as input variable h tells program over what interval c to evolve abundances. c y = the vector of old abundances.this will be returned c unchanged. c iburn = burn flag. if if iburn.le.0, use 8 element c network ending with si28. otherwise use full c 17 element network. c snucmt1 = nuclear energy production multiplier. c snutmt1 = neutrino energy loss multiplier. c snuwmt1 = multiplier on weak neutrino losses c iold = 1 means circumvent alex's updates to c cno neutrino losses c = 0 use the updates cc lburn = 0 use this subroutine c = 1 call bburn instead to get values c from big BURN network. c irprox = 1 use rprox if in tn and yp range (default) c = 0 NEVER use rprox (useful for some Pop III stars) c = 2 use rprox network independent of yp mass fraction c = 3 ALWAYS use rprox network c.... output quantities c dy = the vector of abundance changes that occurs during c the given time step h. c enuc = nuclear energy generation in erg/gm/sec c sdot = net energy generation rate including neutrino c losses (erg/g/sec). c abar1 = the new average atomic weight. c zbar1 = the new average atomic number. c c sneut.....total neutrino energy loss rate (erg/g/sec). c snuc......nuclear energy production rate (erg/g/sec). c snuw......total "weak" neutrino energy loss rate (erg/g/sec). c snubps....total "pair" neutrino energy loss rate (erg/g/sec). c c.................................................................... c c y is defined as mass fraction divided by nucleon number c (moles/g). c c c yc12=c12 ymg=mg24 ya36=ar36 yti=ti44 yni=ni56 c yo16=o16 ysi=si28 yca=ca40 ycr=cr48 y54=fe54 c yne=ne20 ys32=s32 yalp=a yfe=fe52 ypd=p c yn14=n14 yp=pp yprot=photoprotons yneut=neutrons c yhe3=he3 c c pp are proton burning protons c c c c r1616=o16+o16 rcag=c12(ag) rarap=ar36(ap) c r1212=c12+c12 rclpa=cl35(pa) r3a=3a-c12 c rg3a=c12-3a rkpa=k39(pa) rsigp=si28(gp) c roga=o16(ga) roag=o16(ag) rsgp=s32(gp) c rnega=ne20(ga) rneag=ne20(ag) rargp=ar36(gp) c rmgga=mg24(ga) rmgag=mg24(ag) rcagp=ca40(gp) c rsiga=si28(ga) rmgap=mg24(ap) r1216=o16+c12 c rsga=s32(ga) rsiag=si28(ag) rppg=p31(pg) c rarga=ar36(ga) rsiap=si28(ap) ralpg=al27(pg) c rcaga=ca40(ga) rsag=s32(ag) rclpg=cl35(pg) c ralpa=al27(pa) rsap=s32(ap) rkpg=k39(pg) c rppa=p31(pa) rarag=ar36(ag) c c rcaag=ca40(ag) rtiga=ti44(ga) rscpa=sc43(pa) c rcaap=ca40(ap) rtigp=ti44(gp) rscpg=sc43(pg) c rtiag=ti44(ag) rcrga=cr48(ga) rvpa=v 47(pa) c rtiap=ti44(ap) rcrgp=cr48(gp) rvpg=v 47(pg) c rcrag=cr48(ag) rfega=fe52(ga) rmnpa=mn51(pa) c rcrap=cr48(ap) rfegp=fe52(gp) rmnpg=mn51(pg) c rfeag=fe52(ag) rniga=ni56(ga) rcopa=co55(pa) c rfeap=fe52(ap) rnigp=ni56(gp) rcopg=co55(pg) c rcogp=co55(gp) rfepg=fe54(pg) c rcpg=c12(pg) rnpg=n14(pg) ropg=o16(pg) c rnag=n14(ag) rpp=p+p c r34=he3+he4 r33=he3+he3 c rheng=he3(ng) rhegn=he4(gn) rhng=p(ng) c rdgn=2d(gn) rdpg=2d(pg) rhegp=3he(gp) c rpen=p(e,nu) rnep=n(e+,nu) c r52ng=fe52(ng) r53gn=fe53(gn) r53ng=fe53(ng) c r54gn=fe54(gn) c c xlam is mean electron capture rate on all z.gt.2 c c c.... number of screening factors (also set in screen) parameter (nscf=78) c.... common link to sleqs common/abuncom/aq(361),bq(19) c.... common links to rate, etc. common/ratc/rpp,r33,r34,rcpg,rnpg,ropg,r3a,rg3a,rnag,rcag, 1 r1212,r1216,r1616,roga,roag,rnega,rneag,rmgga,rmgag, 2 rsiga,rmgap,ralpa,ralpg,rsigp,rsiag,rsga,rsiap,rppa,rppg, 3 rsgp,rsag,rarga,rsap,rclpa,rclpg,rargp,rarag,rcaga,rarap, 4 rkpa,rkpg,rcagp,rcaag,rtiga,rcaap,rscpa,rscpg,rtigp, 5 rtiag,rcrga,rtiap,rvpa,rvpg,rcrgp,rcrag,rfega,rcrap, 6 rmnpa,rmnpg,rfegp,rfeag,rniga,rfeap,rcopa,rcopg,rnigp, 7 rfepg,rcogp,r52ng,r53gn,r53ng,r54gn,rheng,rhegn,rhng, 8 rdgn,rdpg,rhegp,rpen,rnep,xlam,sumb,fa,fg,r1,s1,t1,u1, 9 v1,w1,x1,c12agmlt,e1mltc12,e2mltc12,r3amult common /splrat/ ralf1,ralf2,r1f54,r2f54,r3f54,r4f54,r5f54, 1 r6f54,r7f54,r8f54 common /yycom/ yp,yhe3,yalp,yn14,yc12,yo16,yne,ymg,ysi,ys32, 1 ya36,yca,yti,ycr,yfe,yni,y54,yprot,yneut common/screenc/scfac(nscf) dimension be(19) dimension a(0:18,0:18),b(0:18) dimension y(19),dy(19) dimension ai(19) dimension zi(19) dimension rr(78) c.... arrays required for mazurek's ni56 electron capture fit dimension rnt(2),rne(2,7),datn(2,6,7),tv(7),rv(6) dimension rfdm(4),rfd0(4),rfd1(4),rfd2(4) dimension tfdm(5),tfd0(5),tfd1(5),tfd2(5) equivalence (aq,a) equivalence (bq,b) equivalence (y(1),yp) equivalence (rr(1),rpp) data rv / 6.D0,7.D0,8.D0,9.D0,10.D0,11.D0 / data tv / 2.D0,4.D0,6.D0,8.D0,10.D0,12.D0,14.D0 / data idweak / 0 / data be / 0.D0,7.71819D0,28.29603D0,104.65998D0,92.16294D0, 1 127.62093D0,160.64788D0,198.2579D0,236.5379D0,271.7825D0, 2 306.7202D0,342.0568D0,375.4772D0,411.469D0,447.708D0, 3 484.003D0,471.7696,0.D0,0.D0 / data ai / 1.D0,3.D0,4.D0,14.D0,12.D0,16.D0,20.D0,24.D0,28.D0, 1 32.D0,36.D0,40.D0,44.D0,48.D0,52.D0,56.D0,54.D0,1.D0,1.D0 / data zi / 1.D0,2.D0,2.D0,7.D0,6.D0,8.D0,10.D0,12.D0,14.D0,16.D0, 1 18.D0,20.D0,22.D0,24.D0,26.D0,28.D0,26.D0,1.D0,0.D0 / c.... electron capture data for ni56 (mazurek 1973) data (((datn(ii,ik,ij),ik=1,6),ij=1,7),ii=1,1) / + -3.98D0, -2.84D0, -1.41D0, 0.20D0, 1.89D0, 3.63D0, + -3.45D0, -2.62D0, -1.32D0, 0.22D0, 1.89D0, 3.63D0, + -2.68D0, -2.30D0, -1.19D0, 0.27D0, 1.91D0, 3.62D0, + -2.04D0, -1.87D0, -1.01D0, 0.34D0, 1.94D0, 3.62D0, + -1.50D0, -1.41D0, -0.80D0, 0.45D0, 1.99D0, 3.60D0, + -1.00D0, -0.95D0, -0.54D0, 0.60D0, 2.06D0, 3.58D0, + -0.52D0, -0.49D0, -0.21D0, 0.79D0, 2.15D0, 3.55D0 / data (((datn(ii,ik,ij),ik=1,6),ij=1,7),ii=2,2) / + -3.68D0, -2.45D0, -0.80D0, 1.12D0, 3.13D0, 5.19D0, + -2.91D0, -2.05D0, -0.64D0, 1.16D0, 3.14D0, 5.18D0, + -1.95D0, -1.57D0, -0.40D0, 1.24D0, 3.16D0, 5.18D0, + -1.16D0, -0.99D0, -0.11D0, 1.37D0, 3.20D0, 5.18D0, + -0.48D0, -0.40D0, 0.22D0, 1.54D0, 3.28D0, 5.16D0, + 0.14D0, 0.19D0, 0.61D0, 1.78D0, 3.38D0, 5.14D0, + 0.75D0, 0.78D0, 1.06D0, 2.07D0, 3.51D0, 5.11D0 / c.... don't take too small a timestep if(h.lt.1.d-99) h=1.d-99 c c.... calculate the constant cubic interpolation parameters c.... for mazurek's ni56 electron capture fit c.... if this is the first pass if(idweak.lt.1) then do k=2,4 rfdm(k)=1.D0/ & ((rv(k-1)-rv(k))*(rv(k-1)-rv(k+1))*(rv(k-1)-rv(k+2))) rfd0(k)=1.D0/ & ((rv(k)-rv(k-1))*(rv(k)-rv(k+1))*(rv(k)-rv(k+2))) rfd1(k)=1.D0/ & ((rv(k+1)-rv(k-1))*(rv(k+1)-rv(k))*(rv(k+1)-rv(k+2))) rfd2(k)=1.D0/ & ((rv(k+2)-rv(k-1))*(rv(k+2)-rv(k))*(rv(k+2)-rv(k+1))) enddo do j=2,5 tfdm(j)=1.D0/ & ((tv(j-1)-tv(j))*(tv(j-1)-tv(j+1))*(tv(j-1)-tv(j+2))) tfd0(j)=1.D0/ & ((tv(j)-tv(j-1))*(tv(j)-tv(j+1))*(tv(j)-tv(j+2))) tfd1(j)=1.D0/ & ((tv(j+1)-tv(j-1))*(tv(j+1)-tv(j))*(tv(j+1)-tv(j+2))) tfd2(j)=1.D0/ & ((tv(j+2)-tv(j-1))*(tv(j+2)-tv(j))*(tv(j+2)-tv(j+1))) enddo idweak=1 endif c.... if conditions are favorable for a rp-process, branch around c.... usual approx network and use rp-network in rprox. c.... range of validity of rprox is roughly .2 < t9 < 2.5. below c.... t9=.2 the cno-cycle may not be beta limited and the usual c.... approx network will do a better job (it also includes beta c.... limitation on the cno cycle). above t9=2.5 photodisintegration c.... reactions begin to impede the rp-process flow and, in any event, c.... alpha capture should dominate at such high temperature c.... and approx be appropriate. c.... the rp-process is an accelerated mechanism for c.... burning hydrogen and helium to heavy elements. only c.... employ the alternate network if hydrogen is present c.... note that yp is the abundance of hydrogen, not the abundance c.... of protons from photodisintegration (yprot) or the background c.... abundance of protons in advanced burning stages (not carried c.... in approx) ehvy=0.D0 epen=0.D0 enep=0.D0 snuni56=0.D0 c c usual approx network c zero matrix c 120 continue call qvset(0.d+0,aq,361) call qvset(0.d+0,bq,19) c c subroutine rate generates thermonuclear reaction rates c call rate(t9,rho,iburn) c.... calculate averages for screening routine abar=0.D0 zbar=0.D0 ytot1=0.D0 z2bar=0.D0 do iy=1,19 ytot1=ytot1+y(iy) z2bar=z2bar+zi(iy)**2*y(iy) abar=abar+ai(iy)*y(iy) zbar=zbar+zi(iy)*y(iy) enddo z2bar=z2bar/ytot1 abar=abar/ytot1 zbar=zbar/ytot1 c.... get screening factors and modify reaction rates c.... c.... screening routine revised as per r.k. wallace july 1982 c.... call screen(t9*1.d+9,rho,abar,zbar,z2bar,ytot1) do i1=1,78 rr(i1)=rr(i1)*scfac(i1) enddo c.... c.... check for beta limitation of he3+he4, 14n(pg), and 16o(pg). c.... assume steady state abundances of 14o and 15o in beta limited c.... cno cycle (ie. divide 15o decay rate by 1.57). c.... if (yalp.gt.0.D0) r34=min(r34,0.896d+0/yalp) if (yp.gt.0.D0) rnpg=min(rnpg,5.68d-03/(yp*1.57d+0)) if (yp.gt.0.D0) ropg=min(ropg,0.0105d+0/yp) c c Nov. 27, 1992 - put in very approximate coding for high rho correction c to nitrogen burning c c Above rho = 5. x 10**5 use a rate for 14n(ag) that can be c accelerated by electron capture to 14c and subsequent alpha c capture (actually to 18o, but assumed 20ne here). See hashimoto c et al apj 307 687 (1986) for discussion. see esp fig 1. use only c the third thru fifth terms of their eq. 1 which should for 0.03 le c t9. do not allow alpha capture to proceed faster than electron c capture though, which in turn must be faster than 14n(ag) snun14=0.0D0 c.... c.... zero electron capture rates and neutrino losses c.... note call to ecap has been removed (see keplerh). c.... revised nov.25, 1981 c.... a call to ecapnuc is issued to get rpen and spen which c.... are used to calculate high density hydrogen burning (see c.... matrix elements for hydrogen burning) c.... xlam=0.D0 rpen=0.D0 rnep=0.D0 ehvy=0.D0 epen=0.D0 enep=0.D0 ye=0.5D0*(1.D0+yp+yhe3-2.d0*y54+yprot-yneut) c.... calculate ni56 electron capture and neutrino loss rates rn56ec=0.D0 snuni56=0.D0 if((t9.lt.2.0D0).or.(rho*ye.lt.1.d+6)) go to 20 tm=t9*1.d+9 if(tm.gt.1.4d+10) tm=1.4d+10 if(tm.lt.2.d+9) tm=2.d+9 t9m=tm*1.d-9 r=log10(rho*ye) r=min(11.d+0,r) r=max(6.d+0,r) jp=min(max(2,int(0.5d+0*t9m)),5) kp=min(max(2,int(r)-5),4) rfm=r-rv(kp-1) rf0=r-rv(kp) rf1=r-rv(kp+1) rf2=r-rv(kp+2) dfacm=rf0*rf1*rf2*rfdm(kp) dfac0=rfm*rf1*rf2*rfd0(kp) dfac1=rfm*rf0*rf2*rfd1(kp) dfac2=rfm*rf0*rf1*rfd2(kp) tfm=t9m-tv(jp-1) tf0=t9m-tv(jp) tf1=t9m-tv(jp+1) tf2=t9m-tv(jp+2) tfacm=tf0*tf1*tf2*tfdm(jp) tfac0=tfm*tf1*tf2*tfd0(jp) tfac1=tfm*tf0*tf2*tfd1(jp) tfac2=tfm*tf0*tf1*tfd2(jp) do jr=1,2 do jd=jp-1,jp+2 rne(jr,jd)=dfacm*datn(jr,kp-1,jd)+dfac0*datn(jr,kp,jd) 1 +dfac1*datn(jr,kp+1,jd)+dfac2*datn(jr,kp+2,jd) enddo rnt(jr)=tfacm*rne(jr,jp-1)+tfac0*rne(jr,jp) 1 +tfac1*rne(jr,jp+1)+tfac2*rne(jr,jp+2) enddo rn56ec=(10.d+0)**rnt(1) snuni56=6.022548d+23*8.18683d-7*yni*(10.d+0)**rnt(2) 20 continue c.... calculate mass sums and mass deficit c.... note call to ratfix has been removed (see keplerh) ytot=0.D0 do i1=1,19 ytot=ytot+ai(i1)*y(i1) enddo xmas=abs(1.d+0-ytot) c.... calculate branching ratios for dummy proton links c.... and special combined rates for alpha photodisintegration and fe54 r1=0.D0 s1=0.D0 t1=0.D0 u1=0.D0 v1=0.D0 w1=0.D0 x1=0.D0 z1=0.D0 cr=0.D0 den=0.D0 ralf1=0.D0 ralf2=0.D0 den1=0.D0 den2=0.D0 r1f54=0.D0 r2f54=0.D0 r3f54=0.D0 r4f54=0.D0 r5f54=0.D0 r6f54=0.D0 r7f54=0.D0 r8f54=0.D0 yneut2=yneut*yneut yprot2=yprot*yprot if (iburn.le.0) go to 60 r1=ralpa/(ralpa+ralpg) s1=rppa/(rppa+rppg) t1=rclpa/(rclpa+rclpg) u1=rkpa/(rkpa+rkpg) v1=rscpa/(rscpa+rscpg) w1=rvpa/(rvpa+rvpg) x1=rmnpa/(rmnpa+rmnpg) den=rhegp*rdgn+yneut*rheng*rdgn+yneut*yprot*rheng*rdpg if(den.lt.1.d-150) go to 70 ralf1=rhegn*rhegp*rdgn/den ralf2=rheng*rdpg*rhng/den c 70 continue den1=(rcogp+yprot*(rcopg+rcopa)) den2=(r53gn+yneut*r53ng) if(den1.gt.1.d-99) den1=1.D0/den1 if(den2.gt.1.d-99) den2=1.D0/den2 c.... don't consider 54fe photodisintegration if temperature is c.... too low or there is still free hydrogen around. c.... (would like to use parameters here someday) if (t9.lt.1.5D0.or.yp.gt.1.d-05) then den1=0.D0 den2=0.D0 endif r1f54=r54gn*r53gn*den2 r2f54=r52ng*r53ng*den2 r3f54=rfepg*rcopg*den1 r4f54=rnigp*rcogp*den1 r5f54=rfepg*rcopa*den1 r6f54=rfeap*rcogp*den1 r7f54=rfeap*rcopg*den1 r8f54=rnigp*rcopa*den1 60 continue c.... calculate the terms of the coupled, linearized rate equations c.... for the abundances (a*dy=b) c c 1) hydrogen c a(1,1)=1.D0/h+6.D0*yp*rpp+2.D0*yc12*rcpg+ & 2.D0*yn14*rnpg+2.D0*yo16*ropg a(1,2)=-4.D0*yhe3*r33+yalp*r34 a(1,3)=yhe3*r34 a(1,4)=2.D0*yp*rnpg a(1,5)=2.D0*yp*rcpg a(1,6)=2.D0*yp*ropg b(1)=-3.D0*yp*yp*rpp+2.D0*yhe3*yhe3*r33-yhe3*yalp*r34 1 -2.D0*yc12*yp*rcpg-2.D0*yn14*yp*rnpg-2.D0*yo16*yp*ropg c c add special matrix elements for high density hydrogen burning c assume 3 protons make he3 at a rate given by that of electron c capture on free protons. call to ecapnuc gives weak rate. c c if etanx is not properly defined (from call to es during c previous cycle) skip such captures on this pass. c a(1,1)=a(1,1)+3.D0*rpen b(1)=b(1)-3.D0*yp*rpen c c 2) he3 c a(2,1)=-2.D0*yp*rpp a(2,2)=1.D0/h+4.D0*yhe3*r33+yalp*r34 a(2,3)=yhe3*r34 b(2)=yp*yp*rpp-2.D0*yhe3*yhe3*r33-yhe3*yalp*r34 c c add high density h-burn matrix elements c a(2,1)=a(2,1)-rpen b(2)=b(2)+yp*rpen c c for now don't compute nucleon capture rate on protons c or neutrons from photodisintegration. c rpen1=rpen rnep1=rnep rpen=0.D0 rnep=0.D0 c c 3) helium (or alpha) c a(3,1)=-yn14*fa*rnpg-yo16*ropg a(3,2)=-2.D0*yhe3*r33-yalp*r34 ca a33 = 1.D0/h-r1*ymg*rmgap-s1*ysi*rsiap-t1*ys32*rsap-u1*ya36* ca 1 rarap+ yne*rneag+ymg*(rmgag+rmgap)+ysi*(rsiag ca 2 +rsiap)+ys32*(rsag+rsap)+ya36*(rarag+rarap)+yc12*rcag ca 3 +yo16*roag ca a33=a33+ ca & 9.0D0*yalp*yalp*r3a-v1*yca*rcaap-w1*yti*rtiap-x1*ycr*rcrap ca 1 +yca*(rcaag+rcaap)+yti*(rtiag+rtiap)+ycr* ca 2 (rcrag+rcrap)+yfe*(rfeag+r6f54+yprot*r7f54)+ralf1 a33 = 1.D0/h+yne*rneag+yc12*rcag+yo16*roag+9.D0*yalp*yalp*r3a a(3,3)=a33-yhe3*r34+yn14*rnag*1.5D0 a(3,4)=-yp*fa*rnpg+yalp*rnag*1.5D0 a(3,5)=-3.D0*rg3a-2.D0*yc12*r1212+yalp*rcag-0.5D0*yo16*r1216 a(3,6)=-roga+yalp*roag-0.5D0*yc12*r1216-yp*ropg a(3,7)=-rnega+yalp*rneag a(3,8)=-rmgga c.... test for high-z bypass if (iburn.le.0) go to 260 a(3,3)=a(3,3) & -r1*ymg*rmgap-s1*ysi*rsiap-t1*ys32*rsap-u1*ya36* & rarap+ymg*(rmgag+rmgap)+ysi*(rsiag & +rsiap)+ys32*(rsag+rsap)+ya36*(rarag+rarap) & -v1*yca*rcaap-w1*yti*rtiap-x1*ycr*rcrap & +yca*(rcaag+rcaap)+yti*(rtiag+rtiap)+ycr* & (rcrag+rcrap)+yfe*(rfeag+r6f54+yprot*r7f54)+ralf1 a(3,6)=a(3,6)-yo16*r1616*2.D0*(0.56D0+0.34D0*s1) a(3,8)=a(3,8)+yalp*(rmgag+rmgap*(1.D0-r1)) a(3,9)=-rsiga-r1*rsigp-s1*yalp*rsiap+yalp*(rsiag+rsiap) a(3,10)=-rsga-s1*rsgp-t1*yalp*rsap+yalp*(rsag+rsap) a(3,11)=-rarga-t1*rargp-u1*yalp*rarap+yalp*(rarag+rarap) a(3,12)=-rcaga-u1*rcagp-v1*yalp*rcaap+yalp*(rcaag+rcaap) a(3,13)=-rtiga-v1*rtigp-w1*yalp*rtiap+yalp*(rtiag+rtiap) a(3,14)=-rcrga-w1*rcrgp-x1*yalp*rcrap+yalp*(rcrag+rcrap) a(3,15)=-rfega-x1*rfegp+yalp*(rfeag+r6f54+yprot*r7f54) a(3,16)=-rniga-yprot*r8f54 a(3,17)=-yprot2*r5f54 a(3,18)=-2.D0*yneut2*yprot*ralf2-yni*r8f54+yfe*yalp*r7f54 1 -2.D0*y54*yprot*r5f54 a(3,0)=-2.D0*yprot2*yneut*ralf2 b1=yhe3*yhe3*r33+yhe3*yalp*r34+yn14*yp*fa*rnpg+yo16*yp*ropg 1 -yalp*yn14*rnag*1.5D0 b2 = -yc12*yc12*r1212-0.56D0*yo16*yo16*r1616-r1*yalp*ymg*rmgap 1 -s1*yalp*ysi*rsiap-t1*yalp*ys32*rsap-u1*yalp*ya36*rarap 2 +yalp*yo16*roag+yalp*yne*rneag+yalp*yc12*rcag+yalp*ymg 3 *(rmgag+rmgap)+yalp*ysi*(rsiag+rsiap)+yalp*ys32*(rsag+rsap) 4 +yalp*ya36*(rarag+rarap)-0.5D0*yc12*yo16*r1216 b3=-v1*yalp*yca*rcaap-w1*yalp*yti*rtiap-x1*yalp*ycr*rcrap 1+yfe*yalp*yprot*r7f54+yalp*yca*(rcaag+rcaap)+yalp*yti*(rtiag+ 2 rtiap)+yalp*ycr*(rcrag+rcrap)+yalp*yfe*(rfeag+r6f54) 3 +3.D0*yalp**3*r3a-s1*yo16*yo16*0.34D0*r1616 b4=3.D0*yc12*rg3a+yo16*roga+yne*rnega+ymg*rmgga 1 +ysi*(rsiga+r1*rsigp) + ys32*(rsga+s1*rsgp) 2 +ya36*(rarga+t1*rargp) + yca*(rcaga+u1*rcagp) 3 +yti*(rtiga+v1*rtigp) + ycr*(rcrga+w1*rcrgp) 4 +yfe*(rfega+x1*rfegp) + yni*(rniga+yprot*r8f54) b5=-yalp*ralf1+yneut2*yprot2*ralf2+y54*yprot2*r5f54 go to 270 c.... high-z element contributions omitted 260 continue b1=yhe3**2*r33+yhe3*yalp*r34+yn14*yp*fa*rnpg+yo16*yp*ropg 1 -yalp*yn14*rnag*1.5D0 ca b2 = -yc12*yc12*r1212-0.56D0*yo16*yo16*r1616 ca 2 +yalp*yo16*roag+yalp*yne*rneag+yalp*yc12*rcag ca 4 -0.5D0*yc12*yo16*r1216 ca b3 =+3.D0*yalp**3*r3a-s1*yo16*yo16*0.34D0*r1616 b2 =-yc12**2*r1212 2 +yalp*yo16*roag+yalp*yne*rneag+yalp*yc12*rcag 4 -0.5D0*yc12*yo16*r1216 b3 =+3.D0*yalp**3*r3a b4=3.D0*yc12*rg3a+yo16*roga+yne*rnega+ymg*rmgga b5=0.D0 270 continue b(3)=b1-b2-b3+b4+b5 c c 4) nitrogen 14 c a(4,1)=-yc12*rcpg+yn14*rnpg-yo16*ropg a(4,3)=yn14*rnag a(4,4)=1.D0/h+yp*rnpg+yalp*rnag a(4,5)=-yp*rcpg a(4,6)=-yp*ropg b(4)=yc12*yp*rcpg-yn14*yp*rnpg+yo16*yp*ropg-yalp*yn14*rnag c c 5) carbon 12 c a(5,1)=yc12*rcpg-yn14*fa*rnpg a(5,3)=yc12*rcag-3.D0*yalp**2*r3a a(5,4)=-yp*fa*rnpg ca a(5,5)=1.D0/h+4.D0*yc12*r1212+rg3a+yalp*rcag+yo16*r1216+yp*rcpg a(5,5)=1.D0/h+4.D0*yc12*r1212+rg3a+yalp*rcag+yp*rcpg ca a(5,6)=yc12*r1216-roga a(5,6)=-roga ca b(5)=-2.D0*yc12*yc12*r1212-yc12*yalp*rcag-yc12*yo16*r1216 ca 1 +yalp**3*r3a-yc12*rg3a+yo16*roga-yc12*yp*rcpg ca 2 +yn14*yp*fa*rnpg b(5)=-2.D0*yc12*yc12*r1212-yc12*yalp*rcag 1 +yalp**3*r3a-yc12*rg3a+yo16*roga-yc12*yp*rcpg 2 +yn14*yp*fa*rnpg IF (iburn.GT.0) THEN a(5,5)=a(5,5)+yo16*r1216 a(5,6)=a(5,6)+yc12*r1216 b(5)=b(5)-yc12*yo16*r1216 ELSE a(5,5)=a(5,5)+0.5D0*yo16*r1216 a(5,6)=a(5,6)+0.5D0*yc12*r1216 b(5)=b(5)-0.5D0*yc12*yo16*r1216 ENDIF c c 6) oxygen 16 c a(6,1)=-yn14*fg*rnpg+yo16*ropg a(6,3)=yo16*roag-yc12*rcag a(6,4)=-yp*fg*rnpg ca a(6,5)=yo16*r1216-yalp*rcag a(6,5)=-yalp*rcag ca a(6,6)=1.D0/h+4.D0*yo16*r1616+roga+yalp*roag+yc12*r1216+yp*ropg a(6,6)=1.D0/h+roga+yalp*roag+yp*ropg a(6,7)=-rnega b(6)=-yo16*yalp*roag 1 +yc12*yalp*rcag-yo16*roga+yne*rnega+yn14*yp*fg*rnpg 2 -yo16*yp*ropg IF (iburn .GT. 0) THEN a(6,5)=a(6,5)+yo16*r1216 a(6,6)=a(6,6)+yc12*r1216+4.D0*yo16*r1616 b(6)=b(6)-2.D0*yo16**2*r1616-yc12*yo16*r1216 ELSE a(6,5)=a(6,5)+0.5D0*yo16*r1216 a(6,6)=a(6,6)+0.5D0*yc12*r1216 b(6)=b(6)-0.5D0*yc12*yo16*r1216 ENDIF c c 7) neon 20 c a(7,4)=- yalp*rnag a(7,5)=-2.D0*yc12*r1212 a(7,6)=-yalp*roag a(7,7)=1.D0/h+rnega+yalp*rneag a(7,8)=-rmgga a(7,3)=yne*rneag-yo16*roag-yn14*rnag b(7)=yc12*yc12*r1212-yne*yalp*rneag+yo16*yalp*roag 1 -yne*rnega+ymg*rmgga+ yn14*yalp*rnag c c 8) magnesium 24 c ca a(8,3)=ymg*(rmgag+rmgap*(1.D0-r1))-yne*rneag ca a(8,8)=1.D0/h+rmgga+yalp*(rmgag+rmgap*(1.D0-r1)) ca if(iburn.gt.0) a(8,9)=-rsiga-r1*rsigp ca b(8)=-ymg*yalp*(rmgag+rmgap*(1.D0-r1))+0.5D0*yc12*yo16*r1216 ca 1 +yne*yalp*rneag-ymg*rmgga+ysi*(rsiga+r1*rsigp) a(8,3)=-yne*rneag a(8,5)=-yo16*r1216*0.5D0 a(8,6)=-yc12*r1216*0.5D0 a(8,7)=-yalp*rneag a(8,8)=1.D0/h+rmgga b(8)=0.5D0*yc12*yo16*r1216+yne*yalp*rneag-ymg*rmgga IF (iburn.gt.0) THEN a(8,3)=a(8,3)+ymg*(rmgag+rmgap*(1.D0-r1)) a(8,8)=a(8,8)+yalp*(rmgag+rmgap*(1.D0-r1)) a(8,9)=-rsiga-r1*rsigp b(8)=b(8)-ymg*yalp*(rmgag+rmgap*(1.D0-r1))+ & ysi*(rsiga+r1*rsigp) ENDIF c.... test for high-z bypass if (iburn.le.0) go to 130 c c 9) silicon 28 c a(9,5)=-yo16*r1216*0.5D0 a(9,6)=-1.12D0*yo16*r1616-0.68D0*s1*yo16*r1616-0.5D0*yc12*r1216 a(9,8)=-yalp*(rmgag+(1.D0-r1)*rmgap) a(9,9)=1.D0/h+r1*rsigp+rsiga+yalp*(rsiag+rsiap)-s1*yalp*rsiap a(9,10)=-s1*rsgp-rsga a(9,3)=-ymg*(rmgag+(1.D0-r1)*rmgap)+ysi*(rsiag+rsiap*(1.D0-s1)) b(9)=ymg*yalp*(rmgag+rmgap*(1.D0-r1))+0.56D0*yo16**2*r1616 1 -ysi*yalp*(rsiag+rsiap*(1.D0-s1))+0.34D0*yo16**2*s1*r1616 2 +0.5D0*yc12*yo16*r1216-ysi*(r1*rsigp+rsiga)+ys32*(rsga+s1 3 *rsgp) c c 10) sulfur 32 c a(10,3)=-ysi*rsiag+(s1-1.D0)*ysi*rsiap+ys32*(rsag+rsap) 1 -t1*ys32*rsap a(10,6)=-0.2D0*yo16*r1616+(s1-1.D0)*0.68D0*yo16*r1616 a(10,9) =-yalp*rsiag+(s1-1.D0)*yalp*rsiap a(10,10)=1.D0/h+rsga+s1*rsgp+yalp*(rsag+rsap)-t1*yalp*rsap a(10,11)=-t1*rargp-rarga b(10)=ysi*yalp*(rsiag+rsiap*(1.D0-s1))+0.34D0*yo16*yo16*r1616 1 *(1.D0-s1)-ys32*yalp*(rsag+rsap*(1.D0-t1))+0.1D0*yo16* 2 yo16*r1616-ys32*(rsga+s1*rsgp)+ya36*(rarga+t1*rargp) c c 11) argon 36 c a(11,3)=(t1-1.D0)*ys32*rsap-ys32*rsag+ya36*(rarag+rarap) 1 -u1*ya36*rarap a(11,10)=(t1-1.D0)*yalp*rsap-yalp*rsag a(11,11)=1.D0/h+rarga+t1*rargp+yalp*(rarag+rarap)-u1*yalp*rarap a(11,12)=-u1*rcagp-rcaga b(11)=ys32*yalp*(rsag+rsap*(1.D0-t1))-ya36*yalp*(rarag+rarap* 1 (1.D0-u1))-ya36*(rarga+t1*rargp)+yca*(rcaga+rcagp*u1) c c 12) calcium 40 c a(12,3)=(u1-1.D0)*ya36*rarap-ya36*rarag+yca*(rcaag+rcaap) 1 -v1*yca*rcaap a(12,11)=(u1-1.D0)*yalp*rarap-yalp*rarag a(12,12)=1.D0/h+rcaga+u1*rcagp+yalp*(rcaag+rcaap)-v1*yalp*rcaap a(12,13)=-v1*rtigp-rtiga b(12)=ya36*yalp*(rarag+rarap*(1.D0-u1))-yca*yalp*(rcaag+rcaap* 1 (1.D0-v1))-yca*(rcaga+rcagp*u1)+yti*(rtiga+rtigp*v1) c c 13) titanium 44 c a(13,3)=(v1-1.D0)*yca*rcaap-yca*rcaag+yti*(rtiag+rtiap) 1 -w1*yti*rtiap a(13,12)=(v1-1.D0)*yalp*rcaap-yalp*rcaag a(13,13)=1.D0/h+rtiga+v1*rtigp+yalp*(rtiag+rtiap)-w1*yalp*rtiap a(13,14)=-w1*rcrgp-rcrga b(13)=yca*yalp*(rcaag+rcaap*(1.D0-v1))-yti*yalp*(rtiag+rtiap* 1 (1.D0-w1))-yti*(rtiga+v1*rtigp)+ycr*(rcrga+w1*rcrgp) c c 14) chromium 48 c a(14,3)=(w1-1.D0)*yti*rtiap-yti*rtiag+ycr*(rcrag+rcrap) 1 -x1*ycr*rcrap a(14,13)=(w1-1.D0)*yalp*rtiap-yalp*rtiag a(14,14)=1.D0/h+rcrga+w1*rcrgp+yalp*(rcrag+rcrap)-x1*yalp*rcrap a(14,15)=-x1*rfegp-rfega b(14)=yti*yalp*(rtiag+rtiap*(1.D0-w1))-ycr*yalp*(rcrag+rcrap* 1 (1.D0-x1))-ycr*(rcrga+w1*rcrgp)+yfe*(rfega+x1*rfegp) c c 15) iron 52 c a(15,3)=(x1-1.D0)*ycr*rcrap-ycr*rcrag+yfe*(rfeag+r6f54 1 +yprot*r7f54) a(15,14)=(x1-1.D0)*yalp*rcrap-yalp*rcrag a(15,15)=yalp*(rfeag+r6f54+yprot*r7f54)+yneut2*r2f54 1 +rfega+x1*rfegp+1.D0/h a(15,16)=-rniga-yprot*r8f54 a(15,17)=-yprot2*r5f54-r1f54 a(15,18)=-2.D0*y54*yprot*r5f54+yfe*yalp*r7f54-yni*r8f54 a(15,0)=2.D0*yfe*yneut*r2f54 b(15)=ycr*yalp*(rcrag+(1.D0-x1)*rcrap)-yfe*(rfega+x1*rfegp) 1 -yfe*yalp*(rfeag+r6f54)+yni*rniga+y54*yprot2*r5f54 2 -yfe*yalp*yprot*r7f54+yni*yprot*r8f54 3 -yfe*yneut2*r2f54+y54*r1f54 c c 16) nickel 56 c a(16,3)=-yfe*(yprot*r7f54+rfeag) a(16,15)=-yalp*(rfeag+yprot*r7f54) a(16,16)=rniga+r4f54+yprot*r8f54+1.D0/h+rn56ec a(16,17)=-yprot2*r3f54 a(16,18)=-yfe*yalp*r7f54+yni*r8f54-2.D0*y54*r3f54*yprot b(16)=yfe*yalp*(rfeag+yprot*r7f54)-yni*(rniga+r4f54 1 +yprot*r8f54)+y54*yprot2*r3f54-rn56ec*yni c c 17) iron 54 coupled by a double neutron link to fe52 c and a double proton link to ni56 c a(17,3)=-yfe*r6f54 a(17,15)=-yneut2*r2f54-yalp*r6f54 a(17,16)=-r4f54-56.D0*rn56ec/54.D0 a(17,17)=r1f54+yprot2*(r3f54+r5f54)+1.D0/h a(17,18)=2.D0*y54*yprot*(r3f54+r5f54) a(17,0)=-2.D0*yfe*yneut*r2f54 b(17)=-y54*r1f54+yfe*yneut2*r2f54-y54*yprot2*(r3f54+r5f54) 1 +yni*r4f54+yfe*yalp*r6f54+56.D0*rn56ec*yni/54.D0 c c 18) protons from alpha photodisintegration,fe54, c and weak interactions c a(18,3)=-2.D0*yfe*r6f54-2.D0*ralf1 a(18,15)=-2.D0*yalp*r6f54 a(18,16)=-2.D0*r4f54 a(18,17)=2.D0*yprot2*(r3f54+r5f54) a(18,18)=4.D0*yprot*yneut2*ralf2+rpen+4.D0*y54*yprot* 1 (r3f54+r5f54)+1.D0/h a(18,0)=-rnep+4.D0*yprot2*yneut*ralf2 b(18)=-2.D0*yprot2*yneut2*ralf2+2.D0*yalp*ralf1-yprot*rpen 1 +yneut*rnep-2.D0*y54*yprot2*(r3f54+r5f54) 2 +2.D0*yni*r4f54+2.D0*yfe*yalp*r6f54 c c 0/19) neutrons from alpha photodisintegration,54fe, c and weak interactions c a(0,0)=4.D0*yfe*yneut*r2f54+4.D0*yprot2*yneut*ralf2+rnep+1.D0/h a(0,3)=-2.D0*ralf1 a(0,15)=2.D0*yneut2*r2f54 a(0,17)=-2.D0*r1f54 a(0,18)=4.D0*yprot*yneut2*ralf2-rpen b(0)=-2.D0*yfe*yneut2*r2f54-2.D0*yprot2*yneut2*ralf2-yneut*rnep 1 +2.D0*yalp*ralf1+2.D0*y54*r1f54+yprot*rpen c c now invert the matrix to obtain the differential vector c of abundances b. c.... rotate matrix for more efficient inversion in the 19x19 case call qvrota(aq,361) call qvrota(bq,19) call sleqs(19,19,4,4,16,19) call qvrota(bq,19) go to 250 c.... invert repacked, unrotated matrix in 8x8 case 130 continue lf=0 do i1=1,8 lf1=lf+1 lf21=lf+21 do iq=0,7 aq(lf1+iq)=aq(lf21+iq) enddo lf=lf+19 bq(i1)=bq(i1+1) enddo call sleqs(8,19,0,0,0,0) do i1=8,1,-1 bq(i1+1)=bq(i1) enddo bq(1)=0.D0 c c.... compute energy generation using renormalized dy's to c.... remove potential glitches due to mass nonconservation c.... also calculate the new average atomic weitght and number c 250 continue do i1=1,18 dy(i1)=b(i1) enddo dy(19)=b(0) edot=0.D0 dytot=0.D0 abar1=0.D0 enbar1=0.D0 zbar1=0.D0 do i1=1,19 dytot=dytot+ai(i1)*dy(i1) yi1=max(0.d+0,y(i1)+dy(i1)) abar1=ai(i1)*yi1+abar1 zbar1=zi(i1)*yi1+zbar1 enbar1=yi1+enbar1 enddo abar1=abar1/enbar1 zbar1=zbar1/enbar1 yfac=dytot/ytot do i1=1,19 edot=edot+(dy(i1)-yfac*y(i1))*be(i1) enddo enuc=edot*1.602d-06*6.02259d+23/h c c call sneut for neutrino losses (only if t9.gt 0.1) c REMOVED by alex (UTC 20020320025140) c 110 continue sneut=0.D0 snubps=0.0D0 shvy=0.D0 spen=0.D0 snep=0.D0 snuh=0.D0 snuw=0.D0 c REMOVED by alex (UTC 20020320025141) c if (t9.lt.0.10D0) go to 27 temp=t9*1.D9 iolde=iold*2 call sneutr(temp,rho,zbar1,abar1,snubps,iolde) c.... do not use the if rprox was called c.... REASON: rates were not computed. c.... ADDED by alex (UTC 20060220024415) if (iirprox .eq. 1) goto 140 do ji=4,17 shvy=shvy+ehvy*(y(ji)+dy(ji)) enddo spen=epen*(y(18)+dy(18)) snep=enep*(y(19)+dy(19)) c REMOVED by alex (UTC 20020320025141) c 27 continue IF (iold .EQ. 0) THEN c.... add neutrino losses for H burning (A. Heger 30 Jul 1997) c.... pp1 snuh=yp*yp*rpp*0.263D0 c.... pp2 reaction be7pg = 3.11d+05*t9**(-2.d0/3.d0)* & exp(-10.262D0*t9**(-1.d0/3.d0)) & +2.53d+03*t9**(-1.5d0) & *exp(-7.306D0/t9)*rho*yp c.... pp3 reaction be7ec=1.50558D-7 if ((t9 .le. 3.0D0).AND.(T9.GT.0.001)) then be7ec = 1.34d-10*t9**(-0.5d0)* & (1.D0-0.537D0*t9**(1.d0/3.d0) & +3.86D0*t9**(2.d0/3.d0)+ & 0.0027D0/t9 & *exp(2.515d-03/t9)) endif c.... pp2 and pp3 added in with branching snuh=snuh+(0.8d0*be7ec+7.2*be7pg)/(be7ec+be7pg) & *yhe3*yalp*r34 c$$$ if (t9 .lt. 0.023D0) then c$$$c.... pp2 c$$$ snuh=snuh+yhe3*yalp*r34*0.8D0 c$$$ else c$$$c.... pp3 c$$$ snuh=snuh+yhe3*yalp*r34*7.2D0 c$$$ endif c.... CNO cycle c.... N13-->C13+e+nu snuh=snuh+yc12*yp*rcpg*0.71D0 c.... O15-->N15+e+nu snuh=snuh+yn14*yp*rnpg*1.0D0 c.... F17-->O17+e+nu snuh=snuh+yo16*yp*ropg*0.94D0 c.... get the units right snuh=snuh*1.602E-6*6.02259d+23 c.... add up weak neurino losses snuw=shvy+spen+snep+snuni56+snun14+snuh c.... account for rest mass difference n-(p+e) c.... (important in weak reactions only) c.... (A. Heger 31 Jul 1997) enuc=enuc-0.782D0*1.602D-6*6.02259d+23*( c.... pp1 & yp*yp*rpp+ c.... pp2&3 & yhe3*yalp*r34+ c.... N13-->C13+e+nu & yc12*yp*rcpg+ c.... O15-->N15+e+nu & yn14*yp*rnpg+ c.... F17-->O17+e+nu & yo16*yp*ropg) ENDIF 140 continue c.... calculate total energy generation and loss rates snuw=snuw*snuwmt1*snutmt1 snubps=snubps*snutmt1 enuc=enuc*snucmt1 sneut=snubps+snuw sdot=enuc-sneut end subroutine leqs(n) implicit real*8 (a-h,o-z) save c.... c.....leqs performs matrix inversion c.... c n-dimension of a common /acl/ amat(10,10),bvec(10) n1=n-1 c find maximum element in each row, and divide do 19 i=1,n r=abs(amat(i,1)) do 16 j=2,n c=abs(amat(i,j)) if (r.lt.c) r=c 16 continue do 17 j=1,n 17 amat(i,j)=amat(i,j)/r bvec(i)=bvec(i)/r 19 continue do 9 j=1,n1 l=j+1 c.....no pivoting do 9 i=l,n r=-amat(i,j) if (r.eq.0.0) go to 9 r=r/amat(j,j) do 7 k=l,n 7 amat(i,k)=amat(i,k)+r*amat(j,k) bvec(i)=bvec(i)+r*bvec(j) 9 continue c.... c the matrix is now in triangular form c start the back substitution and find solution c.... bvec(n)=bvec(n)/amat(n,n) do 13 l=1,n1 i=n-l r=0.0 imax=i+1 do 12 j=imax,n jj=i+n+1-j 12 r=r+amat(i,jj)*bvec(jj) 13 bvec(i)=(bvec(i)-r)/amat(i,i) return end subroutine sleqs(nn,nnp,nbloc,ntrid,lbord1,lbord2) implicit real*8 (a-h,o-z) save c optimize c.... this subroutine solves the sparse linear system a*dy=b, c.... returning the solution for dy in b. nn is the dimension of c.... the matrix to be inverted, nnp is the total dimension c.... of the matrix in which the matrix to be inverted is c.... embedded. an nbloc x nbloc matrix is assumed to reside in c.... the upper left corner of a,accompanied by two full columns c.... at positions lbord1 and lbord2. after this a tridiagonal c.... pattern is assumed to persist for ntrid rows accompanied c.... by a filled column at lbord1. subsequent rows can be general. common/abuncom/a(361),b(19) n=nn np=nnp n1=n-1 n2=n-2 lnb1=(lbord1-1)*np lnb2=(lbord2-1)*np ntridx=nbloc+ntrid c.... no conditioning c.... reduce matrix to upper triangular form taking advantage c.... of matrix sparseness lj=0 do 80 j=1,n1 l=j+1 lf0=j*np ajji=1./a(lj+j) c.... no pivoting do 70 i=l,n if(a(lj+i).eq.0.) go to 70 r=-a(lj+i)*ajji lf=lf0+i ld=j-i c.... do 7 k=l,n c...7 a(i,k)=a(i,k)+r*a(j,k) if(i.gt.nbloc) go to 20 call qma2p(a(lf),a(lf+ld),r,a(lf),nbloc-j,np,np,np) lf=lnb1+i a(lf)=a(lf)+r*a(lf+ld) lf=lnb2+i a(lf)=a(lf)+r*a(lf+ld) go to 40 20 if(i.gt.ntridx) go to 30 a(lf)=a(lf)+r*a(lf+ld) lf=lf+np a(lf)=a(lf)+r*a(lf+ld) lf=lnb1+i a(lf)=a(lf)+r*a(lf+ld) lf=lnb2+i a(lf)=a(lf)+r*a(lf+ld) go to 40 30 call qma2p(a(lf),a(lf+ld),r,a(lf),n-j,np,np,np) 40 b(i)=b(i)+r*b(j) 70 continue lj=lj+np 80 continue c.... use back substitution to find the vector solution npp1=np+1 ln=n1*np li=ln+n b(n)=b(n)/a(li) do 60 l=1,n1 i=n-l li=li-npp1 r=0. lf=ln+i call qvdot1(r,a(lf),b(n),l,0,-np,-1) 60 b(i)=(b(i)-r)/a(li) return end subroutine screen(t,d,abar,zbar,z2bar,ytot1) implicit real*8 (a-h,o-z) save c.... this subroutine calculates screening factors for nuclear reaction c.... rates in the weak, intermediate , and strong regimes. c.... the treatment is based on graboske, dewit, grossman, and cooper c.... ap j. 181,457 (1973) for weak screening and on c.... alastuey and jancovici, ap.j. 226, 1034, 1978, with plasma c.... parameters from itoh, totsuji, setsuo, and dewitt, ap.j. 234, c.... 1079,1979, for strong screening. c.... revised march 22, 1990 to cease screening photodisintegration reactions c.... revised august 10, 1992 to speed up computation c.... c.... last earlier revision 11 nov 1982 . c.... parameter (nscf=78) parameter (nscf1=79) parameter (x13=1.d+0/3.d+0,x14=1.d+0/4.d+0) common/screenc/scfac(nscf) dimension rz1(nscf1),rz2(nscf1),isflag(nscf1),ra1(nscf1), 1 ra2(nscf1),scfac1(nscf1),zhat(nscf1),zhat2(nscf1), 2 gamefac(nscf1),tau12fac(nscf1),cfac(nscf1) data isflag / 1 1,1,1,1,1,1,1,0,1,1, 1,1,1,0,1,0,1,0,1,0, 1,1,1,0,1,0,1,1,1,0, 2 1,0,1,1,1,0,1,0,1,1, 1,0,1,0,1,1,1,0,1,0, 1,1,1,0,1,0,1,1,1,0, 3 1,0,1,1,1,0,1,0,0,0, 0,0,0,0,0,0,1,0,1 / data rz1 / 1 1.,2.,2.,1.,1.,1.,2.,2.,2.,2., 6.,6.,8.,2.,2.,2.,2.,2.,2.,2., 2 2.,2.,1.,1.,2.,2.,2.,2.,1.,1., 2.,2.,2.,2.,1.,1.,2.,2.,2.,2., 3 1.,1.,2.,2.,2.,2.,1.,1.,2.,2., 2.,2.,1.,1.,2.,2.,2.,2.,1.,1., 4 2.,2.,2.,2.,1.,1.,1.,1.,0.,0., 0.,0.,0.,0.,0.,0.,1.,1.,2. / data rz2 / 1 1.,2.,2.,6.,7.,8.,4.,4.,7.,6., 6.,8.,8.,6.,8., 2 8.,10.,10.,12.,12.,12.,12.,13.,13.,14., 14.,14.,14.,15.,15., 3 16.,16.,16.,16.,17.,17.,18.,18.,18.,18., 19.,19.,20.,20.,20., 4 20.,21.,21.,22.,22.,22.,22.,23.,23.,24., 24.,24.,24.,25.,25., 5 26.,26.,26.,26.,27.,27.,26.,26.,26.,26., 26.,26.,2.,2.,1., 6 1.,1.,1.,2. / data ra1 / 1 1.,3.,4.,1.,1.,1.,4.,4.,4.,4., 12.,16.,16.,4.,4., 2 4.,4.,4.,4.,4.,4.,4.,1.,1.,4., 4.,4.,4.,1.,1., 3 4.,4.,4.,4.,1.,1.,4.,4.,4.,4., 1.,1.,4.,4.,4., 4 4.,1.,1.,4.,4.,4.,4.,1.,1.,4., 4.,4.,4.,1.,1., 5 4.,4.,4.,4.,1.,1.,1.,1.,1.,1., 1.,1.,1.,1.,1., 6 1.,1.,1.,4. / data ra2 / 1 1.,3.,3.,12.,14.,16.,8.,8.,14.,12., 12.,12.,16.,12.,16., 2 16.,20.,20.,24.,24.,24.,24.,27.,27.,28., 28.,28.,28.,31.,31., 3 32.,32.,32.,32.,35.,35.,36.,36.,36.,36., 39.,39.,40.,40.,40., 4 40.,43.,43.,44.,44.,44.,44.,47.,47.,48., 48.,48.,48.,51.,51., 5 52.,52.,52.,52.,55.,55.,54.,54.,52.,52., 53.,53.,3.,3.,1., 6 1.,2.,2.,4. / data ipass / 0 / c.... initialize constants if (ipass.gt.0) go to 30 theta=1. x53=5./3. x512=5./12. c.... c.... calculate individual screening factors c.... approx. for strong screening only good for alpha .lt. 1.6 do 110 i1=1,nscf1 z1=rz1(i1) z2=rz2(i1) a1=ra1(i1) a2=ra2(i1) zhat(i1)=(z1+z2)**x53-z1**x53-z2**x53 zhat2(i1)=(z1+z2)**x512-z1**x512-z2**x512 gamefac(i1)=2.**x13*z1*z2/(z1+z2)**x13 tau12fac(i1)=(z1**2*z2**2*a1*a2/(a1+a2))**x13 cfac(i1)=x53*log(z1*z2/(z1+z2)) 110 continue ipass=1 c.... calculate average plasma parameters 30 qlam0=1.88d+8*sqrt(d/(abar*t*t*t)) ztilda=sqrt(z2bar+zbar*theta) qlam0z=qlam0*ztilda gamp=2.27493d+5*(d*zbar*ytot1)**x13/t taufac=4.248719d+3/t**x13 c.... calculate individual screening factors c.... approx. for strong screening only good for alpha .lt. 1.6 do 100 i1=1,nscf1 if(isflag(i1).eq.0) go to 60 gamef=gamp*gamefac(i1) tau12=taufac*tau12fac(i1) alph12=3.*gamef/tau12 c.... c.... limit alph12 to 1.6 to prevent unphysical behavior c.... (h dec. as rho inc.) at high rho. this should really c.... be replaced by a pycnonuclear reaction rate formula. c.... if(alph12.le.1.6) go to 70 alph12=1.6 gamef=1.6*tau12*x13 gamp=gamef/gamefac(i1) 70 h12w=rz1(i1)*rz2(i1)*qlam0z 20 h12=h12w if(gamef.gt.0.3) go to 40 go to 50 40 c=0.896434*gamp*zhat(i1)-3.44740*gamp**x14*zhat2(i1)- 1 0.5551*(log(gamp)+cfac(i1))-2.996 alph122=alph12*alph12 alph123=alph12*alph122 alph124=alph12*alph123 alph125=alph12*alph124 alph126=alph12*alph125 h12=c-(tau12/3.)*(5.*alph123/32.-0.014*alph124 1 -0.0128*alph125)-gamef*(0.0055*alph124 2 -0.0098*alph125+0.0048*alph126) xlgfac=1.-0.0562D0*alph123 if (xlgfac.lt.0.77D0) xlgfac=0.77D0 h12=log(xlgfac)+h12 if(gamef.gt.0.8D0) go to 50 h12=h12w*((0.8D0-gamef)*2.0D0)+h12*((gamef-0.3D0)*2.0D0) 50 if(h12.gt.300.D0) h12=300.D0 if(h12.lt.0.D0) h12=0.D0 scfac1(i1)=exp(h12) go to 100 60 scfac1(i1)=1.D0 100 continue c.... complete triple alpha screening factor scfac1(7)=scfac1(7)*scfac1(nscf1) scfac1(7)=min(scfac1(7),1.0d+130) do i1=1,nscf scfac(i1)=scfac1(i1) ENDDO return end subroutine rate(tt9,rrho,iburn) implicit real*8 (a-h,o-z) save c default for now is ivrate = 1 ivrate = 1 IF (ivrate.EQ.0) THEN c.... rate subroutine by ww95 CALL rate0(tt9,rrho,iburn) ELSE IF (ivrate.EQ.2) THEN c.... rate subroutine with updated and NACRE rates CALL rate2(tt9,rrho,iburn) ELSE c.... rate subroutine updated with rath rates c.... ivrate = 1 Buchmann (1996) C12(a,g) rate c.... ivrate = 3 Kunz et al. 2002 adopted C12(a,g) rate c.... ivrate = 4 Kunz et al. 2002 high C12(a,g) rate c.... ivrate = 5 Kunz et al. 2002 low C12(a,g) rate icagrate=0 IF (ivrate .EQ. 3) icagrate=1 IF (ivrate .EQ. 4) icagrate=2 IF (ivrate .EQ. 5) icagrate=3 CALL rate1(tt9,rrho,iburn,icagrate) ENDIF END c####################################################################### c####################################################################### c####################################################################### c####################################################################### c c.... old rate subroutine subroutine rate0(tt9,rrho,iburn) implicit real*8 (a-h,o-z) save c c this subroutine generates nuclear reaction rate factors c using analytical fitting functions to experimental (fczii) c and theoretical (woosley 76) cross sections. all rates c are current as of may 1976. c.... screening corrections are calculated by subroutine screen c.... and incorporated into the rates in burn. c.... updated march 11, 1988 as per caughlan and fowler and c.... waf (priv. comm.) c.... updated july 10, 1984 to include revisions to rates by c.... harris et al (araa, 21, 165, 1983) and caughlan et al (oap 400) common/ratc/rpp,r33,r34,rcpg,rnpg,ropg,r3a,rg3a,rnag,rcag, 1 r1212,r1216,r1616,roga,roag,rnega,rneag,rmgga,rmgag, 2 rsiga,rmgap,ralpa,ralpg,rsigp,rsiag,rsga,rsiap,rppa,rppg, 3 rsgp,rsag,rarga,rsap,rclpa,rclpg,rargp,rarag,rcaga,rarap, 4 rkpa,rkpg,rcagp,rcaag,rtiga,rcaap,rscpa,rscpg,rtigp, 5 rtiag,rcrga,rtiap,rvpa,rvpg,rcrgp,rcrag,rfega,rcrap, 6 rmnpa,rmnpg,rfegp,rfeag,rniga,rfeap,rcopa,rcopg,rnigp, 7 rfepg,rcogp,r52ng,r53gn,r53ng,r54gn,rheng,rhegn,rhng, 8 rdgn,rdpg,rhegp,rpen,rnep,xlam,sumb,fa,fg,r1,s1,t1,u1, 9 v1,w1,x1,c12agmlt,e1mltc12,e2mltc12,r3amult common /splrat/ ralf1,ralf2,r1f54,r2f54,r3f54,r4f54,r5f54, 1 r6f54,r7f54,r8f54 common /yycom/ yp,yhe3,yalp,yn14,yc12,yo16,yne,ymg,ysi,ys32, 1 ya36,yca,yti,ycr,yfe,yni,y54,yprot,yneut dimension rr(78),rrr(78) equivalence (rr(1),rpp) dimension nrho(52) data nrho / 1,2,3,4,5,6,7,9,10,11,12,13,15,17,19,21,22,23, 1 25,27,28,29,31,33,34,35,37,39,40,41,43,45,46,47,49,51, 2 52,53,55,57,58,59,61,63,64,65,67,69,71,73,75,77 / common /n14pgr/n14pgx c.... if the temperature hasn't changed,just correct the c.... rates for the change in density if(tt9.ne.t9r) go to 50 do i4=1,78 rr(i4)=rrr(i4) enddo rfac=rrho/rho do i1=1,14 i2=nrho(i1) rr(i2)=rr(i2)*rfac enddo r3a=rfac*r3a if(yalp.gt.0.) r34=min(r34,0.896d+0/yalp) if(yp.gt.0.) rnpg=min(rnpg,5.68d-3/yp) if(yp.gt.0.) ropg=min(ropg,1.05d-2/yp) if (iburn.le.0) go to 30 do i1=15,52 i2=nrho(i1) rr(i2)=rr(i2)*rfac enddo if(revx.gt.0.D0)rfepg=rrho*rcogp/revx if(revx.le.0.D0)rfepg=rfac*rfepgx 30 return c.... perform main rate calculation with t9 limited to 10, except c.... in calculating the reverse ratios 50 continue t9r=tt9 t9=t9r if (t9.gt.10.D0) t9=10.D0 rho=rrho t912=sqrt(t9) t913=t9**(1.d+0/3.d+0) t923=t913*t913 t943=t9*t913 t953=t9*t923 t932=t9*t912 t985=t9**1.6D0 t92=t9*t9 t93=t9*t92 t972=t912*t93 t9r32=t9r*sqrt(t9r) t9i=1.D0/t9 t9i13=1.D0/t913 t9i23=1.D0/t923 t9i32=1.D0/t932 t9i12=1.D0/t912 t9ri=1.D0/t9r c c triple alpha reaction rate c rate is for formation of c12 compound nucleus c (has been divided by 6) c c.... rate revised according to caughlan and fowler (nomoto ref.) 1988 c c.... q = -0.092 r2abe=(7.40d+05*t9i32)*exp(-1.0663*t9i) 1 +4.164d+09*t9i23*exp(-13.49*t9i13-(t9/0.098)**2) 2 *(1.+0.031*t913+8.009*t923+1.732*t9+49.883 3 *t943+27.426*t953) c.... q = 7.367 c rbeac=(130.D0*t9i32)*exp(-2.1759D0*t9i) ! 3a res. - 100 keV c rbeac=(130.D0*t9i32)*exp(-4.4969D0*t9i) ! 3a res. + 100 keV c rbeac=(130.D0*t9i32)*exp(-1.0154D0*t9i) ! 3a res. - 200 keV rbeac=(130.D0*t9i32)*exp(-3.3364D0*t9i) 1 +2.510d+07*t9i23*exp(-23.57*t9i13-(t9/0.235)**2) 2 *(1.+0.018*t913+5.249*t923+0.650*t9+19.176 3 *t943+6.034*t953) c.... q = 7.275 if (t9.gt.0.08) then r3a=2.90d-16*(r2abe*rbeac) 1 +0.1*1.35d-07*t9i32*exp(-24.811*t9i) else r3a=2.90d-16*(r2abe*rbeac) 1 *(0.01 + 0.2*(1. + 4.*exp(-(0.025*t9i)**3.263)) 2 /(1. + 4.*exp(-(t9/0.025)**9.227))) 3 +0.1*1.35d-07*t9i32*exp(-24.811*t9i) end if r3a=r3a*(rho*rho)/6. r3a=r3a*r3amult rev=2.00d+20*exp(-84.424*t9ri) rg3a=rev*(t9r**3)*6.*r3a/(rho*rho) c c.... c12 + c12 reaction c.... rate revised according to caughlan and fowler 1988 c.... see cf88 references 47-49 c t9a=t9/(1.+0.0396*t9) t9a13=t9a**(1.d+0/3.d+0) t9a56=t9a**(5.d+0/6.d+0) r1212 = 4.27d+26*t9a56/t932*exp(-84.165/t9a13-2.12d-03*t9**3) r1212 = 0.5*rho*r1212 c c c12 + o16 reaction q = 16.755 c.... rate revised according to caughlan and fowler 1988 c.... this rate only valid for t9 .gt. 0.5 c.... expression for rate from cf88 c c y(nc12)*y(no16)*rc28 is the rate of formation c of the si28 compound nucleus c if (iburn.le.0) go to 20 if (t9.ge.0.5D0) then t9a=t9/(1.+0.055*t9) t9a13=t9a**(1.d+0/3.d+0) t9a23=t9a13*t9a13 t9a56=t9a**(5.d+0/6.d+0) r1216 = 1.72d+31*t9a56*t9i32*exp(-106.594/t9a13) 1 /(exp(-0.18*t9a**2)+1.06d-03*exp(2.562*t9a23)) r1216 = rho*r1216 else r1216 = 0.D0 endif c.... c.... 16o+16o rate updated to cf88 q = 16.542 c.... y16*y16*r32 is rate of formation of 32s compound nucleus c.... r1616 = 7.10d+36*t9i23*exp(-135.93*t9i13-0.629*t923-0.445*t943 1 +0.0103*t9**2) r1616 = rho*r1616*0.5 20 continue c c proton burning (neglect trivial pep) fczii c assume 3p become he3 at rate =rpp c rpp= 0.5*4.01d-15*t9i23*exp(-3.380*t9i13)*(1.+0.123*t913 1 +1.09*t923+0.938*t9)*rho c c he3+he3 goes to he4+2p fczii c r33= 0.5*rho*6.04d+10*t9i23*exp(-12.276*t9i13) 1 *(1.+0.034*t913-0.522*t923-0.124*t9+0.353*t943 2 +0.213*t953) c... c... 3he(ag)7be (listed as he4(he3,g) in cf88) c... q=1.588 c assumed to to to 2 alphas c rate is beta limited (in burn) by decay of b8 c... t9ax=t9/(1.+0.0495*t9) t9a56x=t9ax**(5.d+0/6.d+0) t9am13x=1./(t9ax**0.333333) r34=rho*5.61d+6*t9a56x*t9i32*exp(-12.826*t9am13x) c c 12c(pg)13n assume c12+2p goes to n14 fczii c (no beta limiting, material can go to o14) c term1= 2.04d+07*t9i23*exp(-13.690*t9i13-(t9/1.500)**2) 1 *(1.+0.030*t913+1.19*t923+0.254*t9+2.06*t943+1.12*t953) term2= 1.08d+05*t9i32*exp(-4.925*t9i) 1 +2.15d+05*t9i32*exp(-18.179*t9i) rcpg= rho*(term1+term2) c c 14n(pg)150 this gives rate of formation of (dummy c nucleus) n15. rate is limited (in burn) by the c cf88 beta decay of o14 and o15 (revised 11/27/81) c term1= 4.90d+07*t9i23*exp(-15.228*t9i13-(t9/3.294)**2) 1 *(1.+0.027*t913-0.778*t923-0.149*t9+0.261*t943 2 +0.127*t953) IF (n14pgx.EQ.1) term1=term1*0.54D0 term2= 2.37d+03*t9i32*exp(-3.011*t9i) 1 +2.19d+04*exp(-12.530*t9i) rnpg=rho*(term1+term2) c c a fraction fg of n15 goes (pg) to o16 c a fraction fa of n15 goes (pa) to c12 c c result is n14+2p goes to o16 rate= rnpg*fg c goes to c12+a rate=rnpg*fa c c.... n15(pa) rate revised according to harris et al (fcz iii) c.... terms 1 and 2 are for 15n(pg); terms 3 and 4 are for 15n(pa) term1=9.78d+08*t9i23*exp(-15.251*t9i13-(t9/0.450)**2) 1 *(1.+0.027*t913+0.219*t923+0.042*t9+6.83*t943 2 +3.32*t953) term2= 1.11d+04*t9i32*exp(-3.328*t9i) 1 +1.49d+04*t9i32*exp(-4.665*t9i) 1 +3.80d+06*t9i32*exp(-11.048*t9i) term3=1.08d+12*t9i23*exp(-15.251*t9i13-(t9/0.522)**2) 1 *(1.+0.027*t913+2.62*t923+0.501*t9+5.36*t943+2.60*t953) term4=1.19d+08*t9i32*exp(-3.676*t9i)+5.41d+08*t9i12 1 *exp(-8.926*t9i)+0.1*(4.72d+08*t9i32*exp(-7.721*t9i) 2 +2.20d+09*t9i32*exp(-11.418*t9i)) den= term1+term2+term3+term4 fa= (term3+term4)/den fg= 1.-fa c c 16o(pg)17f o16+2p goes to n14+alpha c rate is beta limited (in burn) by c fczii f17 decay rate c ropg = rho*1.50d+08/(t923*(1.+2.13*(1.-exp(-0.728*t923)))) 1 *exp(-16.692*t9i13) c c 14n(ag)18f n14+1.5*a goes to ne20 c fczii no beta limiting c term1= 7.78d+09*t9i23*exp(-36.031*t9i13-(t9/0.881)**2) 1 *(1.+0.012*t913+1.45*t923+0.117*t9+1.97*t943+0.406*t953) term2= 2.36d-10*t9i32*exp(-2.798*t9i)+2.03*t9i32*exp(-5.054*t9i) 1 +1.15d+04*t9i23*exp(-12.310*t9i) rnag= rho*(term1+term2) c c now do rates for alpha chain c c 12c(ag)16o c cf88 c rcag=e1mltc12*1.04d+08/(t92*(1.+0.0489*t9i23)**2)*exp(-32.120 1 *t9i13-(t9/3.496)**2)+e2mltc12*1.76d+08/(t92*(1.+0.2654*t9i23) 2 **2)*exp(-32.120*t9i13)+1.25d+03*t9i32*exp(-27.499*t9i) 3 +1.43d-02*t92*t93*exp(-15.541*t9i) rcag=rcag*rho*c12agmlt rev=5.13d+10*t9r32*exp(-83.111*t9ri) roga=rev*rcag/rho c c 16o(ag)2one + inverse c c cf88 c term1=9.37d+09*t9i23*exp(-39.757*t9i13 1 -(t9/1.586)**2) term2= 62.1*t9i32*exp(-10.297*t9i) 1 +538.*t9i32*exp(-12.226*t9i)+13.0*t92*exp(-20.093*t9i) roag= rho*(term1+term2) rev= 5.65d+10*t9r32*exp(-54.937*t9ri) rnega= rev*(term1+term2) c c 20ne(ag)24mg + inverse fczii c term1= 4.11d+11*t9i23*exp(-46.766*t9i13-(t9/2.219)**2) 1 *(1.+0.009*t913+0.882*t923+0.055*t9+0.749*t943 2 +0.119*t953) term2= 5.27d+03*t9i32*exp(-15.869*t9i)+6.51d+03*t912 1 *exp(-16.223*t9i) term3=0.1*(42.1*t9i32*exp(-9.115*t9i) 1 +32.0*t9i23*exp(-9.383*t9i)) rneag=rho*(term1+term2+term3)/(1.+5.*exp(-18.960*t9i)) rev= 6.01d+10*t9r32*exp(-108.059*t9ri) rmgga=(rneag/rho)*rev c c.... test for high-z bypass if (iburn.le.0) go to 60 c c 24mg(ag)28si + inverse c revised according to harris et al (fcz iii) c gt9cd =(1.+5.*exp(-15.882*t9i)) rmgag=rho*(4.78d+01*t9i32*exp(-13.506*t9i)+2.38d+03*t9i32 1*exp(-15.218*t9i)+2.47d+02*t932*exp(-15.147*t9i)+0.1*(1.72d-09 2*t9i32*exp(-5.028*t9i)+1.25d-03*t9i32*exp(-7.929*t9i)+2.43d+01*t9i 3*exp(-11.523*t9i)))/gt9cd rev= 6.27d+10*t9r32*exp(-115.862*t9ri) rsiga= rev*rmgag/rho c c c 24mg(ap)27al + inverse cf88 c (given as 27al(pa)24mg) c gt9h=1.+exp(-9.792*t9i)/3.0+2.*exp(-11.773*t9i)/3. term1=1.10d+08*t9i23*exp(-23.261*t9i13-(t9/0.157)**2) 1 *(1.+0.018*t913+12.85*t923+1.61*t9+89.87*t943 2 +28.66*t953) term2=129.*t9i32*exp(-2.517*t9i)+5660.*t972*exp(-3.421*t9i) 1 +0.1*(3.89d-08*t9i32*exp(-0.853*t9i)+8.18d-09*t9i32 2 *exp(-1.001*t9i)) rev=1.81*exp(-18.572*t9ri) rmgap=rev*rho*(term1+term2)/gt9h ralpa=rho*(term1+term2)/gt9h c c 27al(pg)28si cf88 c ralpg=rho*(1.67d+08*t9i23*exp(-23.261*t9i13-(t9/0.155)**2) 1* (1.+0.018*t913+5.81*t923+0.728*t9+27.31*t943+8.71*t953) 2+ 2.20d+00*t9i32*exp(-2.269*t9i) + 1.22d+01*t9i32*exp(-2.491*t9i) 3+ 1.50d+04*t9*exp(-4.112*t9i)+0.1*(6.50d-10*t9i32*exp(-0.853*t9i) 4+ 1.63d-10*t9i32*exp(-1.001*t9i)))/gt9h rev=1.13d+11*t9r32*exp(-134.434*t9ri) rsigp=ralpg*rev/rho c c 28si(ag)32s + inverse whfz78/4 c term=4.82d+22*t9i23*exp(-61.015*t9i13*(1.+6.340d-02*t9 1 +2.541d-03*t92-2.900d-04*t93)) rev= 6.461d+10*t9r32*exp(-80.643*t9ri) rsiag= rho*term rsga= rev*term c c 28si(ap)31p + inverse woosley (76) c (given as 31p(pa)28si) c term=4.16d+13*t9i23*exp(-25.631*t9i13*(1.+2.798d-03*t9 1 +2.763d-03*t92-2.341d-04*t93)) rev= 0.5825*exp(-22.224*t9ri) rppa= rho*term rsiap= rppa*rev c c 31p(pg)32s + inverse woosley (76) c term=1.08d+16*t9i23*exp(-27.042*t9i13*(1.+1.928d-01*t9 1 -1.540d-02*t92+6.444d-04*t93)) rev=3.764d+10*t9r32*exp(-102.865*t9ri) rppg=rho*term rsgp=rev*term c c 32s(ag)36ar + inverse woosley (76) c term=1.16d+24*t9i23*exp(-66.690*t9i13*(1.+4.913d-02*t9 1 +4.637d-03*t92-4.067d-04*t93)) rev= 6.616d+10*t9r32*exp(-77.080*t9ri) rsag= rho*term rarga= rev*term c c 32s(ap)35cl + inverse woosley (76) c (given as 35cl(pa)32s) c term=1.27d+16*t9i23*exp(-31.044*t9i13*(1.+1.041d-01*t9 1 -1.368d-02*t92+6.969d-04*t93)) rev= 1.144*exp(-21.643*t9ri) rsap= rho*rev*term rclpa= rho*term c c 35cl(pg)36ar + inverse woosley (76) c term=4.48d+16*t9i23*exp(-29.483*t9i13*(1.+1.761d-01*t9 1 -1.322d-02*t92+5.245d-04*t93)) rev=7.568d+10*t9r32*exp(-98.722*t9ri) rclpg=rho*term rargp=rev*term c c 36ar(ag)40ca + inverse woosley(76) c term=2.81d+30*t9i23*exp(-78.271*t9i13*(1.+1.458d-01*t9 1 -1.069d-02*t92+3.790d-04*t93)) rev= 6.740d+10*t9r32*exp(-81.711*t9ri) rarag= rho*term rcaga= rev*term c c 36ar(ap)39k + inverse woosley (76) c (given as 39k(pa)36ar) c term=2.76d+13*t9i23*exp(-34.922*t9i13*(1.+4.826d-03*t9 1 -5.534d-03*t92+4.021d-04*t93)) rev= 1.128*exp(-14.959*t9ri) rarap= rho*term*rev rkpa= rho*term c c 39k(pg)40ca + inverse woosley (76) c term=4.09d+16*t9i23*exp(-31.727*t9i13*(1.+1.622d-01*t9 1 -1.119d-02*t92+3.910d-04*t93)) rev=7.600d+10*t9r32*exp(-96.657*t9ri) rkpg=rho*term rcagp=rev*term c c 40ca(ag)44ti + inverse woosley (76) c term=4.66d+24*t9i23*exp(-76.435*t9i13*(1.+1.650d-02*t9 1 +5.973d-03*t92-3.889d-04*t93)) rev= 6.843d+10*t9r32*exp(-59.510*t9ri) rcaag= rho*term rtiga= rev*term c c 40ca(ap)43sc + inverse woosley (76) c (given as 43sc(pa)40ca) c term=4.54d+14*t9i23*exp(-32.177*t9i13*(1.-1.206d-02*t9 1 +7.753d-03*t92-5.071d-04*t93)) rev= 2.229*exp(-40.966*t9ri) rcaap= rho*rev*term rscpa= rho*term c c 43sc(pg)44ti + inverse woosley (76) c term=3.85d+16*t9i23*exp(-33.234*t9i13*(1.+1.023d-01*t9 1 -2.242d-03*t92-5.463d-05*t93)) rev=1.525d+11*t9r32*exp(-100.475*t9ri) rscpg=rho*term rtigp=rev*term c c 44ti(ag)48cr + inverse woosley (76) c term=1.37d+26*t9i23*exp(-81.227*t9i13*(1.+1.066d-01*t9 1 -1.102d-02*t92+5.324d-04*t93)) rev= 6.928d+10*t9r32*exp(-89.289*t9ri) rtiag= rho*term rcrga= rev*term c c 44ti(ap)47v + inverse woosley (76) c (given as 47v(pa)44ti) c term=6.54d+20*t9i23*exp(-66.678*t9i13*(1.+2.655d-02*t9 1 -3.947d-03*t92+2.522d-04*t93)) rev= 1.104*exp(-4.723*t9ri) rtiap= rev*rho*term rvpa= rho*term c c 47v(pg)48cr + inverse woosley (76) c term=2.05d+17*t9i23*exp(-35.568*t9i13*(1.+9.979d-02*t9 1 -2.269d-03*t92-6.662d-05*t93)) rev=7.649d+10*t9r32*exp(-93.999*t9ri) rvpg=rho*term rcrgp=rev*term c c 48cr(ag)52fe + inverse woosley (76) c term=1.04d+23*t9i23*exp(-81.420*t9i13*(1.+6.325d-02*t9 1 -5.671d-03*t92+2.848d-04*t93)) rev= 7.001d+10*t9r32*exp(-92.177*t9ri) rcrag= rho*term rfega= rev*term c c 48cr(ap)51mn + inverse woosley (76) c term=1.83d+26*t9i23*exp(-86.741*t9i13*(1.+1.384d-02*t9 1 +1.081d-03*t92-5.933d-05*t93)) rev= 0.6087*exp(-6.510*t9ri) rcrap= rho*term rmnpa= rev*rho*term c c 51mn(pg)52fe + inverse woosley (76) c term=3.77d+17*t9i23*exp(-37.516*t9i13*(1.+8.922d-02*t9 1 -1.256d-03*t92-9.453d-05*t93)) rev=1.150d+11*t9r32*exp(-85.667*t9ri) rmnpg=rho*term rfegp=rev*term c c 52fe(ag)56ni + inverse woosley (76) c term=1.05d+27*t9i23*exp(-91.674*t9i13*(1.+7.846d-02*t9 1 -7.430d-03*t92+3.723d-04*t93)) rev= 7.064d+10*t9r32*exp(-92.850*t9ri) rfeag= rho*term rniga= rev*term c c 52fe(ap)55co + inverse woosley (76) c term=1.30d+27*t9i23*exp(-91.674*t9i13*(1.+1.367d-02*t9 1 +7.428d-04*t92-3.050d-05*t93)) rev= 0.4597*exp(-9.470*t9ri) rfeap= rho*term rcopa= rev*rho*term c c 55co(pg)56ni + inverse woosley (76) c term=1.21d+18*t9i23*exp(-39.604*t9i13*(1.+9.894d-02*t9 1 -3.131d-03*t92-2.160d-05*t93)) rev= 1.537d+11*t9r32*exp(-83.382*t9ri) rcopg= rho*term rnigp= rev*term c c 54fe(pg)55co + inverse woosley (76) c term=4.51d+17*t9i23*exp(-38.483*t9i13*(1.+9.593d-02*t9 1 -3.445d-03*t92+8.594d-05*t93)) rev= 2.400d+09*t9r32*exp(-58.605*t9ri) rfepg= rho*term rcogp= rev*term revx=rev rfepgx=rfepg c c 52fe(ng)53fe + inverse woosley(76) c tq1=t9/0.348 tq2=t9-0.348 tq10=tq1**0.10 term=9.604d+05*exp(-0.0626*tq2) r52ng=rho*term r53gn=2.43d+09*t9r32*exp(-123.951*t9ri)*term c c 53fe(ng)54fe + inverse woosley(76) c term=1.817d+06*tq10*exp(-0.06319*tq2) r53ng=rho*term r54gn=1.56d+11*t9r32*term*exp(-155.284*t9ri) c c include rates necessary for alpha photodisintegration c c 3he(ng)4he schramm and wagoner (1977) c term=6.62*(1.+905.*t9) rev=2.61d+10*t9r32*exp(-238.81*t9ri) rheng=rho*term rhegn=rev*term c c p(ng)2d schramm and wagoner (1977) c term=4.4d+04*(1.-0.86*t912+0.429*t9) rev=4.71d+09*t9r32*exp(-25.82*t9ri) rhng=rho*term rdgn=rev*term c c 2d(pg)3he cf88 c term=2.24d+03*t9i23*exp(-3.720*t9i13) 1 *(1.+0.112*t913+3.38*t923+2.65*t9) rev=1.63d+10*t9r32*exp(-63.750*t9ri) rdpg=rho*term rhegp=rev*term c c go to 40 c.... zero high-z rates 60 do 80 i=19,78 80 rr(i)=0. r1216=0. r1616=0. c.... remember rates for density extrapolation 40 do 100 i3=1,78 100 rrr(i3)=rr(i3) return end c####################################################################### c####################################################################### c####################################################################### c####################################################################### c c c.... new rate subroutine subroutine rate1(tt9,rrho,iburn,icagrate) c implicit real*8 (a-h,o-z) implicit none save REAL*8 OneThird,TwoThird,OneSixths,FiveSixths PARAMETER ( 1 OneThird =0.333333333333333333333333333333333333333D0, 2 TwoThird =0.666666666666666666666666666666666666667D0, 3 OneSixths =0.166666666666666666666666666666666666667D0, 4 FiveSixths=0.833333333333333333333333333333333333333D0 5 ) REAL*8 MinusOneThird PARAMETER ( 1 MinusOneThird=-OneThird 2 ) c c this subroutine generates nuclear reaction rate factors c using analytical fitting functions to experimental (fczii) c and theoretical (woosley 76) cross sections. all rates c are current as of may 1976. c.... screening corrections are calculated by subroutine screen c.... and incorporated into the rates in burn. c.... Updated 20011215 by Alexander Heger c.... include switch for Kunz et al. (2002) C12(a,g) rates: c.... icagrate = 0 Buchmann (1996) C12(a,g) rate c.... icagrate = 1 Kunz et al. (2002) adopted C12(a,g) rate c.... icagrate = 2 Kunz et al. (2002) high C12(a,g) rate c.... icagrate = 3 Kunz et al. (2002) low C12(a,g) rate c.... updated july, 2000 by T. Rauscher c.... (replaced woosley 1976 by rauscher & thielemann 2000; c.... replaced N=Z alpha capture rates by rath2000 rates c.... renormalized to empirical rates from rauscher, thielemann, c.... goerres, wiescher 2000; c.... replaced schramm & wagoner 1977 by wiescher 1995; c.... updated 12c(a,g) !!) c.... updated march 11, 1988 as per caughlan and fowler and c.... waf (priv. comm.) c.... updated july 10, 1984 to include revisions to rates by c.... harris et al (araa, 21, 165, 1983) and caughlan et al (oap 400) REAl*8 tt9,rrho INTEGER iburn,icagrate REAL*8 rpp,r33,r34,rcpg,rnpg,ropg,r3a,rg3a,rnag,rcag, r1212,r1216 $ ,r1616,roga,roag,rnega,rneag,rmgga,rmgag, rsiga,rmgap,ralpa $ ,ralpg,rsigp,rsiag,rsga,rsiap,rppa,rppg, rsgp,rsag,rarga,rsap $ ,rclpa,rclpg,rargp,rarag,rcaga,rarap, rkpa,rkpg,rcagp,rcaag $ ,rtiga,rcaap,rscpa,rscpg,rtigp, rtiag,rcrga,rtiap,rvpa,rvpg $ ,rcrgp,rcrag,rfega,rcrap, rmnpa,rmnpg,rfegp,rfeag,rniga,rfeap $ ,rcopa,rcopg,rnigp, rfepg,rcogp,r52ng,r53gn,r53ng,r54gn,rheng $ ,rhegn,rhng, rdgn,rdpg,rhegp,rpen,rnep,xlam,sumb,fa,fg,r1,s1 $ ,t1,u1, v1,w1,x1,c12agmlt,e1mltc12,e2mltc12,r3amult REAL*8 ralf1,ralf2,r1f54,r2f54,r3f54,r4f54,r5f54, r6f54,r7f54 $ ,r8f54 REAL*8 yp,yhe3,yalp,yn14,yc12,yo16,yne,ymg,ysi,ys32, ya36,yca,yti $ ,ycr,yfe,yni,y54,yprot,yneut REAL*8 rr,rrr INTEGER nrho REAL*8 t9r,rfac,rho,t9 INTEGER i,i3,i4,i1,i2 REAL*8 term,term1,term2,term3 REAL*8 t912, t913, t923, t943, t953, t932, t985, t92, t93, t972, $ t9r32, t9i, t9i13, t9i23, t9i32, t9i12, t9ri, t9log, t9rlog REAL*8 t9a,t9a13,t9a56, revx,den, rfepgx, r2abe, rbeac, rev, $ t9a23, t9ax,t9a56x, t9am13x, term4, gt9h, $ t9iterm2, rho2 common/ratc/rpp,r33,r34,rcpg,rnpg,ropg,r3a,rg3a,rnag,rcag, 1 r1212,r1216,r1616,roga,roag,rnega,rneag,rmgga,rmgag, 2 rsiga,rmgap,ralpa,ralpg,rsigp,rsiag,rsga,rsiap,rppa,rppg, 3 rsgp,rsag,rarga,rsap,rclpa,rclpg,rargp,rarag,rcaga,rarap, 4 rkpa,rkpg,rcagp,rcaag,rtiga,rcaap,rscpa,rscpg,rtigp, 5 rtiag,rcrga,rtiap,rvpa,rvpg,rcrgp,rcrag,rfega,rcrap, 6 rmnpa,rmnpg,rfegp,rfeag,rniga,rfeap,rcopa,rcopg,rnigp, 7 rfepg,rcogp,r52ng,r53gn,r53ng,r54gn,rheng,rhegn,rhng, 8 rdgn,rdpg,rhegp,rpen,rnep,xlam,sumb,fa,fg,r1,s1,t1,u1, 9 v1,w1,x1,c12agmlt,e1mltc12,e2mltc12,r3amult common /splrat/ ralf1,ralf2,r1f54,r2f54,r3f54,r4f54,r5f54, 1 r6f54,r7f54,r8f54 common /yycom/ yp,yhe3,yalp,yn14,yc12,yo16,yne,ymg,ysi,ys32, 1 ya36,yca,yti,ycr,yfe,yni,y54,yprot,yneut dimension rr(78),rrr(78) equivalence (rr(1),rpp) dimension nrho(52) data nrho / 1,2,3,4,5,6,7,9,10,11,12,13,15,17,19,21,22,23, 1 25,27,28,29,31,33,34,35,37,39,40,41,43,45,46,47,49,51, 2 52,53,55,57,58,59,61,63,64,65,67,69,71,73,75,77 / c.... DATA for Kunz et al (2002) C12(A,G) rate fits REAL*8 acag DIMENSION acag(0:11,3) DATA acag / 1 1.21D8, 6.06D-2, 3.212D+1, 1.7D+0, 7.4D8, 4.70D-1, * 3.212D1, 0.0D0, 0.0D0, 1.53D4, 2.00D6, 3.8534D1, 2 1.35D8, 5.45D-2, 3.212D+1, 1.0D+1, 9.4D8, 4.10D-1, * 3.212D1, 0.0D0, 0.0D0, 1.70D4, 2.22D6, 3.8600D1, 3 3.20D7, 3.50D-2, 3.212D+1, 8.0D-2, 4.6D8, 2.62D-1, * 3.212D1, 0.0D0, 0.0D0, 1.40D4, 1.90D6, 3.8670D1/ integer n14pgx common /n14pgr/n14pgx c.... if the temperature hasn't changed,just correct the c.... rates for the change in density if(tt9.ne.t9r) go to 50 do i4=1,78 rr(i4)=rrr(i4) enddo rfac=rrho/rho do i1=1,14 i2=nrho(i1) rr(i2)=rr(i2)*rfac enddo r3a=rfac*r3a if(yalp.gt.0.D0) r34=min(r34,0.896d+0/yalp) if(yp.gt.0.D0) rnpg=min(rnpg,5.68d-3/yp) if(yp.gt.0.D0) ropg=min(ropg,1.05d-2/yp) if (iburn.le.0) go to 30 do i1=15,52 i2=nrho(i1) rr(i2)=rr(i2)*rfac enddo if(revx.gt.0.D0)rfepg=rrho*rcogp/revx if(revx.le.0.D0)rfepg=rfac*rfepgx 30 return c.... perform main rate calculation with t9 limited to 10, except c.... in calculating the reverse ratios 50 continue t9r=tt9 t9=t9r if (t9.gt.10.D0) t9=10.D0 rho=rrho rho2=rho**2 t912=sqrt(t9) t913=t9**OneThird t923=t913**2 t943=t9*t913 t953=t9*t923 t932=t9*t912 t985=t9**1.6D0 t92=t9**2 t93=t9*t92 t972=t912*t93 t9r32=t9r*sqrt(t9r) t9i=1.D0/t9 t9i13=1.D0/t913 t9i23=1.D0/t923 t9i32=1.D0/t932 t9i12=1.D0/t912 t9ri=1.D0/t9r t9log=log(t9) t9rlog=log(t9r) c c triple alpha reaction rate c rate is for formation of c12 compound nucleus c (has been divided by 6) c c.... rate revised according to caughlan and fowler (nomoto ref.) 1988 c c.... q = -0.092 r2abe=(7.40d+05*t9i32)*exp(-1.0663D0*t9i) 1 +4.164d+09*t9i23*exp(-13.49D0*t9i13-(t9/0.098D0)**2) 2 *(1.+0.031*t913+8.009D0*t923+1.732D0*t9+49.883D0 3 *t943+27.426D0*t953) c.... q = 7.367 c rbeac=(130.D0*t9i32)*exp(-2.1759D0*t9i) ! 3a res. - 100 keV c rbeac=(130.D0*t9i32)*exp(-4.4969D0*t9i) ! 3a res. + 100 keV c rbeac=(130.D0*t9i32)*exp(-1.0154D0*t9i) ! 3a res. - 200 keV rbeac=(130.D0*t9i32)*exp(-3.3364D0*t9i) 1 +2.510d+07*t9i23*exp(-23.57D0*t9i13-(t9/0.235D0)**2) 2 *(1.D0+0.018D0*t913+5.249D0*t923+0.650D0*t9+19.176D0 3 *t943+6.034D0*t953) c.... q = 7.275 if (t9.gt.0.08D0) then r3a=2.90d-16*(r2abe*rbeac) 1 +0.1D0*1.35d-07*t9i32*exp(-24.811D0*t9i) else r3a=2.90d-16*(r2abe*rbeac) 1 *(0.01D0 + 0.2D0*(1.D0 + 4.D0*exp(-(0.025D0*t9i)**3.263D0)) 2 /(1.D0 + 4.D0*exp(-(t9/0.025D0)**9.227D0))) 3 +0.1D0*1.35d-07*t9i32*exp(-24.811D0*t9i) end if r3a=r3a*rho2*OneSixths r3a=r3a*r3amult rev=2.00d+20*exp(-84.424D0*t9ri) rg3a=rev*(t9r**3)*6.D0*r3a/rho2 c c.... c12 + c12 reaction c.... rate revised according to caughlan and fowler 1988 c.... see cf88 references 47-49 c t9a=t9/(1.D0+0.0396D0*t9) t9a13=t9a**OneThird t9a56=t9a**FiveSixths r1212 = 4.27d+26*t9a56/t932*exp(-84.165D0/t9a13-2.12d-03*t93) r1212 = 0.5D0*rho*r1212 c c c12 + o16 reaction q = 16.755 c.... rate revised according to caughlan and fowler 1988 c.... this rate only valid for t9 .gt. 0.5 c.... expression for rate from cf88 c c y(nc12)*y(no16)*rc28 is the rate of formation c of the si28 compound nucleus c if (iburn.le.0) go to 20 if (t9.ge.0.5D0) then t9a=t9/(1.D0+0.055D0*t9) t9a13=t9a**OneThird t9a23=t9a13*t9a13 t9a56=t9a**FiveSixths r1216 = 1.72d+31*t9a56*t9i32*exp(-106.594D0/t9a13) 1 /(exp(-0.18D0*t9a**2)+1.06d-03*exp(2.562D0*t9a23)) r1216 = rho*r1216 else r1216 = 0.D0 endif c.... c.... 16o+16o rate updated to cf88 q = 16.542 c.... y16*y16*r32 is rate of formation of 32s compound nucleus c.... r1616 = 7.10d+36*t9i23*exp(-135.93D0*t9i13-0.629D0*t923 A -0.445D0*t943 1 +0.0103D0*t92) r1616 = rho*r1616*0.5D0 20 continue c c proton burning (neglect trivial pep) fczii c assume 3p become he3 at rate =rpp c rpp= 0.5D0*4.01d-15*t9i23*exp(-3.380D0*t9i13)*(1.D0+0.123D0*t913 1 +1.09D0*t923+0.938D0*t9)*rho c c he3+he3 goes to he4+2p fczii c r33= 0.5D0*rho*6.04d+10*t9i23*exp(-12.276D0*t9i13) 1 *(1.D0+0.034D0*t913-0.522D0*t923-0.124D0*t9+0.353D0*t943 2 +0.213D0*t953) c... c... 3he(ag)7be (listed as he4(he3,g) in cf88) c... q=1.588 c assumed to to to 2 alphas c rate is beta limited (in burn) by decay of b8 c... t9ax=t9/(1.D0+0.0495D0*t9) t9a56x=t9ax**FiveSixths t9am13x=t9ax**MinusOneThird r34=rho*5.61d+6*t9a56x*t9i32*exp(-12.826D0*t9am13x) c c 12c(pg)13n assume c12+2p goes to n14 fczii c (no beta limiting, material can go to o14) c term1= 2.04d+07*t9i23*exp(-13.690D0*t9i13-(t9/1.500D0)**2) *(1.D0 $ +0.030D0*t913+1.19D0*t923+0.254D0*t9+2.06D0*t943+1.12D0*t953) term2= 1.08d+05*t9i32*exp(-4.925D0*t9i) 1 +2.15d+05*t9i32*exp(-18.179D0*t9i) rcpg= rho*(term1+term2) c c 14n(pg)150 this gives rate of formation of (dummy c nucleus) n15. rate is limited (in burn) by the c cf88 beta decay of o14 and o15 (revised 11/27/81) c term1= 4.90d+07*t9i23*exp(-15.228D0*t9i13-(t9/3.294D0)**2) 1 *(1.D0+0.027D0*t913-0.778D0*t923-0.149D0*t9+0.261D0*t943 2 +0.127D0*t953) IF (n14pgx.EQ.1) term1=term1*0.54D0 term2= 2.37d+03*t9i32*exp(-3.011D0*t9i) 1 +2.19d+04*exp(-12.530D0*t9i) rnpg=rho*(term1+term2) c c a fraction fg of n15 goes (pg) to o16 c a fraction fa of n15 goes (pa) to c12 c c result is n14+2p goes to o16 rate= rnpg*fg c goes to c12+a rate=rnpg*fa c c.... n15(pa) rate revised according to harris et al (fcz iii) c.... terms 1 and 2 are for 15n(pg); terms 3 and 4 are for 15n(pa) term1=9.78d+08*t9i23*exp(-15.251D0*t9i13-(t9/0.450D0)**2) 1 *(1.D0+0.027D0*t913+0.219D0*t923+0.042D0*t9+6.83D0*t943 2 +3.32D0*t953) term2= 1.11d+04*t9i32*exp(-3.328D0*t9i) 1 +1.49d+04*t9i32*exp(-4.665D0*t9i) 1 +3.80d+06*t9i32*exp(-11.048D0*t9i) term3=1.08d+12*t9i23*exp(-15.251D0*t9i13-(t9/0.522D0)**2) 1 *(1.D0+0.027D0*t913+2.62D0*t923+0.501D0*t9+5.36*t943 A +2.60D0*t953) term4=1.19d+08*t9i32*exp(-3.676D0*t9i)+5.41d+08*t9i12 1 *exp(-8.926D0*t9i)+0.1D0*(4.72d+08*t9i32*exp(-7.721D0*t9i) 2 +2.20d+09*t9i32*exp(-11.418D0*t9i)) den= term1+term2+term3+term4 fa= (term3+term4)/den fg= 1.D0-fa c c 16o(pg)17f o16+2p goes to n14+alpha c rate is beta limited (in burn) by c fczii f17 decay rate c ropg = rho*1.50d+08/(t923*(1.D0+2.13D0*(1.D0-exp(-0.728D0*t923)))) 1 *exp(-16.692D0*t9i13) c c 14n(ag)18f n14+1.5*a goes to ne20 c fczii no beta limiting c term1= 7.78d+09*t9i23*exp(-36.031D0*t9i13-(t9/0.881D0)**2) 1 *(1.D0+0.012D0*t913+1.45D0*t923+0.117D0*t9+1.97D0*t943 2 +0.406D0*t953) term2= 2.36d-10*t9i32*exp(-2.798D0*t9i) A +2.03D0*t9i32*exp(-5.054D0*t9i) 1 +1.15d+04*t9i23*exp(-12.310D0*t9i) rnag= rho*(term1+term2) c c now do rates for alpha chain c IF (icagrate .EQ. 0) THEN c 12c(ag)16o buchmann (1996) + priv. comm. c (S(300)=146 keV barns) c term1=-3.5738d7/(t92*(1.D0+0.65711d0*t9i23)**2) % *exp(-30.834d0*t9i13-t92/0.40104d0) % +5.0464d8/(t92*(1.D0+0.66456d0*t9i23)**2) % *exp(-31.332d0*t9i13) % -1.6054d0*t9i32*exp(-16.272d0*t9i) % +85028.d0*t9i23*(1.d0+8844.5d0*t913) % *exp(-33.033d0*t9i13) term2=t9i32*(1.45d3*exp(-31.12d0*t9i) % +1.65d4*exp(-37.06d0*t9i) % +5.53d4*exp(-50.56d0*t9i) % +1.07d6*exp(-61.24d0*t9i) % +3.02d6*exp(-68.78d0*t9i)) term3=1.25d+03*t9i32*exp(-27.499d0*t9i) c c limit term1 of the rate to its value at t9=2 c term1=min(term1,3.28d-03) rcag=(term1+term2+term3) ELSE c.... C12(a,g) c.... adopeted, high and low rates by Kunz et al (2002) i=icagrate rcag=acag(0,i)/(t92*(1.D0+acag(1,i)*t9i23)**2* & exp(acag(2,i)*t9i13+(t9i*acag(3,i))**(-2)))+ & acag(4,i)/(t92*(1.D0+acag(5,i)*t9i23)**2* & exp(acag(6,i)*t9i13))+ & acag(7,i)*t9i32*exp(-acag(8,i)*t9i)+ & acag(9,i)*t9i23*(1.D0+acag(10,i)*t913)* & exp(-acag(11,i)*t9i13) ENDIF c.... add in multiplyer, reho-dep., rev. rate rcag=rcag*c12agmlt rev=5.13d+10*t9r32*exp(-83.111d0*t9ri) roga=rev*rcag rcag=rcag*rho c c 16o(ag)2one + inverse c c cf88 c term1=9.37d+09*t9i23*exp(-39.757D0*t9i13 1 -(t9/1.586D0)**2) term2= 62.1D0*t9i32*exp(-10.297D0*t9i) 1 +538.D0*t9i32*exp(-12.226D0*t9i)+13.0D0*t92*exp(-20.093D0*t9i) roag= rho*(term1+term2) rev= 5.65d+10*t9r32*exp(-54.937D0*t9ri) rnega= rev*(term1+term2) c c 20ne(ag)24mg + inverse c if(t9.lt.3.d0) then c c Rauscher, Thielemann, Goerres, Wiescher, NPA 675 (2000) 695 rneag=2.532038433d+04*t9i32* % (2.86d-04*exp(-9.274716d+00*t9i) % +7.52d-03*exp(-1.21724845d+01*t9i) % +7.25d-02*exp(-1.5634256d+01*t9i) % +1.66d-01*exp(-1.5866356d+01*t9i) % +2.45d+00*exp(-1.8628346d+01*t9i) % +1.48d+00*exp(-1.9777241d+01*t9i) % +1.67d+00*exp(-2.2109846d+01*t9i) % +3.06d+00*exp(-2.3421211d+01*t9i) % +3.28d+00*exp(-2.4848626d+01*t9i)) else c c rath(00) renorm (T9=3.) term1=1.333837363d+02+(-2.504361d+00*t9i)+( 7.351683d+01*t9i13) % +(-2.217197d+02*t913)+( 1.314774d+01*t9)+(-7.475602d-01*t953) % +( 9.602703d+01*t9log) rneag=exp(term1) endif rev= 6.013753d+10*t9r32*exp(-1.081075d+02*t9ri) rmgga=rev*rneag rneag=rho*rneag c c.... test for high-z bypass if (iburn.le.0) go to 60 C C 24Mg(ag) + inverse C if(t9.lt.2.d0) then c c Rauscher, Thielemann, Goerres, Wiescher, NPA 675 (2000) 695 rmgag=2.434566168d+04*t9i32* % (2.18d-12*exp(-6.15065d+00*t9i) % +2.17d-12*exp(-6.45238d+00*t9i) % +2.74d-07*exp(-9.527705d+00*t9i) % +4.99d-06*exp(-1.0804255d+01*t9i) % +1.15d-05*exp(-1.1222035d+01*t9i) % +2.95d-05*exp(-1.172105d+01*t9i) % +1.56d-04*exp(-1.269587d+01*t9i) % +4.58d-04*exp(-1.339217d+01*t9i) % +1.10d-01*exp(-1.5214155d+01*t9i) % +2.94d-02*exp(-1.680404d+01*t9i) % +6.00d-02*exp(-1.7767255d+01*t9i) % +4.00d-02*exp(-1.8568d+01*t9i) % +1.40d-01*exp(-1.940356d+01*t9i) % +3.30d-01*exp(-1.9554425d+01*t9i) % +4.32d-01*exp(-2.2594935d+01*t9i) % +3.60d-01*exp(-2.4219635d+01*t9i) % +1.00d+00*exp(-2.850188d+01*t9i) % +3.60d+00*exp(-3.1809305d+01*t9i)) else c c rath(00) renorm (T9=2.) term1=1.428649975d+02+(-3.288633d+00*t9i)+( 1.042707d+02*t9i13) % +(-2.650548d+02*t913)+( 1.391863d+01*t9)+(-6.999523d-01*t953) % +( 1.216164d+02*t9log) rmgag=exp(term1) endif rev= 6.264815d+10*t9r32*exp(-1.158593d+02*t9ri) rsiga=rev*rmgag rmgag=rho*rmgag c c c 24mg(ap)27al + inverse cf88 c (given as 27al(pa)24mg) c gt9h=1.D0+(exp(-9.792D0*t9i)+2.D0*exp(-11.773D0*t9i))/3.D0 term1=1.10d+08*t9i23*exp(-23.261*t9i13-(t9/0.157)**2) 1 *(1.D0+0.018D0*t913+12.85D0*t923+1.61D0*t9+89.87D0*t943 2 +28.66D0*t953) term2=129.D0*t9i32*exp(-2.517D0*t9i) A +5660.D0*t972*exp(-3.421D0*t9i) 1 +0.1D0*(3.89d-08*t9i32*exp(-0.853D0*t9i)+8.18d-09*t9i32 2 *exp(-1.001D0*t9i)) rev=1.81D0*exp(-18.572D0*t9ri) rmgap=rev*rho*(term1+term2)/gt9h ralpa=rho*(term1+term2)/gt9h c c 27al(pg)28si cf88 c ralpg=(1.67d+08*t9i23*exp(-23.261D0*t9i13-(t9/0.155D0)**2) 1* (1.D0+0.018D0*t913+5.81D0*t923 A +0.728D0*t9+27.31D0*t943+8.71D0*t953) 2+ 2.20d+00*t9i32*exp(-2.269D0*t9i) B+ 1.22d+01*t9i32*exp(-2.491D0*t9i) 3+ 1.50d+04*t9*exp(-4.112D0*t9i) C+ 0.1D0*(6.50d-10*t9i32*exp(-0.853D0*t9i) 4+ 1.63d-10*t9i32*exp(-1.001D0*t9i)))/gt9h rev=1.13d+11*t9r32*exp(-134.434D0*t9ri) rsigp=ralpg*rev ralpg=ralpg*rho C C 28Si(ag) + inverse rath(00) renorm C if(t9.lt.3.d0) then c c Rauscher, Thielemann, Goerres, Wiescher, NPA 675 (2000) 695 rsiag=2.3622781d+04*t9i32* % (5.03d-18*exp(-5.628425d+00*t9i) % +9.52d-17*exp(-6.208675d+00*t9i) % +1.64d-12*exp(-8.72696d+00*t9i) % +3.20d-05*exp(-1.6188975d+01*t9i) % +1.60d-02*exp(-1.7929725d+01*t9i) % +1.20d-02*exp(-2.0204305d+01*t9i) % +1.60d-02*exp(-2.218876d+01*t9i) % +5.20d-02*exp(-2.406877d+01*t9i) % +1.52d-01*exp(-2.6076435d+01*t9i) % +5.40d-01*exp(-2.6540635d+01*t9i) % +7.20d-01*exp(-2.9186575d+01*t9i) % +8.30d-01*exp(-2.9441885d+01*t9i) % +6.30d-01*exp(-3.205301d+01*t9i) % +4.50d+01*exp(-4.453999d+01*t9i)) else c c rath(00) renorm (T9=3.) c 9.410800d+01 term1=9.623587d+01+(-3.324446d+00*t9i)+( 5.358524d+01*t9i13) % +(-1.656830d+02*t913)+( 7.199997d+00*t9)+(-2.828443d-01*t953) % +( 7.933873d+01*t9log) rsiag=exp(term1) endif rev= 6.461484d+10*t9r32*exp(-8.062807d+01*t9ri) rsga=rev*rsiag rsiag=rho*rsiag C C 31P(pa) + inverse rath(00) C term1= 1.078204d+02+(-3.592658d+00*t9i)+( 1.339072d+02*t9i13) % +(-2.470038d+02*t913)+( 1.522161d+01*t9)+(-9.210332d-01*t953) % +( 1.171779d+02*t9log) rppa=rho*exp(term1) rev= 5.824570d-01*exp(-2.223422d+01*t9ri) rsiap=rev*rppa C C 31P(pg) + inverse rath(00) C term1= 9.853927d+01+(-3.317683d+00*t9i)+( 1.092536d+02*t9i13) % +(-2.134738d+02*t913)+( 1.299927d+01*t9)+(-7.761410d-01*t953) % +( 9.773750d+01*t9log) rppg=exp(term1) rev= 3.763988d+10*t9r32*exp(-1.028623d+02*t9ri) rsgp=rev*rppg rppg=rho*rppg C C 32S(ag) + inverse rath(00) renorm C term1=-1.915768d+02+( 6.797907d-01*t9i)+(-3.384737d+02*t9i13) % +( 5.501609d+02*t913)+(-3.881261d+01*t9)+( 2.530003d+00*t953) % +(-2.432384d+02*t9log) rsag=exp(term1) rev= 6.616384d+10*t9r32*exp(-7.704228d+01*t9ri) rarga=rev*rsag rsag=rho*rsag C C 35Cl(pa) + inverse rath(00) C term1= 1.033873d+02+(-3.734209d+00*t9i)+( 1.207016d+02*t9i13) % +(-2.314379d+02*t913)+( 1.542642d+01*t9)+(-9.861862d-01*t953) % +( 1.071695d+02*t9log) rclpa=rho*exp(term1) rev= 1.143850d+00*exp(-2.166560d+01*t9ri) rsap=rev*rclpa C C 35Cl(pg) + inverse rath(00) C term1= 1.192145d+02+(-3.909451d+00*t9i)+( 1.359473d+02*t9i13) % +(-2.636043d+02*t913)+( 1.622665d+01*t9)+(-9.716807d-01*t953) % +( 1.202736d+02*t9log) rclpg=exp(term1) rev= 7.568153d+10*t9r32*exp(-9.870788d+01*t9ri) rargp=rev*rclpg rclpg=rho*rclpg C C 36Ar(ag) + inverse rath(00) renorm C term1=-1.289706d+02+( 3.252004d-01*t9i)+(-3.322780d+02*t9i13) % +( 4.687703d+02*t913)+(-2.913671d+01*t9)+( 1.765989d+00*t953) % +(-2.224539d+02*t9log) rarag=exp(term1) rev= 6.740601d+10*t9r32*exp(-8.169568d+01*t9ri) rcaga=rev*rarag rarag=rho*rarag C C 39K(pa) + inverse rath(00) C term1= 7.161177d+01+(-3.724742d+00*t9i)+( 9.635326d+01*t9i13) % +(-1.793475d+02*t913)+( 1.245522d+01*t9)+(-8.604366d-01*t953) % +( 8.715678d+01*t9log) rkpa=rho*exp(term1) rev= 1.127565d+00*exp(-1.494660d+01*t9ri) rarap=rev*rkpa C C 39K(pg) + inverse rath(00) C term1= 1.450013d+02+(-4.520932d+00*t9i)+( 1.660014d+02*t9i13) % +(-3.236890d+02*t913)+( 2.027642d+01*t9)+(-1.221072d+00*t953) % +( 1.467014d+02*t9log) rkpg=exp(term1) rev= 7.600766d+10*t9r32*exp(-9.664228d+01*t9ri) rcagp=rev*rkpg rkpg=rho*rkpg C C 40Ca(ag) + inverse rath(00) renorm C term1=-7.490256d+02+( 1.235853d+01*t9i)+(-1.360082d+03*t9i13) % +( 2.199106d+03*t913)+(-1.330803d+02*t9)+( 7.734556d+00*t953) % +(-1.034036d+03*t9log) rcaag=exp(term1) rev= 6.843156d+10*t9r32*exp(-5.949627d+01*t9ri) rtiga=rev*rcaag rcaag=rho*rcaag C C 43Sc(pa) + inverse rath(00) C term1= 1.359149d+02+(-5.257497d+00*t9i)+( 1.995844d+02*t9i13) % +(-3.470471d+02*t913)+( 1.994009d+01*t9)+(-1.135324d+00*t953) % +( 1.683062d+02*t9log) rscpa=rho*exp(term1) rev= 2.229105d+00*exp(-4.088265d+01*t9ri) rcaap=rev*rscpa C C 43Sc(pg) + inverse rath(00) C term1= 1.648428d+02+(-5.314970d+00*t9i)+( 2.093861d+02*t9i13) % +(-3.894799d+02*t913)+( 2.299542d+01*t9)+(-1.327666d+00*t953) % +( 1.810928d+02*t9log) rscpg=exp(term1) rev= 1.525411d+11*t9r32*exp(-1.003789d+02*t9ri) rtigp=rev*rscpg rscpg=rho*rscpg C C 44Ti(ag) + inverse rath(00) C term1=-9.130837d+02+( 1.594906d+01*t9i)+(-1.694960d+03*t9i13) % +( 2.711843d+03*t913)+(-1.575353d+02*t9)+( 8.856425d+00*t953) % +(-1.292620d+03*t9log) rtiag=exp(term1) rev= 6.928540d+10*t9r32*exp(-8.926181d+01*t9ri) rcrga=rev*rtiag rtiag=rtiag*rho C C 47V(pa) + inverse rath(00) C term1=-3.212816d+02+( 1.836281d+00*t9i)+(-4.304292d+02*t9i13) % +( 7.768205d+02*t913)+(-5.094106d+01*t9)+( 3.056705d+00*t953) % +(-3.355510d+02*t9log) rvpa=rho*exp(term1) rev= 1.103956d+00*exp(-4.734636d+00*t9ri) rtiap=rev*rvpa C C 47V(pg) + inverse rath(00) C term1= 1.946071d+02+(-5.996369d+00*t9i)+( 2.451967d+02*t9i13) % +(-4.601315d+02*t913)+( 2.782109d+01*t9)+(-1.632197d+00*t953) % +( 2.121328d+02*t9log) rvpg=exp(term1) rev= 7.649567d+10*t9r32*exp(-9.399645d+01*t9ri) rcrgp=rev*rvpg rvpg=rvpg*rho C C 48Cr(ag) + inverse rath(00) C term1=-9.291031d+02+( 2.057299d+01*t9i)+(-2.039678d+03*t9i13) % +( 3.077998d+03*t913)+(-1.715707d+02*t9)+( 9.388271d+00*t953) % +(-1.509299d+03*t9log) rcrag=exp(term1) rev= 7.001673d+10*t9r32*exp(-9.212813d+01*t9ri) rfega=rev*rcrag rcrag=rho*rcrag C C 48Cr(ap) + inverse rath(00) C term1=-1.009579d+03+( 1.630895d+01*t9i)+(-1.770721d+03*t9i13) % +( 2.898769d+03*t913)+(-1.729656d+02*t9)+( 9.917646d+00*t953) % +(-1.360558d+03*t9log) rcrap=rho*exp(term1) rev= 6.089616d-01*exp(-6.475311d+00*t9ri) rmnpa=rev*rcrap C C 51Mn(pg) + inverse rath(00) C term1= 2.103038d+02+(-6.679315d+00*t9i)+( 2.786404d+02*t9i13) % +(-5.124387d+02*t913)+( 3.045880d+01*t9)+(-1.768209d+00*t953) % +( 2.386763d+02*t9log) rmnpg=exp(term1) rev= 1.150232d+11*t9r32*exp(-8.565281d+01*t9ri) rfegp=rev*rmnpg rmnpg=rho*rmnpg C C 52Fe(ag) + inverse rath(00) C term1=-1.051633d+03+( 2.255988d+01*t9i)+(-2.240776d+03*t9i13) % +( 3.416331d+03*t913)+(-1.925435d+02*t9)+( 1.063188d+01*t953) % +(-1.666427d+03*t9log) rfeag=exp(term1) rev= 7.064972d+10*t9r32*exp(-9.277798d+01*t9ri) rniga=rev*rfeag rfeag=rho*rfeag C C 52Fe(ap) + inverse rath(00) C term1=-1.059529d+03+( 2.224880d+01*t9i)+(-2.217729d+03*t9i13) % +( 3.406663d+03*t913)+(-1.930014d+02*t9)+( 1.067078d+01*t953) % +(-1.652698d+03*t9log) rfeap=rho*exp(term1) rev= 4.597833d-01*exp(-9.631735d+00*t9ri) rcopa=rev*rfeap C C 55Co(pg) + inverse rath(00) C term1= 2.263128d+02+(-7.208898d+00*t9i)+( 2.993418d+02*t9i13) % +(-5.534111d+02*t913)+( 3.368339d+01*t9)+(-1.991465d+00*t953) % +( 2.561665d+02*t9log) rcopg=exp(term1) rev= 1.536895d+11*t9r32*exp(-8.314624d+01*t9ri) rnigp=rev*rcopg rcopg=rho*rcopg C C 54Fe(pg) + inverse rath(00) C term1= 1.932501d+02+(-6.986446d+00*t9i)+( 2.851875d+02*t9i13) % +(-5.002389d+02*t913)+( 2.907342d+01*t9)+(-1.685965d+00*t953) % +( 2.392321d+02*t9log) rfepg=exp(term1) rev= 2.400157d+09*t9r32*exp(-5.876519d+01*t9ri) rcogp=rev*rfepg rfepg=rho*rfepg revx=rev rfepgx=rfepg C C 52Fe(ng) + inverse rath(00) C term1= 1.222244d+01+(-2.830146d-02*t9i)+( 1.022311d+00*t9i13) % +( 2.649588d+00*t913)+(-2.672443d-01*t9)+( 4.837979d-03*t953) % +(-2.774695d-01*t9log) r52ng=exp(term1) rev= 2.397638d+09*t9r32*exp(-1.239825d+02*t9ri) r53gn=rev*r52ng r52ng=rho*r52ng C C 53Fe(ng) + inverse rath(00) C term1= 1.421715d+01+( 2.013716d-02*t9i)+(-2.135617d+00*t9i13) % +( 3.381802d+00*t913)+( 8.884158d-02*t9)+(-4.017888d-02*t953) % +(-1.887780d+00*t9log) r53ng=exp(term1) rev= 1.535297d+11*t9r32*exp(-1.552450d+02*t9ri) r54gn=rev*r53ng r53ng=rho*r53ng c c include rates necessary for alpha photodisintegration c c c 3He(ng)4He wies (1985) c term1=-0.181429d-10*t9i13 % +0.605380d-10*t913-0.611866d-11*t9+0.514765d-12*t953 term2=-0.205445d-10*t9log rheng=rho*exp(term1+0.138002d+01+0.136916d-12*t9i % +term2) rhegn=t9r32*exp(term1+0.253473d+02-0.238801d+03*t9ri+term2) c c p(ng)2d wies (1985) c term1=-0.202686d-09*t9i13 % +0.690232d-09*t913-0.709921d-10*t9+0.604445d-11*t953 term2=-0.231844d-09*t9log rhng=rho*exp(term1+0.106851d+02+0.150511d-11*t9i % +term2) rdgn=t9r32*exp(term1+0.329456d+02-0.258142d+02*t9ri+term2) c c 2d(pg)3he cf88 c term=2.24d+03*t9i23*exp(-3.720D0*t9i13) 1 *(1.D0+0.112D0*t913+3.38D0*t923+2.65D0*t9) rev=1.63d+10*t9r32*exp(-63.750D0*t9ri) rdpg=rho*term rhegp=rev*term c c go to 40 c.... zero high-z rates 60 continue do i=19,78 rr(i)=0.0D0 enddo r1216=0.D0 r1616=0.D0 c.... remember rates for density extrapolation 40 continue do i3=1,78 rrr(i3)=rr(i3) enddo end c####################################################################### c####################################################################### c####################################################################### c####################################################################### c c c c.... NACRE rate subroutine subroutine rate2(tt9,rrho,iburn) c implicit real*8 (a-h,o-z) implicit none save REAL*8 OneThird,TwoThird,OneSixths,FiveSixths PARAMETER ( 1 OneThird =0.333333333333333333333333333333333333333D0, 2 TwoThird =0.666666666666666666666666666666666666667D0, 3 OneSixths =0.166666666666666666666666666666666666667D0, 4 FiveSixths=0.833333333333333333333333333333333333333D0 5 ) REAL*8 MinusOneThird PARAMETER ( 1 MinusOneThird=-OneThird 2 ) c c this subroutine generates nuclear reaction rate factors c using analytical fitting functions to experimental (fczii) c and theoretical (woosley 76) cross sections. all rates c are current as of may 1976. c.... screening corrections are calculated by subroutine screen c.... and incorporated into the rates in burn. c.... updated july, 2000 by R. Hoffman c.... NACRE rates update CF88 c.... updated july, 2000 by T. Rauscher c.... (replaced woosley 1976 by rauscher & thielemann 2000; c.... replaced N=Z alpha capture rates by rath2000 rates c.... renormalized to empirical rates from rauscher, thielemann, c.... goerres, wiescher 2000; c.... replaced schramm & wagoner 1977 by wiescher 1995; c.... updated 12c(a,g) !!) c.... updated march 11, 1988 as per caughlan and fowler and c.... waf (priv. comm.) c.... updated july 10, 1984 to include revisions to rates by c.... harris et al (araa, 21, 165, 1983) and caughlan et al (oap 400) REAl*8 tt9,rrho INTEGER iburn REAL*8 rpp,r33,r34,rcpg,rnpg,ropg,r3a,rg3a,rnag,rcag, r1212,r1216 $ ,r1616,roga,roag,rnega,rneag,rmgga,rmgag, rsiga,rmgap,ralpa $ ,ralpg,rsigp,rsiag,rsga,rsiap,rppa,rppg, rsgp,rsag,rarga,rsap $ ,rclpa,rclpg,rargp,rarag,rcaga,rarap, rkpa,rkpg,rcagp,rcaag $ ,rtiga,rcaap,rscpa,rscpg,rtigp, rtiag,rcrga,rtiap,rvpa,rvpg $ ,rcrgp,rcrag,rfega,rcrap, rmnpa,rmnpg,rfegp,rfeag,rniga,rfeap $ ,rcopa,rcopg,rnigp, rfepg,rcogp,r52ng,r53gn,r53ng,r54gn,rheng $ ,rhegn,rhng, rdgn,rdpg,rhegp,rpen,rnep,xlam,sumb,fa,fg,r1,s1 $ ,t1,u1, v1,w1,x1,c12agmlt,e1mltc12,e2mltc12,r3amult REAL*8 ralf1,ralf2,r1f54,r2f54,r3f54,r4f54,r5f54, r6f54,r7f54 $ ,r8f54 REAL*8 yp,yhe3,yalp,yn14,yc12,yo16,yne,ymg,ysi,ys32, ya36,yca,yti $ ,ycr,yfe,yni,y54,yprot,yneut REAL*8 rr,rrr INTEGER nrho REAL*8 t9r,rfac,rho,t9 INTEGER i,i3,i4,i1,i2 REAL*8 term,term1,term2,term3 REAL*8 t912, t913, t923, t943, t953, t932, t985, t92, t93, t972, $ t9r32, t9i, t9i13, t9i23, t9i32, t9i12, t9ri, t9log, t9rlog, $ t94, ft9, gt9 REAL*8 t9a,t9a13,t9a56, revx,den, rfepgx, r2abe, rbeac, rev, $ t9a23, t9ax,t9a56x, t9am13x, term4, gt9h, $ t9iterm2, rho2 common/ratc/rpp,r33,r34,rcpg,rnpg,ropg,r3a,rg3a,rnag,rcag, 1 r1212,r1216,r1616,roga,roag,rnega,rneag,rmgga,rmgag, 2 rsiga,rmgap,ralpa,ralpg,rsigp,rsiag,rsga,rsiap,rppa,rppg, 3 rsgp,rsag,rarga,rsap,rclpa,rclpg,rargp,rarag,rcaga,rarap, 4 rkpa,rkpg,rcagp,rcaag,rtiga,rcaap,rscpa,rscpg,rtigp, 5 rtiag,rcrga,rtiap,rvpa,rvpg,rcrgp,rcrag,rfega,rcrap, 6 rmnpa,rmnpg,rfegp,rfeag,rniga,rfeap,rcopa,rcopg,rnigp, 7 rfepg,rcogp,r52ng,r53gn,r53ng,r54gn,rheng,rhegn,rhng, 8 rdgn,rdpg,rhegp,rpen,rnep,xlam,sumb,fa,fg,r1,s1,t1,u1, 9 v1,w1,x1,c12agmlt,e1mltc12,e2mltc12,r3amult common /splrat/ ralf1,ralf2,r1f54,r2f54,r3f54,r4f54,r5f54, 1 r6f54,r7f54,r8f54 common /yycom/ yp,yhe3,yalp,yn14,yc12,yo16,yne,ymg,ysi,ys32, 1 ya36,yca,yti,ycr,yfe,yni,y54,yprot,yneut dimension rr(78),rrr(78) equivalence (rr(1),rpp) dimension nrho(52) data nrho / 1,2,3,4,5,6,7,9,10,11,12,13,15,17,19,21,22,23, 1 25,27,28,29,31,33,34,35,37,39,40,41,43,45,46,47,49,51, 2 52,53,55,57,58,59,61,63,64,65,67,69,71,73,75,77 / c.... if the temperature hasn't changed,just correct the c.... rates for the change in density if(tt9.ne.t9r) go to 50 do 90 i4=1,78 90 rr(i4)=rrr(i4) rfac=rrho/rho do 120 i1=1,14 i2=nrho(i1) rr(i2)=rr(i2)*rfac 120 continue r3a=rfac*r3a if(yalp.gt.0.) r34=min(r34,0.896d+0/yalp) if(yp.gt.0.) rnpg=min(rnpg,5.68d-3/yp) if(yp.gt.0.) ropg=min(ropg,1.05d-2/yp) if (iburn.le.0) go to 30 do 110 i1=15,52 i2=nrho(i1) rr(i2)=rr(i2)*rfac 110 continue if(revx.gt.0.D0)rfepg=rrho*rcogp/revx if(revx.le.0.D0)rfepg=rfac*rfepgx 30 return c.... perform main rate calculation with t9 limited to 10, except c.... in calculating the reverse ratios 50 t9r=tt9 t9=t9r if (t9.gt.10.D0) t9=10.D0 rho=rrho rho2=rho**2 t912=sqrt(t9) t913=t9**OneThird t923=t913**2 t943=t9*t913 t953=t9*t923 t932=t9*t912 t985=t9**1.6D0 t92=t9**2 t93=t9*t92 t94=t9*t93 t972=t912*t93 t9r32=t9r*sqrt(t9r) t9i=1.D0/t9 t9i13=1.D0/t913 t9i23=1.D0/t923 t9i32=1.D0/t932 t9i12=1.D0/t912 t9ri=1.D0/t9r t9log=log(t9) t9rlog=log(t9r) c.... triple alpha reaction rate c.... analytic approximation (23%) c.... (nacre 1999) c rate is for formation of c12 compound nucleus c (has been divided by 6) term1 = 2.43d+09*t9i23*exp(-13.490d0*t9i13-(t9/0.15)**2) term2 = 1.d0+7.45d+01*t9 term3 = 6.09d+05*t9i32*exp(-1.054d0*t9i) r2abe = term1*term2+term3 term1 = 2.76d+07*t9i23*exp(-23.57d0*t9i13-((t9/0.4)**2)) term2 = 1.d0+5.47d0*t9+3.26d+02*t92 term3 = 1.307d+02*t9i32*exp(-3.338d0*t9i) term4 = 2.51d+04*t9i32*exp(-20.307d0*t9i) rbeac = term1*term2+term3+term4 if (t9.le.0.03d0) then r3a = r2abe*rbeac*3.07d-16*(1.d0-29.1d0*t9+1.308d+03*t92) else r3a = r2abe*rbeac*3.44d-16*(1.d0+0.0158d0*(t9**(-0.65d0))) end if r3a = r3a*rho2*OneSixths r3a=r3a*r3amult rev = 2.003d+20*(t9r**3)*exp(-84.415d0*t9i) rg3a = r3a*6.d0*rev/rho2 c c.... c12 + c12 reaction c.... rate revised according to caughlan and fowler 1988 c.... see cf88 references 47-49 c t9a=t9/(1.D0+0.0396D0*t9) t9a13=t9a**OneThird t9a56=t9a**FiveSixths r1212 = 4.27d+26*t9a56/t932*exp(-84.165D0/t9a13-2.12d-03*t93) r1212 = 0.5D0*rho*r1212 c c c12 + o16 reaction q = 16.755 c.... rate revised according to caughlan and fowler 1988 c.... this rate only valid for t9 .gt. 0.5 c.... expression for rate from cf88 c c y(nc12)*y(no16)*rc28 is the rate of formation c of the si28 compound nucleus c if (iburn.le.0) go to 20 if (t9.ge.0.5D0) then t9a=t9/(1.D0+0.055D0*t9) t9a13=t9a**OneThird t9a23=t9a13*t9a13 t9a56=t9a**FiveSixths r1216 = 1.72d+31*t9a56*t9i32*exp(-106.594D0/t9a13) 1 /(exp(-0.18D0*t9a**2)+1.06d-03*exp(2.562D0*t9a23)) r1216 = rho*r1216 else r1216 = 0.D0 endif c.... c.... 16o+16o rate updated to cf88 q = 16.542 c.... y16*y16*r32 is rate of formation of 32s compound nucleus c.... r1616 = 7.10d+36*t9i23*exp(-135.93D0*t9i13-0.629D0*t923 A -0.445D0*t943 1 +0.0103D0*t92) r1616 = rho*r1616*0.5D0 20 continue c c he3+he3 goes to he4+2p fczii c r33= 0.5D0*rho*6.04d+10*t9i23*exp(-12.276D0*t9i13) 1 *(1.D0+0.034D0*t913-0.522D0*t923-0.124D0*t9+0.353D0*t943 2 +0.213D0*t953) c... the fundamental reaction p(p,e+nu)2h c.... forward rate only; divide by two for identical reactants c... q = 1.442 c... c analytical approximation (3%) c (nacre99) rpp=rho*4.08d-15*t9i23*exp(-3.381d0*t9i13)*(1.d0+3.82d0*t9 1 +1.51d0*t92 +1.44d-1*t93-1.14d-2*t94) rpp=rpp*0.5d0 c... 3he(a,g)7be c... q=1.587 c... analytical approximation (19%) min t9 = 0.005 c... (nacre99) c... term=5.46d+6*t9i23*exp(-12.827d0*t9i13) 1 *(1.0d0-0.307d0*t9+8.81d-2*t92-1.06d-2*t93+4.46d-4*t94) r34 =rho*term c... 12c(p,g)13n c... q = 1.943 c... analytical approximation (6%) c... (nacre99) c.. old note c 12c(pg)13n assume c12+2p goes to n14 fczii c (no beta limiting, material can go to o14) term = 2.00d+07*t9i23*exp(-13.692*t9i13-(t9/0.46d0)**2) 1 *(1.0d0+9.89d0*t9-5.98d+01*t92+2.66d+02*t93) 2 +1.00d+05*t9i32*exp(-4.913d0*t9i) 3 +4.24d+05*t9i32*exp(-21.62d0*t9i) rcpg = rho*term c... 14n(p,g)15o cf88 c... q= 7.297 c... analytical approximation (10%) c... (nacre99) c.. old note c c 14n(pg)150 this gives rate of formation of (dummy c nucleus) n15. rate is limited (in burn) by the c cf88 beta decay of o14 and o15 (revised 11/27/81) c term = 4.83d+07*t9i23*exp(-15.231d0*t9i13-(t9/0.8d0)**2) 1 *(1.0d0-2.00d0*t9+3.41d0*t92-2.43d0*t93) 2 +2.36d+03*t9i32*exp(-3.010d0*t9i) 3 +6.72d+03*(t9**0.380d0)*exp(-9.530d0*t9i) rnpg =rho*term c c a fraction fg of n15 goes (pg) to o16 c a fraction fa of n15 goes (pa) to c12 c c result is n14+2p goes to o16 rate= rnpg*fg c goes to c12+a rate=rnpg*fa c c.... n15(pa) rate revised according to (nacre99) c.... terms 1 (and 2) are for 15n(pg); c.... terms 3 (and 4) are for 15n(pa) c... 15n(p,g)16o q = 12.127 analytical approximation (15%) if (t9.le.3.5d0) then term1 = 1.08d+09*t9i23*exp(-15.254d0*t9i13-(t9/0.34d0)**2) 1 *(1.0d0+6.15d0*t9+16.4d0*t92) 2 +9.23d+03*t9i32*exp(-3.597d0*t9i) 3 +3.27d+06*t9i32*exp(-11.024d0*t9i) else term1 = 3.54d+04*(t9**0.095d0)*exp(-2.306d0*t9i) endif c... 15n(p,a)12c q = 4.966 analytical approximation (13%) if (t9.le.2.5d0) then term2 = 1.12d+12*t9i23*exp(-15.253d0*t9i13-(t9/0.28d0)**2) 1 *(1.0d0+4.95d0*t9+1.43d+02*t92) 2 +1.01d+08*t9i32*exp(-3.643d0*t9i) 3 +1.19d+09*t9i32*exp(-7.406d0*t9i) else term2 = 4.17d+07*(t9**0.917d0)*exp(-3.292d0*t9i) endif den = term1+term2 fa = term2/den fg = 1.d0-fa c.. NOTE: no density dependence is included above!!! c c... 16o(p,g)17f c... q = 0.600 c... analytical approximation (10%) c... (nacre99) c old note c c 16o(pg)17f o16+2p goes to n14+alpha c rate is beta limited (in burn) by c fczii f17 decay rate term = 7.37d+07*exp(-16.696d0*t9i13)*(t9**(-0.82d0)) term = term*(1.0d0+2.02d+02*exp(-70.348d0*t9i-0.161d0*t9)) ropg = term*rho c... 14n(a,g)18f c... q = 4.415 c... analytical approximation (5%) c... (nacre99) c.. old note c c 14n(ag)18f n14+1.5*a goes to ne20 c fczii no beta limiting c if (t9.le.2.0d0) then term1 = 2.24d-11*t9i32*exp(-2.750d0*t9i) term2 = 1.68d+00*t9i32*exp(-5.045d0*t9i) term3 = 2.63d+03*(t9**0.364d0)*exp(-10.535d0*t9i) term = (term1+term2+term3) else term = 1.52d+02*(t9**1.567d0)*exp(-6.315d0*t9i) endif term = term*1.0d0-0.340d0*exp(-26.885d0*t9i-0.012d0*t9) rnag = term*rho c c now do rates for alpha chain c c 12c(ag)16o c cf88 c c 12c(ag)16o buchmann (1996) + priv. comm. c (S(300)=146 keV barns) c term1=-3.5738d7/(t92*(1.d0+0.65711d0*t9i23)**2) % *exp(-30.834d0*t9i13-t92/0.40104d0) % +5.0464d8/(t92*(1.d0+0.66456d0*t9i23)**2) % *exp(-31.332d0*t9i13) % -1.6054d0*t9i32*exp(-16.272d0*t9i) % +85028.d0*t9i23*(1.d0+8844.5d0*t913) % *exp(-33.033d0*t9i13) term2=t9i32*(1.45d3*exp(-31.12d0*t9i) % +1.65d4*exp(-37.06d0*t9i) % +5.53d4*exp(-50.56d0*t9i) % +1.07d6*exp(-61.24d0*t9i) % +3.02d6*exp(-68.78d0*t9i)) term3=1.25d+03*t9i32*exp(-27.499d0*t9i) c c limit term1 of the rate to its value at t9=2 c term1=min(term1,3.28d-03) rcag=(term1+term2+term3)*c12agmlt rev=5.13d+10*t9r32*exp(-83.111d0*t9ri) roga=rev*rcag rcag=rcag*rho c... 16o(a,g)20ne c.. q = 4.730 c... analytical approximation (3%) c... (nacre99) term = 2.68d+10*t9i23*exp(-39.760d0*t9i13-(t9/1.6d0)**2) 1 +5.111d+01*t9i32*exp(-10.32d0*t9i) 2 +6.161d+02*t9i32*exp(-12.20d0*t9i) 3 +4.10d-01*(t9**2.966d0)*exp(-11.900d0*t9i) roag = rho*term ft9=1.d0+ 2.38d+00*(t9** 5.15d-01) & *exp(-6.87d+01*t9i+ 1.53d-02*t9) gt9=1.d0+ 1.02d+01*(t9**(-6.02d-01)) & *exp(-1.99d+01*t9i+ 8.95d-02*t9) rnega = ft9/gt9*term*5.653d+10*t932*exp(-54.886d0*t9ri) c... 20ne(a,g)24mg c... q = 9.316 MeV c... analytical approximation (12%) c... nacre99 if (t9.le.1.0d0) then term = 8.72d0*(t9**(-0.532d0))*exp(-8.995d0*t9i) else term = 3.74d+02*(t9**2.229d0)*exp(-12.681d0*t9i) end if term = term*(1.0d0-7.787d0*exp(-19.821d0*t9i-0.114d0*t9)) rneag = rho*term ft9=1.d0+ 1.02d+01*(t9**(-6.02d-01)) & *exp(-1.99d+01*t9i+ 8.95d-02*t9) gt9=1.d0+ 8.23d+00*(t9**(-4.74d-01)) & *exp(-1.65d+01*t9i+ 7.87d-02*t9) rmgga = ft9/gt9*term*6.01d+10*t932*exp(-108.11d0*t9ri) c!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! c.... nacre is older than the following one; suggest to use the more c recent one. (tr) c c 20ne(ag)24mg + inverse c c if(t9.lt.3.d0) then c c Rauscher, Thielemann, Goerres, Wiescher, NPA 675 (2000) 695 c rneag=2.532038433d+04*t9i32* c % (2.86d-04*exp(-9.274716d+00*t9i) c % +7.52d-03*exp(-1.21724845d+01*t9i) c % +7.25d-02*exp(-1.5634256d+01*t9i) c % +1.66d-01*exp(-1.5866356d+01*t9i) c % +2.45d+00*exp(-1.8628346d+01*t9i) c % +1.48d+00*exp(-1.9777241d+01*t9i) c % +1.67d+00*exp(-2.2109846d+01*t9i) c % +3.06d+00*exp(-2.3421211d+01*t9i) c % +3.28d+00*exp(-2.4848626d+01*t9i)) c else c c rath(00) renorm (T9=3.) c term1=1.333837363d+02+(-2.504361d+00*t9i)+( 7.351683d+01*t9i13) c % +(-2.217197d+02*t913)+( 1.314774d+01*t9)+(-7.475602d-01*t953) c % +( 9.602703d+01*t9log) c rneag=exp(term1) c endif c rev= 6.013753d+10*t9r32*exp(-1.081075d+02*t9ri) c rmgga=rev*rneag c rneag=rho*rneag c!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! c.... test for high-z bypass if (iburn.le.0) go to 60 C C 24Mg(ag) + inverse C if(t9.lt.2.d0) then c c Rauscher, Thielemann, Goerres, Wiescher, NPA 675 (2000) 695 rmgag=2.434566168d+04*t9i32* % (2.18d-12*exp(-6.15065d+00*t9i) % +2.17d-12*exp(-6.45238d+00*t9i) % +2.74d-07*exp(-9.527705d+00*t9i) % +4.99d-06*exp(-1.0804255d+01*t9i) % +1.15d-05*exp(-1.1222035d+01*t9i) % +2.95d-05*exp(-1.172105d+01*t9i) % +1.56d-04*exp(-1.269587d+01*t9i) % +4.58d-04*exp(-1.339217d+01*t9i) % +1.10d-01*exp(-1.5214155d+01*t9i) % +2.94d-02*exp(-1.680404d+01*t9i) % +6.00d-02*exp(-1.7767255d+01*t9i) % +4.00d-02*exp(-1.8568d+01*t9i) % +1.40d-01*exp(-1.940356d+01*t9i) % +3.30d-01*exp(-1.9554425d+01*t9i) % +4.32d-01*exp(-2.2594935d+01*t9i) % +3.60d-01*exp(-2.4219635d+01*t9i) % +1.00d+00*exp(-2.850188d+01*t9i) % +3.60d+00*exp(-3.1809305d+01*t9i)) else c c rath(00) renorm (T9=2.) term1=1.428649975d+02+(-3.288633d+00*t9i)+( 1.042707d+02*t9i13) % +(-2.650548d+02*t913)+( 1.391863d+01*t9)+(-6.999523d-01*t953) % +( 1.216164d+02*t9log) rmgag=exp(term1) endif rev= 6.264815d+10*t9r32*exp(-1.158593d+02*t9ri) rsiga=rev*rmgag rmgag=rho*rmgag c c... 24mg(ap)27al + inverse (given as 27al(p,a)24mg) c... q = 1.600 c... analytic approximation (11%) c... nacre99 if (t9.le.6.0d0) then term = 4.76d+10*t9i23*exp(-23.265d0*t9i13-(t9/0.15d0)**2) 1 *(1.0d0-22.3d0*t9+126.7d0*t92) 2 + 9.65d-11*t9i32*exp(-0.834d0*t9i) 3 + 2.09d-03*t9i32*exp(-2.269d0*t9i) 4 + 1.17d-02*t9i32*exp(-3.273d0*t9i) 5 + 2.84d+04*exp(-5.623d0*t9i)*t9i 6 + 1.38d+06*t9*exp(-10.01d0*t9i) else term = 6.02d+05*(t9**1.862d0)*exp(-14.352d0*t9i) endif term = term*(1.0d0+3.471d0*exp(-11.091d0*t9i-0.129d0*t9)) ralpa = rho*term ft9=1.d0+ 9.52d-01*(t9**(-4.51d-01)) & *exp(-1.08d+01*t9i+ 1.69d-01*t9) gt9=1.d0+ 8.23d+00*(t9**(-4.74d-01)) & *exp(-1.65d+01*t9i+ 7.87d-02*t9) rmgap = ft9/gt9*rho*term*1.809d0*exp(-18.574d0*t9ri) c... 27al(p,g)28si c... q = 11.585 c... analytic approximation (10%) c... nacre99 if (t9.le.6.0d0) then term = 2.51d-11*t9i32*exp(-0.839d0*t9i) 1 + 4.82d+01*(t9**(-0.2d0))*exp(-2.223d0*t9i) 2 + 1.76d+03*(t9**1.12d0)*exp(-3.196d0*t9i) 3 + 3.25d+04*(t9**0.251d0)*exp(-5.805d0*t9i) else term = 1.62d+05*(t9**0.549d0)*exp(-17.222d0*t9i) end if term = term*(1.d0-0.669d0*exp(-10.426*t9i-8.0d-03*t9)) ralpg = rho*term ft9=1.d0+ 9.52d-01*(t9**(-4.51d-01)) & *exp(-1.08d+01*t9i+ 1.69d-01*t9) gt9=1.d0+ 9.87d+00*(t9**(-5.50d-01)) & *exp(-2.16d+01*t9i+ 7.73d-02*t9) rsigp =ft9/gt9*term*1.134d+11*t932*exp(-134.44d0*t9ri) C C 28Si(ag) + inverse rath(00) renorm C if(t9.lt.3.d0) then c c Rauscher, Thielemann, Goerres, Wiescher, NPA 675 (2000) 695 rsiag=2.3622781d+04*t9i32* % (5.03d-18*exp(-5.628425d+00*t9i) % +9.52d-17*exp(-6.208675d+00*t9i) % +1.64d-12*exp(-8.72696d+00*t9i) % +3.20d-05*exp(-1.6188975d+01*t9i) % +1.60d-02*exp(-1.7929725d+01*t9i) % +1.20d-02*exp(-2.0204305d+01*t9i) % +1.60d-02*exp(-2.218876d+01*t9i) % +5.20d-02*exp(-2.406877d+01*t9i) % +1.52d-01*exp(-2.6076435d+01*t9i) % +5.40d-01*exp(-2.6540635d+01*t9i) % +7.20d-01*exp(-2.9186575d+01*t9i) % +8.30d-01*exp(-2.9441885d+01*t9i) % +6.30d-01*exp(-3.205301d+01*t9i) % +4.50d+01*exp(-4.453999d+01*t9i)) else c c rath(00) renorm (T9=3.) term1=9.410800d+01+(-3.324446d+00*t9i)+( 5.358524d+01*t9i13) % +(-1.656830d+02*t913)+( 7.199997d+00*t9)+(-2.828443d-01*t953) % +( 7.933873d+01*t9log) rsiag=exp(term1) endif rev= 6.461484d+10*t9r32*exp(-8.062807d+01*t9ri) rsga=rev*rsiag rsiag=rho*rsiag C C 31P(pa) + inverse rath(00) C term1= 1.078204d+02+(-3.592658d+00*t9i)+( 1.339072d+02*t9i13) % +(-2.470038d+02*t913)+( 1.522161d+01*t9)+(-9.210332d-01*t953) % +( 1.171779d+02*t9log) rppa=rho*exp(term1) rev= 5.824570d-01*exp(-2.223422d+01*t9ri) rsiap=rev*rppa C C 31P(pg) + inverse rath(00) C term1= 9.853927d+01+(-3.317683d+00*t9i)+( 1.092536d+02*t9i13) % +(-2.134738d+02*t913)+( 1.299927d+01*t9)+(-7.761410d-01*t953) % +( 9.773750d+01*t9log) rppg=exp(term1) rev= 3.763988d+10*t9r32*exp(-1.028623d+02*t9ri) rsgp=rev*rppg rppg=rho*rppg C C 32S(ag) + inverse rath(00) renorm C term1=-1.915768d+02+( 6.797907d-01*t9i)+(-3.384737d+02*t9i13) % +( 5.501609d+02*t913)+(-3.881261d+01*t9)+( 2.530003d+00*t953) % +(-2.432384d+02*t9log) rsag=exp(term1) rev= 6.616384d+10*t9r32*exp(-7.704228d+01*t9ri) rarga=rev*rsag rsag=rho*rsag C C 35Cl(pa) + inverse rath(00) C term1= 1.033873d+02+(-3.734209d+00*t9i)+( 1.207016d+02*t9i13) % +(-2.314379d+02*t913)+( 1.542642d+01*t9)+(-9.861862d-01*t953) % +( 1.071695d+02*t9log) rclpa=rho*exp(term1) rev= 1.143850d+00*exp(-2.166560d+01*t9ri) rsap=rev*rclpa C C 35Cl(pg) + inverse rath(00) C term1= 1.192145d+02+(-3.909451d+00*t9i)+( 1.359473d+02*t9i13) % +(-2.636043d+02*t913)+( 1.622665d+01*t9)+(-9.716807d-01*t953) % +( 1.202736d+02*t9log) rclpg=exp(term1) rev= 7.568153d+10*t9r32*exp(-9.870788d+01*t9ri) rargp=rev*rclpg rclpg=rho*rclpg C C 36Ar(ag) + inverse rath(00) renorm C term1=-1.289706d+02+( 3.252004d-01*t9i)+(-3.322780d+02*t9i13) % +( 4.687703d+02*t913)+(-2.913671d+01*t9)+( 1.765989d+00*t953) % +(-2.224539d+02*t9log) rarag=exp(term1) rev= 6.740601d+10*t9r32*exp(-8.169568d+01*t9ri) rcaga=rev*rarag rarag=rho*rarag C C 39K(pa) + inverse rath(00) C term1= 7.161177d+01+(-3.724742d+00*t9i)+( 9.635326d+01*t9i13) % +(-1.793475d+02*t913)+( 1.245522d+01*t9)+(-8.604366d-01*t953) % +( 8.715678d+01*t9log) rkpa=rho*exp(term1) rev= 1.127565d+00*exp(-1.494660d+01*t9ri) rarap=rev*rkpa C C 39K(pg) + inverse rath(00) C term1= 1.450013d+02+(-4.520932d+00*t9i)+( 1.660014d+02*t9i13) % +(-3.236890d+02*t913)+( 2.027642d+01*t9)+(-1.221072d+00*t953) % +( 1.467014d+02*t9log) rkpg=exp(term1) rev= 7.600766d+10*t9r32*exp(-9.664228d+01*t9ri) rcagp=rev*rkpg rkpg=rho*rkpg C C 40Ca(ag) + inverse rath(00) renorm C term1=-7.490256d+02+( 1.235853d+01*t9i)+(-1.360082d+03*t9i13) % +( 2.199106d+03*t913)+(-1.330803d+02*t9)+( 7.734556d+00*t953) % +(-1.034036d+03*t9log) rcaag=exp(term1) rev= 6.843156d+10*t9r32*exp(-5.949627d+01*t9ri) rtiga=rev*rcaag rcaag=rho*rcaag C C 43Sc(pa) + inverse rath(00) C term1= 1.359149d+02+(-5.257497d+00*t9i)+( 1.995844d+02*t9i13) % +(-3.470471d+02*t913)+( 1.994009d+01*t9)+(-1.135324d+00*t953) % +( 1.683062d+02*t9log) rscpa=rho*exp(term1) rev= 2.229105d+00*exp(-4.088265d+01*t9ri) rcaap=rev*rscpa C C 43Sc(pg) + inverse rath(00) C term1= 1.648428d+02+(-5.314970d+00*t9i)+( 2.093861d+02*t9i13) % +(-3.894799d+02*t913)+( 2.299542d+01*t9)+(-1.327666d+00*t953) % +( 1.810928d+02*t9log) rscpg=exp(term1) rev= 1.525411d+11*t9r32*exp(-1.003789d+02*t9ri) rtigp=rev*rscpg rscpg=rho*rscpg C C 44Ti(ag) + inverse rath(00) C term1=-9.130837d+02+( 1.594906d+01*t9i)+(-1.694960d+03*t9i13) % +( 2.711843d+03*t913)+(-1.575353d+02*t9)+( 8.856425d+00*t953) % +(-1.292620d+03*t9log) rtiag=exp(term1) rev= 6.928540d+10*t9r32*exp(-8.926181d+01*t9ri) rcrga=rev*rtiag rtiag=rtiag*rho C C 47V(pa) + inverse rath(00) C term1=-3.212816d+02+( 1.836281d+00*t9i)+(-4.304292d+02*t9i13) % +( 7.768205d+02*t913)+(-5.094106d+01*t9)+( 3.056705d+00*t953) % +(-3.355510d+02*t9log) rvpa=rho*exp(term1) rev= 1.103956d+00*exp(-4.734636d+00*t9ri) rtiap=rev*rvpa C C 47V(pg) + inverse rath(00) C term1= 1.946071d+02+(-5.996369d+00*t9i)+( 2.451967d+02*t9i13) % +(-4.601315d+02*t913)+( 2.782109d+01*t9)+(-1.632197d+00*t953) % +( 2.121328d+02*t9log) rvpg=exp(term1) rev= 7.649567d+10*t9r32*exp(-9.399645d+01*t9ri) rcrgp=rev*rvpg rvpg=rvpg*rho C C 48Cr(ag) + inverse rath(00) C term1=-9.291031d+02+( 2.057299d+01*t9i)+(-2.039678d+03*t9i13) % +( 3.077998d+03*t913)+(-1.715707d+02*t9)+( 9.388271d+00*t953) % +(-1.509299d+03*t9log) rcrag=exp(term1) rev= 7.001673d+10*t9r32*exp(-9.212813d+01*t9ri) rfega=rev*rcrag rcrag=rho*rcrag C C 48Cr(ap) + inverse rath(00) C term1=-1.009579d+03+( 1.630895d+01*t9i)+(-1.770721d+03*t9i13) % +( 2.898769d+03*t913)+(-1.729656d+02*t9)+( 9.917646d+00*t953) % +(-1.360558d+03*t9log) rcrap=rho*exp(term1) rev= 6.089616d-01*exp(-6.475311d+00*t9ri) rmnpa=rev*rcrap C C 51Mn(pg) + inverse rath(00) C term1= 2.103038d+02+(-6.679315d+00*t9i)+( 2.786404d+02*t9i13) % +(-5.124387d+02*t913)+( 3.045880d+01*t9)+(-1.768209d+00*t953) % +( 2.386763d+02*t9log) rmnpg=exp(term1) rev= 1.150232d+11*t9r32*exp(-8.565281d+01*t9ri) rfegp=rev*rmnpg rmnpg=rho*rmnpg C C 52Fe(ag) + inverse rath(00) C term1=-1.051633d+03+( 2.255988d+01*t9i)+(-2.240776d+03*t9i13) % +( 3.416331d+03*t913)+(-1.925435d+02*t9)+( 1.063188d+01*t953) % +(-1.666427d+03*t9log) rfeag=exp(term1) rev= 7.064972d+10*t9r32*exp(-9.277798d+01*t9ri) rniga=rev*rfeag rfeag=rho*rfeag C C 52Fe(ap) + inverse rath(00) C term1=-1.059529d+03+( 2.224880d+01*t9i)+(-2.217729d+03*t9i13) % +( 3.406663d+03*t913)+(-1.930014d+02*t9)+( 1.067078d+01*t953) % +(-1.652698d+03*t9log) rfeap=rho*exp(term1) rev= 4.597833d-01*exp(-9.631735d+00*t9ri) rcopa=rev*rfeap C C 55Co(pg) + inverse rath(00) C term1= 2.263128d+02+(-7.208898d+00*t9i)+( 2.993418d+02*t9i13) % +(-5.534111d+02*t913)+( 3.368339d+01*t9)+(-1.991465d+00*t953) % +( 2.561665d+02*t9log) rcopg=exp(term1) rev= 1.536895d+11*t9r32*exp(-8.314624d+01*t9ri) rnigp=rev*rcopg rcopg = rho*rcopg C C 54Fe(pg) + inverse rath(00) C term1= 1.932501d+02+(-6.986446d+00*t9i)+( 2.851875d+02*t9i13) % +(-5.002389d+02*t913)+( 2.907342d+01*t9)+(-1.685965d+00*t953) % +( 2.392321d+02*t9log) rfepg=exp(term1) rev= 2.400157d+09*t9r32*exp(-5.876519d+01*t9ri) rcogp=rev*rfepg rfepg=rho*rfepg revx=rev rfepgx=rfepg C C 52Fe(ng) + inverse rath(00) C term1= 1.222244d+01+(-2.830146d-02*t9i)+( 1.022311d+00*t9i13) % +( 2.649588d+00*t913)+(-2.672443d-01*t9)+( 4.837979d-03*t953) % +(-2.774695d-01*t9log) r52ng=exp(term1) rev= 2.397638d+09*t9r32*exp(-1.239825d+02*t9ri) r53gn=rev*r52ng r52ng=rho*r52ng C C 53Fe(ng) + inverse rath(00) C term1= 1.421715d+01+( 2.013716d-02*t9i)+(-2.135617d+00*t9i13) % +( 3.381802d+00*t913)+( 8.884158d-02*t9)+(-4.017888d-02*t953) % +(-1.887780d+00*t9log) r53ng=exp(term1) rev= 1.535297d+11*t9r32*exp(-1.552450d+02*t9ri) r54gn=rev*r53ng r53ng=rho*r53ng c c include rates necessary for alpha photodisintegration c c c 3He(ng)4He wies (1985) c term1=-0.181429d-10*t9i13 % +0.605380d-10*t913-0.611866d-11*t9+0.514765d-12*t953 term2=-0.205445d-10*t9log rheng=rho*exp(term1+0.138002d+01+0.136916d-12*t9i % +term2) rhegn=t9r32*exp(term1+0.253473d+02-0.238801d+03*t9ri+term2) c c p(ng)2d wies (1985) c term1=-0.202686d-09*t9i13 % +0.690232d-09*t913-0.709921d-10*t9+0.604445d-11*t953 term2=-0.231844d-09*t9log rhng=rho*exp(term1+0.106851d+02+0.150511d-11*t9i % term2) rdgn=t9r32*exp(term1+0.329456d+02-0.258142d+02*t9ri+term2) c c 2d(pg)3he cf88 c term=2.24d+03*t9i23*exp(-3.720D0*t9i13) 1 *(1.D0+0.112D0*t913+3.38D0*t923+2.65D0*t9) rev=1.63d+10*t9r32*exp(-63.750D0*t9ri) rdpg=rho*term rhegp=rev*term c c go to 40 c.... zero high-z rates 60 continue do i=19,78 rr(i)=0.0D0 enddo r1216=0.D0 r1616=0.D0 c.... remember rates for density extrapolation 40 continue do i3=1,78 rrr(i3)=rr(i3) enddo end subroutine qvrota(w,n) implicit real*8 (a-h,o-z) save c.... reverse vector c.... do not use if n greater than 500 dimension w(1),z(500) do 100 k=1,n 100 z(k)=w(n+1-k) do 200 k=1,n 200 w(k)=z(k) return end subroutine qvset(c,w,n) implicit real*8 (a-h,o-z) save dimension w(1) do 100 k=1,n w(k) = c 100 continue return end subroutine qvdot(e,u,v,n) implicit real*8 (a-h,o-z) save dimension u(1),v(1) e=0.d0 do 100 k=1,n e = e + u(k)*v(k) 100 continue return end subroutine qvdot1(e,u,v,n,ie,iu,iv) implicit real*8 (a-h,o-z) save dimension u(0:0),v(0:0) e=0.d0 do 100 k=0,n-1 e = e + u(k*iu)*v(k*iv) 100 continue return end subroutine qma2p(w,u,f,y,n,i1,i2,i3) implicit real*8 (a-h,o-z) save dimension w(1),u(1),y(1) ii1 = 1 ii2 = 1 ii3 = 1 do 100 k=1,n w(ii1) = (u(ii2)*f)+y(ii3) ii1 = ii1 + i1 ii2 = ii2 + i2 ii3 = ii3 + i3 100 continue return end subroutine sneutr(temp_x,den_x,zbar_x,abar_x,total_x,iold_x) implicit real*8 (a-h,o-z) save c defaults for xnumu12 and axion are zero xnumu12 = 0.0d0 axion = 0.0d0 xnumu12_x=xnumu12 axion_x=axion CALL sneutrx( & temp_x, & den_x, & zbar_x, & abar_x, & total_x, & xnumu12_x, & axion_x, & iold_x) end subroutine sneutrx(temp,den,zbar,abar,total,xnumu12,axion,iold) c subroutine sneut_newx(temp,den,zbar,abar, c 1 pair,plas,phot,brem,recomb,pmu,total,wpl) implicit none save c.. c..this routine computes neutrino losses from the analytic fits of c..itoh et al. apjs 102, 411, 1996 c.. c.. 20070116 (alex) c.. add neutrino magnetic moment c.. addition provided by Alexander Friedland and Maurizio Giannotti c.. for now only correction on plasma neutrinos c.. c.. 20070411 (alex) c.. add axion energy loss c.. addition provided by Alexander Friedland and Maurizio Giannotti c.. c..input is the temperature temp, density den, c..mean number of nucleons abar, and mean charge zbar. c..xnumu12 is neutrino magnetic moment in units of 1E-12 Bor magnetrons c..axion is the axion mass in eV c.. c..output is the c..pair neutrino contributions pair, c..plasma neutrino contributions plas, c..photoneutrino contributions phot, c..bremstrahlung neutrino contributions brem, c..recombination neutrino contributions recomb 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 c... addition for nuetrino magnetic moment REAL*8 xnumu12, pmu, wpl, pmupl REAL*8 pmu2, tl75, tl752, tl7510, rl3, rl7 c... additions for axions REAL*8 axion, paxion, zee, fofzee, axcomp 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... some variable for debugging INTEGER ifirst data ifirst/1/ 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.... if iold eq 1 call old subroutine IF (iold.NE.0) THEN t9o=1.D-9*temp rmu=abar/zbar rho=den total=q RETURN ENDIF 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*xl2 + 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..bremsstrahlung neutrino section c..for reactions like e- + (z,a) => e- + (z,a) + nu + nubar c.. n + n => n + n + nu + nubar c.. n + p => n + p + nu + nubar c..equation 4.3 den6 = den * 1.0d-6 tfermi = 5.9302d9*(sqrt(1.0d0+1.018d0*(den6/xmue)**twoth)-1.0d0) c.."weak" degenerate electrons only if (temp .gt. 0.3d0 * tfermi) then t8 = temp * 1.0d-8 t832 = t8 * sqrt(t8) t83 = t8 * t8 * t8 t86 = t83 * t83 t8m2 = 1.0d0/(t8 * t8) t8m5 = t8m2 * t8m2/t8 c.. equation 5.3 eta = rm/(7.05d6 * t832 + 5.12d4 * t83) etam1 = 1.0d0/eta etam2 = etam1 * etam1 c.. equation 5.2 xnum = 1.0d0/(23.5d0 + 6.83d4*t8m2 + 7.81d8*t8m5) xden = 1.26d0*(1.0d0+etam1)/(1.0d0+1.47d0*etam1+3.29d-2*etam2) fbrem = xnum + xden c.. equation 5.9 xnum = 1.0d0/( (1.0d0 + rm*1d-9) 1 * (230.0d0 + 6.7d5*t8m2 + 7666d9*t8m5) ) c00 = 7.75d5*t832 + 247.0d0*t8**(3.85d0) c01 = 4.07d0 + 0.0240d0 * t8**(1.4d0) c02 = 4.59d-5 * t8**(-0.110d0) xden = 1.0d0/(c00/rm + c01 + c02 * den**(0.656d0)) gbrem = xnum + xden c.. equation 5.1 brem = 0.5738d0*zbar*ye*t86*den * (tfac4*fbrem - tfac5*gbrem) c.. liquid metal with c12 parameters (not too different for other elements) c.. equation 5.18 and 5.16 else t8 = temp * 1.0d-8 t86 = t8**6 u = fac3 * (log10(den) - 3.0d0) gamma = 2.275d-1 * zbar * zbar/t8 * (den6/abar)**third gm1 = 1.0d0/gamma gm13 = gm1**third gm23 = gm13**2 c.. equation 5.25 and 5.26 v = -0.05483d0- 0.01946d0*gm13 + 1.86310d0*gm23 - 0.78873d0*gm1 w = -0.06711d0+ 0.06859d0*gm13 + 1.74360d0*gm23 - 0.74498d0*gm1 c.. compute the expensive trig functions of equation 5.21 only once cos1 = cos(u) cos2 = cos(2.0d0*u) cos3 = cos(3.0d0*u) cos4 = cos(4.0d0*u) cos5 = cos(5.0d0*u) sin1 = sin(u) sin2 = sin(2.0d0*u) sin3 = sin(3.0d0*u) sin4 = sin(4.0d0*u) c.. equation 5.21 fb = 0.5d0 * 0.17946d0 + 0.00945d0*u + 0.34529d0 1 - 0.05821d0*cos1 - 0.04969d0*sin1 2 - 0.01089d0*cos2 - 0.01584d0*sin2 3 - 0.01147d0*cos3 - 0.00504d0*sin3 4 - 0.00656d0*cos4 - 0.00281d0*sin4 5 - 0.00519d0*cos5 c.. equation 5.22 ft = 0.5d0 * 0.06781d0 - 0.02342d0*u + 0.24819d0 1 - 0.00944d0*cos1 - 0.02213d0*sin1 2 - 0.01289d0*cos2 - 0.01136d0*sin2 3 - 0.00589d0*cos3 - 0.00467d0*sin3 4 - 0.00404d0*cos4 - 0.00131d0*sin4 5 - 0.00330d0*cos5 c.. equation 5.23 gb = 0.5d0 * 0.00766d0 - 0.01259d0*u + 0.07917d0 1 - 0.00710d0*cos1 + 0.02300d0*sin1 2 - 0.00028d0*cos2 - 0.01078d0*sin2 3 + 0.00232d0*cos3 + 0.00118d0*sin3 4 + 0.00044d0*cos4 - 0.00089d0*sin4 5 + 0.00158d0*cos5 c..equation 5.24 gt = -0.5d0 * 0.00769d0 - 0.00829d0*u + 0.05211d0 1 + 0.00356d0*cos1 + 0.01052d0*sin1 2 - 0.00184d0*cos2 - 0.00354d0*sin2 3 + 0.00146d0*cos3 - 0.00014d0*sin3 4 + 0.00031d0*cos4 - 0.00018d0*sin4 5 + 0.00069d0*cos5 c.. equation 5.19 and 5.20 fliq = v*fb + (1.0d0 - v)*ft gliq = w*gb + (1.0d0 - w)*gt c.. equation 5.17 brem = 0.5738d0*zbar*ye*t86*den * (tfac4*fliq - tfac5*gliq) end if c..equation 6.7 and table 12 zeta = 1.579d5 * zbar * zbar/temp if (nu .ge. -20.0 .and. nu .lt. 0.0) then a1 = 1.51d-2 a2 = 2.42d-1 a3 = 1.21d0 b = 3.71d-2 c = 9.06e-1 d = 9.28d-1 f1 = 0.0d0 f2 = 0.0d0 f3 = 0.0d0 else if (nu .ge. 0.0 .and. nu .le. 10.0) then a1 = 1.23d-2 a2 = 2.66d-1 a3 = 1.30d0 b = 1.17d-1 c = 8.97e-1 d = 1.77d-1 f1 = -1.20d-2 f2 = 2.29d-2 f3 = -1.04d-3 end if c..equation 6.13 and 6.14 if (nu .ge. -20.0 .and. nu .le. 10.0) then z = zeta/(1.0d0 + f1*nu + f2*nu2 + f3*nu3) bigj = (a1/z + a2*z**(-2.25) + a3*z**(-4.55)) * exp(nu) 1 / (1.0d0 + b*exp(c*nu)*(1.0d0 + d*z)) c.. equation 6.5 recomb = tfac6 * 2.649d-18 * ye* zbar**13 * den * bigj 1 / (exp(zeta + nu) + 1.0d0) end if c... compute loss due to neutrino magnetic moment c... (Friedland and Giannotti 20070116) c..Definition of the zero temperature plasma frequency (in K) wpl = 28.7D0 * sqrt(rm) * & (1.0d0+(1.019d-6*rm)**twoth)**(-0.25D0) c.. emission due to magnetic dipole moment c.. plasma neutrinos pmupl=3.18d7 * (xnumu12 / wpl)**2 * qv c.. pair neutrinos c.. limit T range to log10 T = 7.5 .. 10 c.. limit rho range to log10 rho = 3.0 .. 10 tl75=log10(temp)-7.5 IF (tl75 .GT. 2.5D0) tl75=2.5D0 IF (tl75 .LT. 0.0D0) tl75=0.0D0 tl752=tl75**2 tl7510=tl752**5 rl3=log10(den*zbar/abar)-3.D0 IF (rl3 .GT. 7.D0) rl3=7.0D0 IF (rl3 .LT. 0.D0) rl3=0.0D0 rl7=rl3-4.D0 pmu2 = pair * xnumu12**2 & * 1.D-4/(2.959105D0 + 2.151588D-3*rl3**5 & + (3.629246D-1 + 4.66293972D-1*exp(1.0057943D0*rl7))*tl752 & + (9.983571D-3 + 9.00337D-4*exp(1.155642D0*rl7))*tl7510) c.. sum up neutrino contributions pmu = pmupl + pmu2 c... addition for axion cooling c... axion is axion mass in eV c... composition factor (for H/He comp. and for heaveier) if (zbar .LT. 2.D0) then axcomp=1.D0+zbar/abar else axcomp=(1.D0+zbar)*zbar/abar endif axcomp=axcomp*2.D0 zee = 8.28d19*axcomp* den/temp**3 fofzee = (9.02379852065215d0 - & 6.492190469605808d0*Log(zee/(1.d0 + 0.97d0*zee))) & /(2.d0*Pi*(1d0 + 0.15625d0*zee)) c... paxion is in erg/g/s paxion = 4.95d-31 * axcomp*fofzee*(temp**2*axion)**2 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 brem = brem*deni recomb = recomb*deni pmu = pmu*deni c..... temporary test of photo neutrino enhancement c..... needs to be refined later. c$$$ phot = phot * (1.D0 + xnumu12*0.02D0) c$$$ IF ((ifirst .EQ. 1) .AND. (xnumu12 .GT. 0.D0)) THEN c$$$ PRINT*,'[SNEUTR] WARNING: Using fake phot enhancement.' c$$$ ENDIF ifirst=0 c..the total neutrino loss rate total = plas + pair + phot + brem + recomb + pmu + paxion end