subroutine BNUKE(dl,tl,cmp,dxdt,eps,ieqep,epptot, 1 ieqecn,ecno,e3a,eac,enu,ideriv,dlelt,dleld) implicit double precision (a-h,o-z) dimension cmp(15),dxdt(15) dimension dumv(15) c c driver for nuclear reactions c c if ieqep is 1, use equilibrium rate for pp energy generation, c if ieqecn is 1, " " " " CNO " " c if (tl.gt.6.0d0) then call BRATES(dl,tl,cmp,dxdt,eps,ieqep,epptot, 1 ieqecn,ecno,e3a,eac,enu) endif c If energy production rate is smaller than 1.d-8 then assume it is zero c and so are its derivatives. Also do this if the temperature input c is below a NUCLEAR TCUT of, say, 1.d6 if (dabs(eps).lt.1.d-8.or.tl.le.6.0d0) then eps=1.d-99 epptot=eps ecno=eps e3a=eps eac=eps enu=-eps dlelt=0.d0 dleld=0.d0 return endif c c do numerical derivatives if ideriv is equal to 1 if (ideriv.eq.1) then fac=dlog(10.d0) call BRATES(dl+0.0001d0,tl,cmp,dumv,epsd,ieqep,d1, 1 ieqecn,d2,d3,d4,d5) call BRATES(dl,tl+0.0001d0,cmp,dumv,epst,ieqep,d1, 1 ieqecn,d2,d3,d4,d5) dleld=(epsd-eps)/0.0001d0/eps/fac dlelt=(epst-eps)/0.0001d0/eps/fac endif c return end c subroutine BRATES(dl,tl,cmp,dxdt,eps,ieqep,epptot, 1 ieqecn,ecno,e3a,eac,enu) implicit double precision (a-h,o-z) dimension cmp(15),dxdt(15),rate(16),scrn(50) dimension svna(50),q(50),P1(50),amu(15),qeff(50) c Common block CONST contains fundamental physical constants COMMON/CONST/CMS,CRS,CLS,CG,CSB,CPI,CP4,CCL,CME,CH,CNA,CEV,CK,CSYR c c elements of the COMP array; mass fractions! c c # isotope # isotope c --- ------- --- ------- c 1 H1 6 C12 c 2 H2 7 C13 c 3 He3 8 N14 c 4 He4 9 N15 c 5 Li7 10 O16 c c Reactions: c c Proton-Proton CNO c ============= === c 1 - H(p,vb)D 4 - C12(p,g)N13(b,v)C13 c D(p,g)He3 5 - C13(p,g)N14 c 6 - N14(p,g)O15(b,v)N15 c 7 - N15(p,a)C12 c ppI 2 - He3(He3,2p)He4 c He burning c ppII 3 - He3(He4,g)Be7 ========== c 8 - He4(2a,g)C12 c 9 - C12(a,g)O16 c DATA (AMU(i),i=1,10)/1.007825d0,2.014102d0,3.016029d0, 4 4.002603d0,7.016004d0,12.00000d0,13.00335d0, 8 14.00307d0,15.00011d0,15.99492d0/ DATA (P1(I),I=1,9)/2.9645d23,3.3101d22,4.9885d22, 1 4.9795d22,4.5952d22,4.2672d22,3.9835d22, 2 1.5676d21,1.2546d22/ c effective Q values (in MeV) from Bahcall & Ulrich 1988 and CZ88 c (excluding neutrino energies, of course). Note these are all c Qs to completion of the rapid reactions that occasionally follow data (qeff(i),i=1,9)/6.664d0,12.860d0,18.98d0, 4 3.457d0,7.551d0,9.054d0,4.966d0, 8 7.275d0,7.162d0/ c c Compute nuclear reaction rates (per gram per second) c cme=9.1096d-28 cln=dlog(10.d0) camu=1.d0/cna rho=dexp(cln*dl) t9=dexp(cln*tl)/1.d9 t913=dexp(1.d0/3.d0*dlog(t9)) t923=t913*t913 t943=t923*t923 t953=t923*t9 t912=dsqrt(t9) t932=t9*t912 fmue=2.d0/(1.d0+cmp(1)) c c compute cross sections (Na)for various reactions c c PP - from Bahcall 1989, mostly, but also C&Z88 c c p(p,e+nu)d c svna(1)=4.006d-15/t923*dexp(-3.381d0/t913)* 1 (1.d0+0.1232d0*t913+1.09d0*t923+0.938*t9) c c he3(he3,2p)he4 (Bahcall) (PPI) c svna(2)=5.567d10/t923*dexp(-12.276d0/t913)* 1 (1.d0+0.034d0*t913-0.062d0*t923-0.015d0*t9) c c He3(He4,g)Be7 (C&Z 88) c t9a=t9/(1.d0+0.0495d0*t9) t9a13=t9**(0.3333d0) t9a56=t9a**(0.8333d0) svna(3)=5.61d6*t9a56/t932*dexp(-12.826/t9a13) c c CN0 cycle c c C12(p,g)N13(b,v)C13 (Bahcall) c svna(4)=2.11d7/t923*dexp(-13.690d0/t913-(t9/1.500d0)**2)* 1 (1.d0+0.030d0*t913+0.068d0*t923+0.142d0*t9 2 +1.39d0*t943+0.76d0*t953) 3 +1.08d5/t932*dexp(-4.925d0/t9) 4 +2.15d5/t932*dexp(-18.179d0/t9) c c C13(p,g)N14 (FCZ = Bahcall) c svna(5)=8.01e7/t923*dexp(-13.717d0/t913-(t9/2.d0)**2)* 1 (1.d0+0.030d0*t913+0.958d0*t923+0.204d0*t9 2 +1.39d0*t943+0.753d0*t953) 3 +1.35d6/t932*dexp(-5.978d0/t9) 4 +2.66d5/t932*dexp(-11.987d0/t9) 5 +2.26d6/t932*dexp(-13.463d0/t9) c c N14(p,g)O15(b,v)N15 (FCZ = Bahcall) c svna(6)=5.08d7/t923*dexp(-15.228d0/t913-(t9/3.09d0)**2)* 1 (1.d0+0.027d0*t913-0.778d0*t923-0.149d0*t9 2 +0.261d0*t943+0.127d0*t953) 3 +2.28d3/t932*dexp(-3.011d0/t9) 4 +1.65d4*t913*dexp(-12.007d0/t9) c c N15(p,a)C12 c svna(7)=1.191d12/t923*dexp(-15.251d0/t913-(t9/0.139)**2)* 1 (1.d0+0.0273d0*t913+1.971d0*t923+0.372d0*t9) c c Triple alpha ... and beyond c c He4(2a,g)C12 (FCZ) c svna(8)=2.79d-8/T9**3*dexp(-4.4027d0/t9) c c C12(a,g)O16 c svna(9)=9.03d7/T9**2*(1.d0+0.621d0*T923)**2/ 8 (1.d0+0.047d0/T923)**2* 8 dexp(-32.120d0/t913-(T9/5.863d0)**2) c c Reaction rates per gram (= Na X1*X2*rho*Na/a1/a2/(/2 if a1=a2) c rate(1)=svna(1)*rho*cmp(1)*cmp(1)*p1(1) rate(2)=svna(2)*rho*cmp(3)*cmp(3)*p1(2) rate(3)=svna(3)*rho*cmp(3)*cmp(4)*p1(3) rate(4)=svna(4)*rho*cmp(6)*cmp(1)*p1(4) rate(5)=svna(5)*rho*cmp(7)*cmp(1)*p1(5) rate(6)=svna(6)*rho*cmp(8)*cmp(1)*p1(6) rate(7)=svna(7)*rho*cmp(9)*cmp(1)*p1(7) rate(8)=svna(8)*rho*cmp(4)**3*rho*p1(8) rate(9)=svna(9)*rho*cmp(4)*cmp(6)*p1(9) c c find screening corrections call SCREEN(dl,tl,cmp,scrn) c and correct rates with screening corrections do 10 i=1,9 rate(i)=rate(i)*scrn(i) 10 continue c c rate of change of composition=d(mass fraction)/d(time) c Note: the correction for the amount of mass converted into photons c is included here, since these are rates per gram and not c rates per nucleus. c c Another Note: If the rates are almost identical, the dxdt's can and c may be spurious. Set them to zero if they are smaller than c 1 part in 1e12 of the individual rates c c Deuterium is assumed to burn in equilibrium. c Ditto for Li7, Be7, and N15 [rate(6)=rate(7)] c chyd=2.d0*rate(2) bhyd=3.d0*rate(1)+rate(3)+rate(4)+rate(5)+2.d0*rate(6) dxH=chyd-bhyd dxdt(1)=dxH/cna*amu(1) c c dxHe3= rate(1) - 2.0*rate(2) - rate(3) c che3=rate(1) bhe3=2.d0*rate(2)+rate(3) dxhe3=che3-bhe3 if (dabs(dxhe3).lt.1.d-10*dabs(che3+bhe3)) dxhe3=0.d0 dxdt(3)=dxHe3/cna*amu(3) c c note that in PPII, consume 1 He4 but make 2 He4 c note also that rate(7) is effectively set to rate(6) (n15 in equilibrium) c che4=rate(2)+rate(6)+2.d0*rate(3) bhe4=rate(3)+3.d0*rate(8)+rate(9) dxhe4=che4-bhe4 if (dabs(dxhe4).lt.1.d-10*dabs(che4+bhe4)) dxhe4=0.d0 dxdt(4)=dxHe4/cna*amu(4) c c dxC12= rate(7)+rate(8)-rate(4) c note also that rate(7) is effectively set to rate(6) (n15 in equilibrium) c cC12=rate(6)+rate(8) bC12=rate(4) dxC12=cC12-bC12 if (dabs(dxC12).lt.1.d-10*dabs(cC12+bC12)) dxC12=0.d0 dxdt(6)=dxC12/cna*amu(6) c c dxc13= rate(4)-rate(5) c cC13=rate(4) bC13=rate(5) dxC13=cC13-bC13 if (dabs(dxC13).lt.1.d-10*dabs(cC13+bC13)) dxC13=0.d0 dxdt(7)=dxC13/cna*amu(7) c c dxn14= rate(5)-rate(6) c cN14=rate(5) bN14=rate(6) dxN14=cN14-bN14 if (dabs(dxN14).lt.1.d-10*(cN14+bN14)) dxN14=0.d0 dxdt(8)=dxN14/cna*amu(8) c c dxN15=rate(6)-rate(7) c dxN15=0.d0 dxdt(9)=dxN15/cna*amu(9) c c O16 is only produced in C12(ag)016 c dxO16=rate(8) dxdt(10)=dxO16/cna*amu(10) c c Energy production rates c c p-p with deuterium in equilibrium, c and PPII and PPIII from Bahcall; i.e. use the branching ratio c for Be7+e and Be7+p to compute the contributions of PPII and PPIII c in terms of the He3+He4 rate c epsum=0.d0 c epp=cev*1.d6*qeff(1)*rate(1) epp1=cev*1.d6*qeff(2)*rate(2) c epp2lim here is assuming all PPII, no PPIII epp2lim=cev*1.d6*qeff(3)*rate(3) c but compute PPII/PPIII branching ratio, following Bahcall 1989 b7p=3.126571d5*cmp(1)*dexp(-10.26202/t913) b7e=1.752d-10/t912*(1.d0+0.004d0*(1.d3*t9-16.d0)) b7e=b7e/fmue f2=b7e/(b7e+b7p) f3=b7p/(b7e+b7p) c qpp3pp2 is ratio of q for PPIII to q for PPII qpp3pp2=0.689410d0 epp2=cev*1.d6*rate(3)*f2*qeff(3) epp3=cev*1.d6*rate(3)*f3*qeff(3)*qpp3pp2 epptot=epp+epp1+epp2+epp3 c p-p in complete equilibrium (at solar conditions) safecmp=cmp(1)+1.d-5 qppeq=13.094d0*(1.d0+1.412d8*(1.d0/safecmp-1d0)*dexp(-5.d0/t913)) eppeq=cev*1.d6*qppeq*rate(1) if (ieqep.eq.1) epptot=eppeq c CN0 ecno=cev*1.d6*(qeff(4)*rate(4)+qeff(5)*rate(5)+ 1 (qeff(6)+qeff(7))*rate(6)) ecnoeq=cev*1.d6*(qeff(4)+qeff(5)+qeff(6)+qeff(7))* 1 scrn(6)*svna(6)*rho*(cmp(6)+cmp(7)+cmp(8)+cmp(9)) 2 *cmp(1)*p1(6) if (ieqecn.eq.1) then ecno=ecnoeq endif c Helium burning e3a=cev*1.d6*qeff(8)*rate(8) eac=cev*1.d6*qeff(9)*rate(9) ehe=e3a+eac c c neutrino losses c call neutrino(dl,tl,cmp(1),cmp(4),cmp(6),cmp(10), 1 1.d0-cmp(1)-cmp(4), 1 enplas,enph,enpair,enrec,enbr,enu) c total rate (nuclear prodction minus thermal neutrino emission) eps=epptot+ecno+ehe-enu return end c subroutine SCREEN(dl,tl,cmp,scrn) implicit double precision (a-h,o-z) c c Compute electron screening correction a-la c Graboske et al. 1973 (Ap.J. 181,437) c dimension cmp(15),scrn(50),a(10),zq(10),z158(10),fl12(16) dimension zetw(50),zeti(50),zets5(50),zets4(50),zets2(50) c The zet's are the zeta parameters for the reactions, from Table 4. c For strong screening, the sets are broken down by the exponents of the Z's data (zetw(I),I=1,9)/2.d0,8.d0,8.d0,12.d0,12.d0,14.d0,14.d0, 9 16.d0,24.d0/ data (zeti(I),I=1,9)/1.630d0,5.917d0,5.917,8.302d0,8.302d0, 6 9.520d0,9.520d0,11.206d0,16.19d0/ data (zets5(i),I=1,9)/1.18d0,3.73d0,3.73d0,4.81d0,4.81d0, 6 5.38d0,5.38d0,6.56d0,9.01d0/ data (zets4(i),I=1,9)/0.52d0,1.31d0,1.31d0,1.49d0,1.49d0, 6 1.62d0,1.62d0,2.03d0,2.58d0/ data (zets2(i),i=1,9)/-0.41d0,-0.65d0,-0.65d0,-0.64d0, 6 -0.64d0,-0.66d0,-0.66d0,-0.81d0,-0.89d0/ c zq is the charge on species i data (zq(i),I=1,10)/1.d0,1.d0,2.d0,2.d0,3.d0,6.d0, 7 6.d0,7.d0,7.d0,8.d0/ c z158 is the charge on species i to the 1.58 power data (z158(i),I=1,10)/1.d0,1.d0,2.9897d0,2.9897d0,5.6735d0, 6 16.962d0,16.962d0,21.640d0,21.640d0, 9 26.723d0/ c a is the nuclear mass of species i, in amu data a/1.007825d0,2.014102d0,3.016029d0,4.002603d0,7.016929d0, 6 12.00000d0,13.03354d0,14.00307d0,15.00011d0,15.99491d0/ c c************************ T=dexp(2.302585093d0*tl) rho=dexp(2.302585093d0*dl) c Now some quantities defined by Graboske and Salpeter'54) ztwid2=0.d0 zbar=0.d0 z2bar=0.d0 z158b=0.d0 oomui=0.d0 do 10 i=1,10 if (cmp(i).lt.1.d-8) go to 10 zbar=zbar+zq(i)*cmp(i)/a(i) z2bar=z2bar+zq(i)*zq(i)*cmp(i)/a(i) z158b=z158b+z158(i)*cmp(i)/a(i) oomui=oomui+cmp(i)/a(i) 10 continue fmui=1.d0/oomui c ztwid computed in the nondegenerate gas case ztwid2=z2bar+zbar ztwid=dsqrt(ztwid2) fl0=1.879d8*dsqrt(rho/T**3*oomui) c compute some logs of these quantities, and some others, here c to speed things up third=1.d0/3.d0 dlfl0=dlog(fl0) dlzbar=dlog(zbar) dlztw=dlog(ztwid) c for intermediate screening etab=z158b/(dexp(0.58d0*dlztw+0.28d0*dlzbar)) sifac=0.380d0*etab*dexp(0.86d0*dlfl0) c strong screening flag - set ssf4 .lt. 0 ssf4=-1.d0 c code below will reset it if strong screening needed c c Now compute screening for each reaction rate: do 30 i=1,9 c Compute the quantity lambda(12) for each reaction c Note: fl12 is the weak screening exponential factor, and also c the discriminator between weak and intermediate, and c intermediate and strong, screening. fl12(i)=zetw(i)/2.d0*ztwid*fl0 if (fl12(i).le.0.1d0) then c weak screening scrn(i)=dexp(fl12(i)) c If beyond weak screening, but not into strong screening, c compute intermediate screening value elseif (fl12(i).le.5.d0) then c Intermediate screening: sinter=sifac*zeti(i) endif c compute strong screening if fl12 > 2 if (fl12(i).gt.2.d0) then c constants for strong screening; compute 'em once if (ssf4.lt.0.d0) then ssf4=0.316d0*dexp(third*dlzbar) ssf2=0.737d0/zbar*dexp(-2.d0*third*dlfl0) ssf0=0.624d0*dexp(third*dlzbar)*dexp(2.d0*third*dlfl0) endif zetab=zets5(i)+ssf4*zets4(i)+ssf2*zets2(i) strong=ssf0*zetab endif c Now, must choose intermediate if 0.1 < fl12 < 2.0, strong if >5, c and between the two between 2 and 5 if (fl12(i).le.2.d0.and.fl12(i).gt.0.1d0) scrn(i)=dexp(sinter) if (fl12(i).ge.5.d0) scrn(i)=dexp(strong) if (fl12(i).gt.2.d0.and.fl12(i).lt.5.d0) then scrn(i)=dexp(dmin1(sinter,strong)) endif 30 continue return end c subroutine CNEQ(dl,tl,dtime,cmp) IMPLICIT DOUBLE PRECISION (A-H,O-Z) C c Takes input mass fractions of C12,C13,N14,N15 and c computes equilibrium abundances from CN burning given the c values of svna from XSECT, screened by SCREEN c c N15 is set equal to its equilibrium abundance, while the others c are compute on approach to equilibrium c c Conserves total mass fraction of C12, C13, N14 and N15 C DIMENSION svna(50),q(50),P1(50),scrn(50),cmp(15),amu(15) c Common block CONST contains fundamental physical constants COMMON/CONST/CMS,CRS,CLS,CG,CSB,CPI,CP4,CCL,CME,CH,CNA,CEV,CK,CSYR DATA (P1(I),I=1,16)/2.9645d23,2.9667d23,2.9668d23,3.3101d22, 5 4.9885d22,1.5645d26,8.5167d22,8.5156d22, 9 1.9812d23,9.9339d22,4.9795d22,4.5952d22, 3 4.2672d22,3.9835d22,1.5676d21,1.2546d22/ DATA AMU/1.007825d0,2.014102d0,3.016029d0,4.002603d0, 5 6.015125d0,7.016004d0,7.016929d0,12.00000d0, 9 13.00335d0,14.00307d0,15.00011d0,15.99492d0, 3 1.000000d0,1.000000d0,0.0005458d0/ c c Be sure to conserve total mass fraction of material! c xc12=cmp(8) xc13=cmp(9) xn14=cmp(10) xn15=cmp(11) xsumin=xc12+xc13+xn14+xn15 c Get reaction rates... c call XSECT(TL,SVNA,Q) call SCREEN(dl,tl,cmp,scrn) C Now calculate the EQUILIBRIUM abundances with respect to N14: c The statements below are true when the CN cycle is running in complete c equilibrium, with no contribution from the ON reactions fac=svna(13)*scrn(13)*p1(13) x12ox14=fac/(svna(11)*scrn(11)*p1(11)) x13ox14=fac/(svna(12)*scrn(12)*p1(12)) x15ox14=fac/(svna(14)*scrn(14)*p1(14)) xn14e=xsumin/(1.d0+x12ox14+x13ox14+x15ox14) xc12e=x12ox14*xn14e xc13e=x13ox14*xn14e xn15e=x15ox14*xn14e c The equilibrium abundances retain the total mass fraction of the C c and N isotopes, but now in equilibrium ratios c c The time scales for approach to equilibrium ratios: t12m1=10.d0**dl/cna*p1(11)*svna(11)*scrn(11)*amu(8)*cmp(1) t12=1.d0/t12m1 t13m1=10.d0**dl/cna*p1(12)*svna(12)*scrn(12)*amu(9)*cmp(1) t13=1.d0/t13m1 t14m1=10.d0**dl/cna*p1(13)*svna(13)*scrn(13)*amu(10)*cmp(1) t14=1.d0/t14m1 c c print '(" Time scales for CN:",1p3e12.5," yr")',t12/3.1556d7, c 1 t13/3.1556d7,t14/3.1556d7 c Compute the values after time dtime of approach to equilibrium c values from the input values - if the exponential factors get too c large, keep them from becoming obnoxious by setting to something innocuous fac12=dtime/t12 if (fac12.gt.100.d0) fac12=100.d0 fac13=dtime/t13 if (fac13.gt.100.d0) fac13=100.d0 fac14=dtime/t14 if (fac14.gt.100.d0) fac14=100.d0 c fac15=dtime/t15 c if (fac15.gt.100.d0) fac15=100.d0 xc12=xc12e-(xc12e-xc12)*dexp(-fac12) xc13=xc13e-(xc13e-xc13)*dexp(-fac13) xn14=xn14e-(xn14e-xn14)*dexp(-fac14) xn15=xn15e c c Make sure sum adds up to xsumin xsum=xc12+xc13+xn14+xn15 cmp(8)=xc12*xsumin/xsum cmp(9)=xc13*xsumin/xsum cmp(10)=xn14*xsumin/xsum cmp(11)=xn15*xsumin/xsum c RETURN END C subroutine PPEQ(dl,tl,dtime,cmp) IMPLICIT DOUBLE PRECISION (A-H,O-Z) C c Takes input mass fractions of hydrogen, and He4 and c computes equilibrium abundances of deuterium and He3 c and beryllium c from PPI burning given the values of svna from XSECT. e.g. Clayton, p. 369 c DIMENSION svna(50),q(50),P1(50),cmp(15),scrn(50) c Common block CONST contains fundamental physical constants COMMON/CONST/CMS,CRS,CLS,CG,CSB,CPI,CP4,CCL,CME,CH,CNA,CEV,CK,CSYR DATA (P1(I),I=1,16)/2.9645d23,2.9667d23,2.9668d23,3.3101d22, 5 4.9885d22,1.5645d26,8.5167d22,8.5156d22, 9 1.9812d23,9.9339d22,4.9795d22,4.5952d22, 3 4.2672d22,3.9835d22,1.5676d21,1.2546d22/ c c Be sure to conserve total mass fraction of material! c rho=10.d0**dl t9=10.d0**(tl-9.d0) x=cmp(1) xh2=cmp(2) xhe3=cmp(3) xhe4=cmp(4) xtot=x+xh2 ytot=xhe3+xhe4 c Get reaction rates... c call XSECT(TL,SVNA,Q) call SCREEN(dl,tl,cmp,scrn) c C NOW CALCULATE THE EQUILIBRIUM ABUNDANCES: c c Deuterium: want p(p,e+)d - d(p,g)he3 to be zero fac=p1(1)*svna(1)*scrn(1)/(p1(3)*svna(3)*scrn(3)) xh2=fac*xtot/(1.d0+fac) x=xtot-xh2 cmp(1)=x cmp(2)=xh2 c c Helium 3: want p(d,g)He3 - He3(He3,2p)He4/2 to be zero c but with deuterium in equilibrium, this is the same c as p(p,e+)d - He3(He3,2p)He4/2 equalling zero c NOTE: THIS ASSUMPTION IS ONLY VALID IF T is SUFFICIENTLY HIGH c FOR HE3 to come into EQUILIBRIUM - See Clayton, p.377 c c Equilibrium value of He3: fac=dsqrt(p1(1)*svna(1)*scrn(1)/(2.d0*p1(4)*svna(4)*scrn(4))) xhe30=X*fac c Lifetime for He3 against destruction upon reaching this equilibrium temp=2.d0*p1(4)*svna(4)*scrn(4)*3.01602945d0 t33=1.d0/(rho/cna*xhe30*temp) c If input value of He3 is less than equilibrium, compute the c value of He3 upon approach to equilibrium from zero c (Clayton, 5-24) over timescale dtime c If dtime/t33 gets too big, set it to something innocuous exfac33=dtime/t33 if (exfac33.gt.100.d0) exfac33=100.d0 if (xhe3.lt.xhe30) then fac=(xhe30-xhe3)/(xhe30+xhe3) xhe3=xhe30*(dexp(2.d0*exfac33)-fac)/ 1 (dexp(2.d0*exfac33)+fac) else c If input value of He3 is greater than equilibrium fac=(xhe3-xhe30)/(xhe3+xhe30) xhe3=xhe30*(1.d0+dexp(-2.d0*exfac33)*fac)/ 1 (1.d0-dexp(-2.d0*exfac33)*fac) endif c Convert helium abundances to mass fractions, conserving total helium c but taking from He4 to get He3 cmp(3)=xhe3 xhe4=ytot-xhe3 cmp(4)=xhe4 c Be7: fmue=2.d0/(1.d0+cmp(1)) fac=(p1(5)*svna(5)*scrn(5))/ 1 (svna(8)*scrn(8)*p1(8)) xbe7=fac*cmp(3)*cmp(4)/cmp(1) c conserve composition; call the source of the Be7 He4 cmp(4)=cmp(4)-xbe7 cmp(7)=xbe7 c Lithium 7 - assume Be7 in equilibrium (rate(6) fast): fac=(svna(5)*scrn(5)*p1(5))/(svna(7)*scrn(7)*p1(7)) xli7=fac*cmp(3)*cmp(4)/cmp(1) c call the source of the Li7 He4 cmp(4)=cmp(4)-xli7 cmp(6)=xli7 c RETURN END