program bzams implicit double precision (a-h,o-z) parameter(jmax=1999,kmaxx=2000,nmax=50) c c A program to integrate a homogeneous model of a star, c c Common block STRUC contains model structure common/STRUC/h(jmax,30),comp(jmax,15),sr(jmax),nj,modno c Common block GLOB contains global model parameters: xat0 is the initial X common/GLOB/fms,xat0,z,alfa,time,dtime,pcl,tcl,fls,teffl c Common block CONST contains fundamental physical constants COMMON/CONST/CMS,CRS,CLS,CG,CSB,CPI,CP4,CCL,CME,CH,CNA,CEV,CK,CSYR common/TORHS/cmp(15),sp,st,ishell c Common block EQRATE contains info for equilibrium hydrogen burning common/EQRATE/ieqep,ieqecn dimension v1(2),v2(2),delv1(2),delv2(2),f(4),dv1(2),dv2(2) dimension v(4),ucorr(4) dimension dydx(4),y2(4),y1(4),fmconv(10),iconv(10) dimension y1err(4),y2err(4) c common block carrying name of opacity directory character*40 opacdir common/opacdirectory/opacdir c Common block for shooting routines common /caller/x1,x2,xf,nvar,nn2 common /path/ kmax,kount,dxsav,xp(KMAXX),yp(NMAX,KMAXX) equivalence (v1(1),v(1)),(v2(1),v(3)) logical check external RHS data comp/29985*0.d0/ kmax=0 c open (3,file='bzams.sum') open (4,file='bzams.out') open (20,file='bzams.puls') open (30,file='first.mod') c c Call initialization routine START that reads in fundamental physical c constants into the block /CONST/, and the parameters for the model. c It returns the requested fit point, and loads the rest in common c blocks. NOTE: when parameters used as arguments in START, creepy c things happen with the central point. I dunno. c call START(fmfit,done) c modno=1 c srfit=fms-fmfit nj=jmax c Set first and last mesh points as 0.999M* and 1.995d-10M* srbot=0.999d0*fms c srtop=5.d-10*fms srtop=5.d-9*fms c The helium abundance Y=1.d0-xat0-z xCNO=0.5556*z c Load the compostion array. do 10 i=1,nj comp(i,1)=xat0 comp(i,2)=0.d0 c Set He3 to .001*Hetot comp(i,3)=0.0010d0*Y comp(i,4)=0.9990d0*Y c CNO relative abundances from Ross-Aller(1976), as quoted c in Bahcall and Ulrich (1982), with c12/c13=90. if (fms.le.1.5d0) then comp(i,6)= xCNO*0.289d0 comp(i,8)=xCNO*0.070d0 comp(i,10)=xCNO*0.641d0 c Isotope abundances for Carbon as quoted by Boothroyd and Sackmann (1990) comp(i,7)=comp(i,8)/90.d0 comp(i,6)=comp(i,8)*(1.d0-1.d0/90.d0) else c if CNO burning in the star, then use quasi-CNO Equilibrium values comp(i,6)=xCNO*0.00686d0 comp(i,7)=xCNO*0.00214d0 comp(i,8)=xCNO*0.9899d0 comp(i,9)=0.d0 comp(i,10)=xCNO*0.0011d0 endif 10 continue c c Set up the shooting algorithm: initial guesses at surface and central c properties from input values of logL,logTe,logPc, and logTc c v1(1)=dlog10(fls) v1(2)=teffl v2(1)=pcl v2(2)=tcl c And shoot away! x1=srtop x2=srbot xf=srfit nvar=4 nn2=2 ucorr(1)=0.90d0 ucorr(2)=0.90d0 ucorr(3)=0.90d0 ucorr(4)=0.90d0 if (done.lt.1.d0) then call NEWT(v,4,check,done) if (check) then print *,' shootf failed!!!' stop endif endif c c Final values obtained! Now the fun part. Integrate the model c with these final values and save them! c flsl=v1(1) fls=10.d0**v1(1) teffl=v1(2) pcl=v2(1) pc=10.d0**pcl tcl=v2(2) tc=10.d0**tcl print 102,flsl,teffl,pcl,tcl 102 format('Final Values: log(L/Lo) =',f10.5/ 1 ' log Te =',f10.5/ 2 ' log Pc =',f10.5/ 3 ' log Tc =',f10.5) c c compute central conditions c call LOAD2(srbot,v2,y2) c c Compute conditions the next step out using central expansions c computed by LOAD2. NOTE: LOAD2 must be called BEFORE FILHC, c as FILHC alters the abundances in the COMP array at the zone. c Doing it the other way around results in effective burning to c complete equilibrium. c sr(1)=srbot flnR=y2(1) flnP=y2(2) flnT=y2(3) flr=y2(4) call FILHC(1,flnR,flnP,flnT,flr,sr(1)) c c now integrate outwards to the fitting point, saving things as you go c i=1 do while (sr(i).gt.srfit) call RHS(sr(i),y2,dydx) c calculate the maximum step allowed so that the principal variables c don't change by more than the tolerances. dsr=dabs(0.10d0/(dydx(1)+1.d-128)) dsp=dabs(0.15d0/(dydx(2)+1.d-128)) dst=dabs(0.10d0/(dydx(3)+1.d-128)) if (dabs(y2(4)).gt.0.01*fls) then dsl=dabs(0.10d0*y2(4)/(dydx(4)+1.d-128)) else dsl=100.d0 endif ds=-dmin1(dsr,dsp,dst,dsl) temp=sr(i)+ds if (temp.lt.srfit) ds=srfit-sr(i) call rkck(y2,dydx,4,sr(i),ds,y2,y2err,rhs) i=i+1 sr(i)=sr(i-1)+ds flnR=y2(1) flnP=y2(2) flnT=y2(3) flr=y2(4) c load the composition array with the appropriate values obtained by c RHS ishell=0 do 101 jj=1,15 comp(i,jj)=cmp(jj) 101 continue call FILHC(i,flnR,flnP,flnT,flr,sr(i)) enddo print '(a,3f11.6,1pe13.5)',' Inside-out R,P,T,L',flnR/dlog(10.d0), 1 flnP/dlog(10.d0),flnT/dlog(10.d0),flr ninter=i c c Now integrate down from the surface to the fitting point, saving c things as you go. c First use the atmosphere integrator to obtain conditions at the c photosphere. c fl=10.d0**flsl teff=10.d0**teffl reff=dsqrt(cls*fl/cp4/csb/teff**4)/crs do 15 i=1,15 cmp(i)=comp(nj,i) 15 continue call LOAD1(srtop,v1,y1) i=nj sr(nj)=srtop c Save the topmost point in the model flnr=y1(1) flnp=y1(2) flnt=y1(3) flr=y1(4) call FILHC(nj,flnr,flnp,flnt,flr,srtop) c and integrate downwards do while (sr(i).lt.srfit) call RHS(sr(i),y1,dydx) c Calculate the maximum step allowed so that the principal variables c don't change by more than the tolerances. dsr=dabs(0.10d0/(dydx(1)+1.d-128)) dsp=dabs(0.15d0/(dydx(2)+1.d-128)) dst=dabs(0.10d0/(dydx(3)+1.d-128)) dsl=dabs(0.10d0*y1(4)/(dydx(4)+1.d-128)) ds=dmin1(dsr,dsp,dst,dsl) temp=sr(i)+ds if (temp.gt.srfit) ds=srfit-sr(i) c and make the downwards step call rkck(y1,dydx,4,sr(i),ds,y1,y1err,RHS) i=i-1 sr(i)=sr(i+1)+ds flnR=y1(1) flnP=y1(2) flnT=y1(3) flr=y1(4) c load the composition array with the appropriate values as determined by c RHS do 201 jj=1,15 comp(i,jj)=cmp(jj) 201 continue call FILHC(i,flnR,flnP,flnT,flr,sr(i)) enddo print '(a,3f11.6,1pe13.5)',' Outside-in R,P,T,L',flnR/dlog(10.d0), 1 flnP/dlog(10.d0),flnT/dlog(10.d0),flr nenv=nj-i c c Shift down the COMP, SR, and H arrays so that the envelope c lies on top of the core, and figure out the total number of mesh c points. c do 42 i=1,nenv sr(ninter+i)=sr(nj-nenv+i) do 43 ii=1,30 h(ninter+i,ii)=h(nj-nenv+i,ii) 43 continue do 44 ii=1,15 comp(ninter+i,ii)=comp(nj-nenv+i,ii) 44 continue 42 continue nj=ninter+nenv c c Find convection boundaries and such. c call CONVEC(fmconv,iconv,nconv) do 62 i=1,nconv print 162, i,iconv(i),fmconv(i)/fms write (4,162) i,iconv(i),fmconv(i)/fms 162 format(' Convection zone number',i3,' Near zone no.',i4, 1 ' at Fractional mass:',1pe12.5) 62 continue if (nconv.eq.0) write (4,*) ' No convective zones!' c c And print out model details c write (3,163) fms,flsl,teffl,dlog10(reff),pcl,tcl, 1 h(1,5)/dlog(10.d0),fmconv(1)/fms,fmconv(2)/fms,xat0,z,alfa 163 format(f6.3,f8.4,f7.4,2f8.4,f7.4,4f6.3,f7.4,f5.2) c fac=1.d0/dlog(10.d0) c write (4,142) write (4,152) 152 format(" Zone 1-Mr/Msun lg(r/Rs) log P log T L/Ls" 1 ," log d log K log eps Del DelAd DelRad" 2 ," X Y ") write (4,142) do 52 i=1,nj write(4,142) i,sr(i),fac*h(i,1),fac*h(i,2),fac*h(i,3), 1 h(i,4),fac*h(i,5),dlog10(h(i,12)),dlog10(h(i,15)), 1 h(i,8),h(i,6),h(i,7),comp(i,1),comp(i,4) 142 format(i4,1pe11.3,0p,f8.3,f9.3,f8.3,1pe11.3,0p5f8.3,f12.3,2f9.5) 52 continue c c Call the routine that writes a model out in a format that can be c read by pulsation codes... c call PULSOUT c c And the routine that writes out the model in binary form for c use in the evolution code. call MODOUT c stop end c subroutine START(fmfit,done) implicit double precision (a-h,o-z) c c Initialization routine that reads in fundamental physical c constants into the block /CONST/, atomic data for EOS, and the c parameters for the model.from file ZAMS.PAR C c Common block GLOB contains global model parameters: xat0 is the initial X common/GLOB/fms,xat0,z,alfa,time,dtime,pcl,tcl,fls,teffl c Common block EQRATE contains info for equilibrium hydrogen burning common/EQRATE/ieqep,ieqecn c Common block CTRL contains most of the control parameters for passing common/CTRL/dllenv,dltenv,tcut,dtfac,sstep(7),tstep(7), 1 nmod,nsav,nprint,nzp,itstep 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 SET PHYSICAL AND ASTRONOMICAL CONSTANTS C CMS=1.989D33 CRS=6.9598D10 CLS=3.86D33 CG=6.672320D-8 CSB=5.67051D-5 CPI=3.14159265358979D0 CP4=CPI*4.D0 CCL=2.99792458D10 CME=9.1093897d-28 CH=6.6260755D-27 CNA=6.022136D23 CEV=1.60217733D-12 CK=1.380658D-16 CSYR=3.15576d7 c c Call subroutine that loads atomic data for EOS (complete EOS) c c call RATDAT c c Read in properties of the model from ZAMS.PAR c open (2,file='bzams.par',status='OLD') c print *, ' Enter mass of model in solar masses' read (2,*) fms c print *, ' Enter X, Z' read (2,*) xat0,z c print *,' Enter (1,1) for equilibrium pp,CN energy generation' read (2,*) ieqep,ieqecn c print *,' Enter desired mixing length' read (2,*) alfa c print *,' Enter guessed log(L/Lo), log(Te)' read (2,*) flsl,teffl fls=10.d0**flsl c print *,' Enter guessed central log pressure, log temperature' read (2,*) pcl,tcl c print *,' Enter mass at fitting point' read (2,*) fmfit c print *,' Enter quality of fit at fit point' read (2,*) done c should also read in tcut and zoning parameters, but thats an unnecessary c whistle right now. Just set tcut=5.5 for now. (12/28/93) tcut=5.5 c Set dtime to be the thermal time scale r=dsqrt(10.d0**(flsl-4.d0*teffl)*cls/(cp4*csb))/crs dtime=2.d7*3.1556d7*fms**2/10.d0**flsl/r c return end c subroutine CONVEC(fmconv,iconv,nconv) implicit double precision (a-h,o-z) parameter(jmax=1999) c c Finds convection zone boundaries, working from center to surface c Takes as 'Input' the structure of a model in the H array c As output, produces an array FMCONV with masses of convection boundaries c and an array ICONV with the next exterior mass zone # c for the transition. c c Common block STRUC contains model structure common/STRUC/h(jmax,30),comp(jmax,15),sr(jmax),nj,modno c Common block GLOB contains global model parameters: xat0 is the initial X common/GLOB/fms,xat0,z,alfa,time,dtime,pcl,tcl,fls,teffl c dimension fmconv(10),iconv(10) c c Note: H(i,6) is the adiabatic gradient, c while H(i,7) is the radiative gradient c nconv=0 temp0=h(1,6)-h(1,7) do 10 i=2,nj temp1=h(i,6)-h(i,7) test=temp1*temp0 if (test.lt.0.d0) then c Found a convective/radiative boundary! Interpolate to find precise c mass at which the transition occurs nconv=nconv+1 fac=temp0/(temp0-temp1) ds=fac*(sr(i)-sr(i-1)) fmconv(nconv)=fms-(sr(i-1)+ds) iconv(nconv)=i endif temp0=temp1 10 continue c c When through, nconv is the number of transitions c return end c subroutine PULSOUT implicit double precision (a-h,o-z) parameter(jmax=1999) c c This short subroutine takes as input the model array H and writes c out an ascii onto unit 20 that can be read by the PULS pulsation c code c c The variables used are: r, Mr, Lr, T, d, P, eps, kappa, c Cv, Xd, Xt, del, delad, and Y c c Common block STRUC contains model structure common/STRUC/h(jmax,30),comp(jmax,15),sr(jmax),nj,modno c Common block GLOB contains global model parameters: xat0 is the initial X common/GLOB/fms,xat0,z,alfa,time,dtime,pcl,tcl,fls,teffl 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 dimension xpuls(jmax),rpuls(jmax) dimension gpuls(jmax),sl2(jmax),bv2(jmax),upuls(jmax),vpuls(jmax) write (20,*) fms write (20,*) comp(nj,1),comp(nj,4) write (20,*) nj do 10 j=1,nj rpuls(j)=dexp(h(j,1))*crs r2=rpuls(j)*rpuls(j) r3=r2*rpuls(j) gpuls(j)=cg*(fms-sr(j))*cms/r2 gamma1=h(j,11)/(1.d0-h(j,10)*h(j,6)) xpuls(j)=dlog(rpuls(j))-h(j,2) pscale=dexp(h(j,2)-h(j,5))/gpuls(j) upuls(j)=4.d0*cpi*r3*dexp(h(j,5))/((fms-sr(j))*cms) vpuls(j)=rpuls(j)/pscale vpuls(j)=1.d0/(1.d0+vpuls(j)) sl2(j)=gpuls(j)*pscale*gamma1/r2 bv2(j)=-h(j,10)*gpuls(j)*(h(j,8)-h(j,6))/h(j,11)/pscale 10 continue do 20 j=1,nj write (20,3000) xpuls(j),rpuls(j),gpuls(j),dexp(h(j,5)) 20 continue write (20,*) nj do 30 j=1,nj write (20,3000) xpuls(j),gpuls(j)/rpuls(j), 1 gpuls(j)/sl2(j)/rpuls(j), rpuls(j)*bv2(j)/gpuls(j) write (20,3001) upuls(j),vpuls(j) 30 continue 3000 format(1p4d14.6) 3001 format(1p2d14.6) return end c subroutine MODOUT implicit double precision (a-h,o-z) parameter(jmax=1999) c c This short subroutine takes as input the model array H and writes c out a model file onto unit 30 that can be read by the evolution code. c c Common block STRUC contains model structure common/STRUC/h(jmax,30),comp(jmax,15),sr(jmax),nj,modno c Common block GLOB contains global model parameters: xat0 is the initial X common/GLOB/fms,xat0,z,alfa,time,dtime,pcl,tcl,fls,teffl 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 Write out global parameters write (30,100) fms,xat0,z,alfa,time,dtime,pcl,tcl 1 ,fls,teffl,nj,modno 100 format(1p10e23.15,0p,i8,i8) c Write out array of surface masses write (30,'(1p3e23.15)') (sr(j),j=1,nj) c Write out short array of structural quantities write (30,'(1p3e23.15)') ((H(j,i),j=1,nj),i=1,5) write (30,'(1p3e23.15)') ((H(j,i),j=1,nj),i=19,20) c Write out array of compositional quantities write (30,'(1p3e23.15)') ((COMP(j,i),j=1,nj),i=1,15) c Done! return end c subroutine LOAD2(s2,v2,y2) c c central conditions c implicit double precision (a-h,o-z) parameter(jmax=1999) c c Load up starting vector for shooting at center c c Common block STRUC contains model structure common/STRUC/h(jmax,30),comp(jmax,15),sr(jmax),nj,modno c Common block GLOB contains global model parameters: xat0 is the initial X common/GLOB/fms,xat0,z,alfa,time,dtime,pcl,tcl,fls,teffl c Common block EQRATE contains info for equilibrium hydrogen burning common/EQRATE/ieqep,ieqecn c Common block CONST contains fundamental physical constants COMMON/CONST/CMS,CRS,CLS,CG,CSB,CPI,CP4,CCL,CME,CH,CNA,CEV,CK,CSYR common/TORHS/cmp(15),sp,st,ishell c dimension cmp(15),dxdt(15) dimension dxdt(15) dimension v2(2),y2(4) c pc=10.d0**v2(1) tc=10.d0**v2(2) x=comp(1,1) y=comp(1,4) xc12=comp(1,8) c Determine constitutive quantities at central point tcut=5.5d0 call EOS(Pc,tc,x,y,xc12,tcut,q,cp,ga,dc,xd,xt,u,cv,eta,fmu) c fmr=fms-s2 do 10 i=1,15 cmp(i)=comp(1,i) 10 continue c Establish equilibrium abundances for species participating Hydrogen c burning, and load these into abundance array at center. c call XMASS(dlog10(dc),dlog10(Tc),fmr,dtime,cmp) x=cmp(1)+cmp(2) y=cmp(3)+cmp(4) xc12=cmp(8) xo16=cmp(12) call opacity(dlog10(dc),dlog10(Tc),x,y,Z,XC12,XO16, 1 fk,0,DLKLT,DLKLD) call BNUKE(dlog10(dc),dlog10(tc),cmp,dxdt,eps,ieqep,epptot, 1 ieqecn,ecno,e3a,eac,enu,0,dlelt,dleld) c c Use central expansions to determine next-step-out quantities c c Radius c ccont=3.d0*cms/(cp4*crs**3) y2(1)=(1.d0/3.d0)*dlog(fmr*ccont/dc) c c Pressure c chse=cg/2.d0*(cp4/3.d0)**(1.d0/3.d0)*cms**(2.d0/3.d0) temp=pc-chse*dc**(4.d0/3.d0)*fmr**(2.d0/3.d0) y2(2)=dlog(temp) c c Luminosity c cfl=cms/cls c print *,'eps=',eps y2(4)=cfl*fmr*eps c c Temperature c c find the radiative gradient temp=3.D0*CLS/(16.D0*CP4*CG*CSB*CMS) GR=temp*y2(4)*fK*pc/fmr/tc**4 c check for convective stability IF(GR.LE.GA) THEN c Stable against convection; set gradient to radiative gradient GRAD=GR ELSE c Convectively unstable; use ADIABATIC temperature gradient for core c convection. To be fancy, could call TGRAD here, but nah. grad=ga endif temp=tc-chse*dc**(4.d0/3.d0)*tc/pc*grad*fmr**(2.d0/3.d0) y2(3)=dlog(temp) c c print '(a3,0p6f12.5)','y2s',dlog10(pc),dlog10(tc), c 1 (y2(i)/dlog(10.d0),i=1,3),y2(4) c print *,ieqep,ieqecn,dtime return end c subroutine LOAD1(s1,v1,y1) implicit double precision (a-h,o-z) parameter(jmax=1999,nmax=50,kmaxx=2000) c c Surface conditions: c Load up starting vector for shooting down from photosphere c c Common block STRUC contains model structure common/STRUC/h(jmax,30),comp(jmax,15),sr(jmax),nj,modno c Common block GLOB contains global model parameters: xat0 is the initial X common/GLOB/fms,xat0,z,alfa,time,dtime,pcl,tcl,fls,teffl c Common block CONST contains fundamental physical constants COMMON/CONST/CMS,CRS,CLS,CG,CSB,CPI,CP4,CCL,CME,CH,CNA,CEV,CK,CSYR common/TORHS/cmp(15),sp,st,ishell c Common block EQRATE contains info for equilibrium hydrogen burning common/EQRATE/ieqep,ieqecn dimension v1(2),y1(4) c common/path/kmax common/path/kmax,kount,dxsav,xp(KMAXX),yp(NMAX,KMAXX) external RHS,RKQS kmax=0 c fl=10.d0**v1(1) teff=10.d0**v1(2) c x=comp(nj,1)+comp(nj,2) y=comp(nj,3)+comp(nj,4) XC12=comp(nj,8) XO16=comp(nj,12) c c First the atmospheric integration to obtain the pressure and mass c at the photosphere. call ATMINT(fl,teff,alfa,x,y,z,XC12,XO16,fms,fmenv,Ro,peff) reff=dsqrt(cls*fl/cp4/csb/teff**4)/crs y1(1)=dlog(reff) y1(2)=dlog(peff) y1(3)=dlog(Teff) y1(4)=fl do 10 i=1,15 cmp(i)=comp(nj,i) 10 continue c Zero out the entropy terms sp=0.d0 st=0.d0 c c Now integrate downwards to interior mass point with Sr=S1 h1=fmenv/10.d0 call ODEINT(y1,4,fmenv,s1,1.d-2,h1,0.d0,nok,nbad,RHS,rkqs) return end c subroutine SCORE(xf,y,f) implicit double precision (a-h,o-z) dimension f(4),y(4) c c Dummy subroutine to map array y into array f, used in SHOOTF c do 10 i=1,4 f(i)=y(i) 10 continue print 100,xf,(y(i),i=1,4) 100 format (' Score:',f8.4,4f12.5) return end c subroutine FILHC(J,flnR,flnP,flnT,flr,srin) implicit double precision (a-h,o-z) parameter(jmax=1999) c c Subroutine to load the H and COMP arrays with the auxiliary quantities c for a model at zone j. c c Common block STRUC contains model structure common/STRUC/h(jmax,30),comp(jmax,15),sr(jmax),nj,modno c Common block GLOB contains global model parameters: xat0 is the initial X common/GLOB/fms,xat0,z,alfa,time,dtime,pcl,tcl,fls,teffl c Common block EQRATE contains info for equilibrium hydrogen burning common/EQRATE/ieqep,ieqecn c Common block CONST contains fundamental physical constants COMMON/CONST/CMS,CRS,CLS,CG,CSB,CPI,CP4,CCL,CME,CH,CNA,CEV,CK,CSYR dimension cmp(15),dcdt(15),cmpin(15) c fmr=fms-srin R=dexp(flnR) P=dexp(flnP) T=dexp(flnT) X=comp(j,1) XH2=comp(j,2) XHe3=comp(j,3) XHe4=comp(j,4) XC12=comp(j,8) XC13=comp(j,9) XN14=comp(j,10) XN15=comp(j,11) XO16=comp(j,12) H(j,1)=flnR H(j,2)=flnP H(j,3)=flnT H(j,4)=flr do 9 i=1,15 cmpin(i)=comp(j,i) cmp(i)=cmpin(i) 9 continue c re-establish X and Y as the sum of all isotopes X=cmp(1)+cmp(2) Y=cmp(3)+cmp(4) c c Call equation of state routine to obtain density and thermo derivatives c tcut=5.5d0 call EOS(P,T,x,y,XC12,tcut,q,cp,ga,d,xd,xt,u,cv,eta,fmu) c Find equilibrium abundances of X,XH2,XHe3,and XHe4 with PPEQ. c call XMASS(dlog10(d),dlog10(T),fmr,dtime,cmp) do 10 i=1,15 comp(j,i)=cmp(i) 10 continue c and load into the H array H(j,5)=dlog(d) H(j,6)=ga H(j,9)=q H(j,10)=xt H(j,11)=xd H(j,18)=fmu c c Call opacity routine to get opacity and its derivatives c call opacity(dlog10(d),dlog10(T),x,y,Z,XC12,XO16, 1 fk,1,DLKLT,DLKLD) c and load into H array H(j,12)=fk H(j,13)=dlklt H(j,14)=dlkld c c Call nuclear burning routine to get epsilon and its derivatives c call BNUKE(dlog10(d),dlog10(T),cmp,dcdt,eps,ieqep,epptot, 1 ieqecn,ecno,e3a,eac,enu,1,dlelt,dleld) c and load into H array H(j,15)=eps H(j,16)=dlelt H(j,17)=dleld c c find the radiative gradient c Crad=3.D0*CLS/(16.D0*CP4*CG*CSB*CMS) GR=Crad*flr*fK*P/fmr/T**4 H(j,7)=gr c And the true temperature gradient: check for convective stability IF(GR.LE.GA) THEN c stable against convection; set gradient to radiative gradient GRAD=GR ELSE c Convectively unstable; use mixing length theory to evaluate c the temperature gradient... if T is less than 1e7 K if (t.lt.1.d7) then call tgrad(fmr,R,T,P,d,alfa,infmlt,q,cp,fk,ga,gr,grad) else grad=ga endif endif H(j,8)=grad c c Zero out the entropy terms for now c H(j,19)=0.d0 H(j,20)=0.d0 c c Done! c return end c c Subroutine SHOOTF is now called funcv for use with NEWT c Here it is modified so that the subroutine RHS is used in place of c DERIVS SUBROUTINE funcv(n,v,f) IMPLICIT DOUBLE PRECISION (A-H,O-Z) c c Routine for use with newt to solve a two point boundary value problem for c nvar coupled ODEs by shooting from X1 and X2 to a fitting point XF c Initial values for the NVAR ode's at X1 (X2) are generated from the N2 (N1) c coefficients V1 (V2), using the user supplied routine LOAD1 (LOAD2). c The coefficients v1(v2) should be stored in a single array V(N1+N2) in the c main program by an EQUIVALENCE statement of the form c EQUIVALENCE (v1(1),v(1)),(v2(1),v(n2+1)) c The input parameter n=n1+n2=nvar. c The routine integrates the ode's to XF using the Runge-Kutta method with c tolerance EPS, initial step size H1, and minimum step size HMIN. At XF it c calls the user supplied subroutine SCORE to evaluate the NVAR functions F c that ought to match at XF. The differences F are returned on output. NEWT c uses a globally convergent Newton's method to adjust the values of V until c the functions F are zero. The user supplied subroutine DERIVS(x,y,dydx) c supplies derivative information to the ODE integrator. The common block c CALLER receives its values from the main program so that funcv can have c the syntax required by newt. Set nn2=n2 in the main program. The common c block PATH is for compatibility with ODEINT c Numerical Recipes, p. 753. c c Uses derivs, load1, load2, odeint, rkqs, score INTEGER n,nvar,nn2,kmax,kount,KMAXX,NMAX DOUBLE PRECISION f(n),v(n),x1,x2,xf,dxsav,xp,yp,EPS PARAMETER (NMAX=50,KMAXX=2000,EPS=1.d-8) COMMON /caller/ x1,x2,xf,nvar,nn2 COMMON /path/ kmax,kount,dxsav,xp(KMAXX),yp(NMAX,KMAXX) INTEGER i,nbad,nok DOUBLE PRECISION h1,hmin,f1(NMAX),f2(NMAX),y(NMAX) EXTERNAL rhs,rkqs,bsstep c kmax=0 c h1=(x2-x1)/100.d0 c h1=x1 h1=x1/2.d0 hmin=1.d-10 call load1(x1,v,y) c path from X1 to XF with best trial values V1. call odeint(y,nvar,x1,xf,EPS,h1,hmin,nok,nbad,rhs,rkqs) call score(xf,y,f1) call load2(x2,v(nn2+1),y) c path from X2 to XF with best trial values V2. h1=-(x2-xf)/100.d0 call odeint(y,nvar,x2,xf,EPS,h1,hmin,nok,nbad,rhs,rkqs) call score(xf,y,f2) do 11 i=1,n f(i)=f1(i)-f2(i) 11 continue return end c subroutine atmint(fl,Teff,alfa,x,y,z,c12,o16,fmtot,fmenv,Ro,peff) implicit double precision (a-h,o-z) c c Integrate an atmosphere with given L/Lo (fl) and Teff. c from tau=0 to tau=2/3 to obtain pressure at photosphere c and mass (in solar masses) above the photosphere, and the c radius (in solar radii) at the top of the atmosphere c c INPUT: fl = L/Lsun at surface c Teff = Effective temperature c alfa = mixing length to pressure scale height ratio c x = hydrogen mass fraction c y = helium4 mass fraction c c12 = C12 mass fraction c o16 = O16 mass fraction c fmtot = total mass of star, solar units c c OUTPUT: c fmenv = mass between tau=0 and photosphere, solar units c Ro = estimate of radius at tau=0, solar units c Peff = pressure at tau=2/3 c c Uses Numerical Recipes integrator, and a subroutine (DIFATM) c to provide r.h.s.'s to continuity, hse, radiative diffusion, c and optical depth equation. Allows for convective transport, too. c c Assumes that the atmosphere has a constant gravity with c R=Reff. The T(tau) relation is set by the equation of radiative c diffusion or convection, and turns out to be equivalent to the Eddington c approximation when radiative. c 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 Common block ATMIN contains things needed to evaluate atmosphere c structure equations in DIFATM, but not in calling arguments. c Implemented to preserve structure of Numerical Recipes routines. common/atmin/xc,yc,zc,cc,oc,flc,fmc,reff,grav,alfac c Array struc(4) contains r/Rsun, lnP, lnT, and S (=1-Mr/Msun) dimension struc(4),dydx(4) external deratm c c load common block atmin for transferring variable to r.h.s. subroutine xc=x yc=y zc=z cc=c12 oc=o16 flc=fl fmc=fmtot alfac=alfa c c Conditions at tau=0 (Eddington approximation) c c To=Te/2^1/4 c Po= radiation pressure at that T, times 1.1 c So=0.d0 by definition c Ro is about equal to Reff; it doesn't matter at this level of approximation c we keep track of it though to estimate the `true' Ro. c struc(4)=0.d0 To=Teff/(2.d0**0.25d0) struc(3)=dlog(To) Po=1.1d0*4.d0*csb/3.d0/ccl*To**4 struc(2)=dlog(po) reff=dsqrt(cls*fl/cp4/csb/Teff**4)/crs struc(1)=reff tau=1.d-10 dtau0=1.d-10 taup=2.d0/3.d0 grav=cg*fmtot*cms/reff**2/crs**2 c c Okay, now integrate downwards to the photosphere c iastep=0 do while (tau.lt.taup) call DERATM(tau,struc,dydx) c Integrate down a step, keeping all changes to less than 1.d-3 eps=1.d-4 call BSSTEP(struc,dydx,4,tau,dtau0,eps,struc,dtau, 1 dtaunxt,deratm) dtau0=dtaunxt c Check to see if we'll pass the photosphere at the next step. If so c shorten up on the step to ensure that we don't overshoot. if (tau+dtau0.gt.taup) dtau0=taup-tau istep=istep+1 enddo c c Evaluate the radius at the top of the atmosphere c Ro=Reff/struc(1)*Reff peff=dexp(struc(2)) fmenv=struc(4) return end c subroutine DERATM(tau,struc,dydx) implicit double precision (a-h,o-z) dimension struc(4),dydx(4) c Common block CTRL contains most of the control parameters for ISUEVO common/CTRL/dllenv,dltenv,tcut,dtfac,sstep(7),tstep(7), 1 nmod,nsav,nprint,nzp,itstep c Common block ATMIN contains things needed to evaluate atmosphere c structure equations but not in calling arguments. Implemented to c preserve structure of Numerical Recipes routines. common/ATMIN/x,y,z,c12,o16,fl,fmtot,reff,grav,alfa 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 Common block MLT contains switch for ML parameters COMMON/MLT/infmlt c c Differential equations used to integrate c the atmosphere down from tau=0 to tau=2/3 to find the c photospheric pressure. For use in the Numerical Recipes routines c such as RK4, RKQC, or BSSTEP c r=struc(1)*crs p=dexp(struc(2)) t=dexp(struc(3)) tcut=5.5d0 call EOS(p,t,x,y,c12,tcut,q,cp,ga,d,xd,xt,u,cv,eta,fmu) call OPACITY(dlog10(d),dlog10(t),x,y,Z,c12,o16,fk,0 1 ,DLKLT,DLKLD) c Definition of optical depth gives dr/dtau dydx(1)=-1.d0/fk/d/crs c H.S.E. gives d(lnP)/dtau dydx(2)=grav/fk/p c Radiation diffusion or convection gives d(lnT)/dtau condl=3.d0*cls*fl/(16.d0*cp4*crs**2*csb) c Radiative temperature gradient: gr=condl*fK*P/(grav*reff**2*T**4) c check for convective stability IF(gr.LE.GA) THEN c stable against convection; set gradient to radiative gradient GRAD=GR ELSE c Convectively unstable; use mixing length theory to evaluate c the temperature gradient... call tgrad(fmtot,reff,t,p,d,alfa,infmlt,q,cp,fk,ga,gr,grad) endif c dydx(3)=grav*grad/p/fk dydx(4)=cp4*Reff**2/fk*crs**2/cms c print *,t,p,dydx(1),dydx(2),dydx(3),dydx(4),fk return end c subroutine bsstep(y,dydx,nv,x,htry,eps,yscal,hdid,hnext,DERIVS) implicit double precision (a-h,o-z) c c Bulirsch-Stoer step with monitoring of local truncation error to c ensure accuracy and adjust stepsize. c INPUTS: c Y(NV) dependant variable vector of length NV c X value of independent variable c DYDX(NV) deriv of Y w.r.t. X at the starting value of X c HTRY first try at the step size in X c EPS required accuracy c YSCAL(NV) vector that scales EPS c c on OUTPUT: c Y(NV) values at end of step c X value of X at end of step c HDID size of step in X c HNEXT estimated next step size c c subroutine DERIVS supplies the derivatives of Y(NV) with respect to X c PARAMETER (NMAX=50,KMAXX=8,IMAX=KMAXX+1,SAFE1=0.25d0,SAFE2=0.7d0, 1 REDMAX=1.d-5,REDMIN=0.7d0,TINY=1.d-30,SCALMX=0.1d0) DIMENSION Y(NV),DYDX(NV),YSCAL(NV),YERR(NMAX), 1 a(IMAX),alf(KMAXX,KMAXX),err(KMAXX), 2 ysav(NMAX),yseq(NMAX),nseq(IMAX) LOGICAL first,reduct SAVE a,alf,epsold,first,kmax,kopt,nseq,xnew EXTERNAL derivs DATA first/.true./,epsold/-1.d0/ DATA NSEQ/2,4,6,8,10,12,14,16,18/ c a new tolerance, so reinitialize if (eps.ne.epsold) then hnext=-1.d29 xnew=-1.d29 eps1=SAFE1*eps c compute work coefficients A(k) a(1)=nseq(1)+1 do 11 k=1,KMAXX a(k+1)=a(k)+nseq(k+1) 11 continue c compute alfa(k,q) do 13 iq=2,KMAXX do 12 k=1,iq-1 alf(k,iq)=eps1**((a(k+1)-a(iq+1))/ 1 ((a(iq+1)-a(1)+1.d0)*(2*k+1))) 12 continue 13 continue epsold=eps c Determine the optimal row number for convergence do 14 kopt=2,KMAXX-1 if (a(kopt+1).gt.a(kopt)*alf(kopt-1,kopt)) goto 1 14 continue 1 kmax=kopt endif H=HTRY C save the starting value of x,y(nv) DO 15 I=1,NV YSAV(I)=Y(I) 15 CONTINUE c A new stepsize or a new integration; re-establish the order window if (h.ne.hnext.or.x.ne.xnew) then first=.true. kopt=kmax endif reduct=.false. C EVALUATE THE SEQUENCE OF MODIFIED MIDPOINT INTEGRATIONS 2 DO 17 k=1,kmax xnew=x+h if (xnew.eq.x) STOP 'step size underflow in BSSTEP' CALL MMID(ysav,dydx,nv,x,h,nseq(k),yseq,DERIVS) C SQUARED BELOW SINCE ERROR SERIES IS EVEN XEST=(H/NSEQ(k))**2 C PERFORM EXTRAPOLATION CALL RZEXTR(k,XEST,YSEQ,Y,YERR,NV) c CALL PZEXTR(k,XEST,YSEQ,Y,YERR,NV) C Compute normalized error estimate eps(k) if (k.ne.1) then ERRMAX=tiny DO 16 i=1,NV ERRMAX=DMAX1(ERRMAX,DABS(YERR(i)/YSCAL(i))) 16 CONTINUE C SCALE error RELATIVE TO TOLERANCE ERRMAX=ERRMAX/EPS km=k-1 err(km)=(errmax/SAFE1)**(1.d0/(2*km+1)) endif if (k.ne.1.and.(k.ge.kopt-1.or.first))then c converged IF (ERRMAX.LT.1.d0) goto 4 C check for possible stepsize reduction if (k.eq.kmax.or.k.eq.kopt+1) then red=SAFE2/err(km) goto 3 else if (k.eq.kopt) then if (alf(kopt-1,kopt).lt.err(km))then red=1.d0/err(km) goto 3 endif else if (kopt.eq.kmax) then if (alf(km,kmax-1).lt.err(km))then red=alf(km,kmax-1)*SAFE2/err(km) goto 3 endif else if (alf(km,kopt).lt.err(km))then red=alf(km,kopt-1)/err(km) goto 3 endif endif 17 continue c Reduce stepsize by at least REDMIN and at most REDMAX 3 red=dmin1(red,REDMIN) red=dmax1(red,REDMAX) h=h*red reduct=.true. c try again goto 2 c Succesful step taken! 4 x=xnew hdid=h first=.false. wrkmin=1.d35 c compute optimal row for convergence and corresponding stepsize do 18 kk=1,km fact=dmax1(err(kk),SCALMX) work=fact*a(kk+1) if (work.lt.wrkmin) then scale=fact wrkmin=work kopt=kk+1 endif 18 continue hnext=h/scale c Check for possible order increase but not if stepsize was just reduced if (kopt.ge.k.and.kopt.ne.kmax.and..not.reduct) then fact=dmax1(scale/alf(kopt-1,kopt),SCALMX) if (a(kopt+1)*fact.le.wrkmin)then hnext=h/fact kopt=kopt+1 endif endif return end C SUBROUTINE MMID(Y,DYDX,NVAR,XS,HTOT,NSTEP,YOUT,DERIVS) IMPLICIT DOUBLE PRECISION (A-H,O-Z) C c Modified midpoint step. Dependent variable vector Y(NVAR) and c its derivative vector DYDX(NVAR) are input at XS. Also input is c HTOT, the total step to be made, and NSTEP, the number of substeps c to be used. The output is returned as YOUT(NVAR), which need not be c a distinct array from Y(NVAR); if it is distinct, however, then Y(NVAR) c and DYDX(NVAR) are returned undamaged c PARAMETER(NMAX=50) external DERIVS DIMENSION Y(NVAR),DYDX(NVAR),YOUT(NVAR),YM(NMAX),YN(NMAX) C STEPSIZE THIS TRIP H=HTOT/NSTEP DO 11 I=1,NVAR YM(I)=Y(I) C FIRST STEP YN(I)=Y(I)+H*DYDX(I) 11 CONTINUE X=XS+H C WILL USE YOUT FOR TEMPORARY STORAGE OF DERIVATIVES CALL DERIVS(X,YN,YOUT) H2=2.D0*H C GENERAL STEP DO 13 N=2,NSTEP DO 12 I=1,NVAR SWAP=YM(I)+H2*YOUT(I) YM(I)=YN(I) YN(I)=SWAP 12 CONTINUE X=X+H CALL DERIVS(X,YN,YOUT) 13 CONTINUE C LAST STEP DO 14 I=1,NVAR YOUT(I)=0.5d0*(YM(I)+YN(I)+H*YOUT(I)) 14 CONTINUE RETURN END c subroutine pzextr(iest,xest,yest,yz,dy,nv) implicit double precision (a-h,o-z) c c Uses polynomial extrapolation to evaluate nv functions at x=0 by fitting a c polynomial to a sequence of estimates with progressively smaller values c x=xest, and corresponding function vectors yest(nv). This call is number c iest in the sequence of calls. Extrapolated function values are output c as yz(nv), and their estimated error is output as dy(nv). c parameter(IMAX=13,NMAX=50) dimension dy(nv),yest(nv),yz(nv),d(NMAX),qcol(NMAX,IMAX), 1 x(IMAX) SAVE qcol,x c Save current independent variable x(iest)=xest do 11 j=1,nv dy(j)=yest(j) yz(j)=yest(j) 11 continue c Store first estimate in first column if (iest.eq.1) then do 12 j=1,nv qcol(j,1)=yest(j) 12 continue else do 13 j=1,nv d(j)=yest(j) 13 continue do 15 k1=1,iest-1 delta=1.d0/(x(iest-k1)-xest) f1=xest*delta f2=x(iest-k1)*delta c Propagate tableau 1 diagonal more do 14 j=1,nv q=qcol(j,k1) qcol(j,k1)=dy(j) delta=d(j)-q dy(j)=f1*delta d(j)=f2*delta yz(j)=yz(j)+dy(j) 14 continue 15 continue do 16 j=1,nv qcol(j,iest)=dy(j) 16 continue endif return END c SUBROUTINE RZEXTR(IEST,XEST,YEST,YZ,DY,NV) IMPLICIT DOUBLE PRECISION (A-H,O-Z) C C Uses diagonal rational function extrapolation to evaluate NV functions c at X=0 by fitting a diagonal rational functino to a sequence of estimates c with progressively smaller values X=XEST, and corresponding function c vectors YEST. This call is number IEST in the sequence of calls. c Extrapolated function values are output as YZ, and their estimated error is c output as DY. PARAMETER (IMAX=13,NMAX=50) C MAXIMUM EXPECTED VALUE OF NUSE IS NCOL, OF NV IS NMAX, OF IEST IS IMAX DIMENSION dy(nv),yest(nv),yz(nv),d(NMAX,IMAX),fx(IMAX),x(IMAX) SAVE d,x C SAVE CURRENT INDEPENDENT VARIABLE X(IEST)=XEST IF (IEST.EQ.1) THEN DO 11 J=1,NV YZ(J)=YEST(J) D(J,1)=YEST(J) DY(J)=YEST(J) 11 CONTINUE ELSE do 12 k=1,iest-1 fx(k+1)=x(iest-k)/xest 12 continue C EVALUATE NEXT DIAGONAL IN TABLEAU DO 14 J=1,NV YY=YEST(J) V=D(J,1) C=YY D(J,1)=YY DO 13 K=2,iest B1=FX(K)*V B=B1-C IF (B.NE.0.D0) THEN B=(C-V)/B DDY=C*B C=B1*B C CARE NEEDED TO AVOID DIVISION BY ZERO... ELSE DDY=V ENDIF if (k.ne.iest) v=d(j,k) D(J,K)=DDY YY=YY+DDY 13 CONTINUE DY(J)=DDY YZ(J)=YY 14 CONTINUE ENDIF RETURN END c subroutine difeq(k,k1,k2,jsf,is1,isf,indexv, 1 ne,s,nsi,nsj,y,nyj,nyk) implicit double precision (a-h,o-z) c c Subroutine to supply coefficients for SOLVDE for the c differential equations of stellar structure equations c k is the point at which the difference equations c are to be evaluated c k1 and k2 label the first and last points c in the mesh. c is1 : starting row (equation) that needs to be filled c isf : final row (equation) that needs to be filled c jsf : column containing the r.h.s. of the difference c equations; c y : y(j,k)value of variable j at meshpoint k c s : s(i,j), where i labels the equation, and j is the c index of the variable. j=1,n refers to the values c at meshpoint k-1, j=n+1,2n refers to the values c at meshpoint k. c c Note: variables postscripted with a 0 are for zone k-1 c variables postscripted with a 1 are for zone k c c Arrays for difference equation solution: parameter(jmax=1999) dimension y(nyj,nyk) dimension indexv(nyj),s(nsi,nsj) c Common block CTRL contains most of the control parameters for passing common/CTRL/dllenv,dltenv,tcut,dtfac,sstep(7),tstep(7), 1 nmod,nsav,nprint,nzp,itstep 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 Common block STRUC contains model structure common/STRUC/h(jmax,30),comp(jmax,15),sr(jmax),nj,modno c Common block GLOB contains global model parameters: xat0 is the initial X common/GLOB/fms,xat0,z,alfa,time,dtime,pcl,tcl,fls,teffl c Common block BETA contains offsets for difference equations common/BETA/br,bp,bt,bl c Common block TORHS contains things for use in computing DYDX c things exported are compostion variables cmp, sp and st are changes c in P and T from last model to this one common/TORHS/cmp(15),sp,st,ishell c Common block EQRATE contains info for equilibrium hydrogen burning common/EQRATE/ieqep,ieqecn c dimension yy(4),dydx0(4),dydx1(4),dcdt(15) c c h(j,1)=lnR c h(j,2)=lnP c h(j,3)=lnT c h(j,4)=Fl c h(j,5)=ln(rho) c h(j,6)=delad c h(j,7)=delrad c h(j,8)=del c h(j,9)=q =-dln(rho)/dln(T)|P=Xt/Xd c h(j,10)=Xt =dln(P)/dln(T)|rho c h(j,11)=Xd =dln(P)/dln(rho)|T c h(j,12)=kappa c h(j,13)=dln(kappa)/dln(T)|rho c h(j,14)=dln(kappa)/dln(rho)|T c h(j,15)=epsilon c h(j,16)=dln(epsilon)/dln(T)|rho c h(j,17)=dln(epsilon)/dln(rho)|T c h(j,18)=mu c h(j,19)=delta(lnP) from previous model to this model c h(j,20)=delta(lnT) " " " " " " c c CENTRAL BOUNDARY CONDITIONS: really trivial, as most of the work done c for central convergence is at zone 1 c if (k.eq.k1) then call DIFCBC(nyj,nyk,y,nsi,nsj,jsf,indexv,s) c print '(i5,1p4e13.5)',1,0.d0,0.d0,s(3,jsf),s(4,jsf) return endif c c SURFACE BOUNDARY CONDITIONS c if (k.gt.k2) then call DIFOBC(nyj,nyk,y,nsi,nsj,jsf,indexv,s) c print '(i5,1p4e13.5)',k,s(1,jsf),s(2,jsf),0.d0,0.d0 return endif c c EQUATIONS FOR THE INTERIOR POINTS c sr1=sr(k) sr0=sr(k-1) c fmr is Mr fmr1=fms-sr1 fmr0=fms-sr0 r1=dexp(y(1,k)) r0=dexp(y(1,k-1)) p1=dexp(y(2,k)) p0=dexp(y(2,k-1)) t1=dexp(y(3,k)) t0=dexp(y(3,k-1)) fl1=y(4,k) fl0=y(4,k-1) c careful here, if fl1 or fl0 is exactly zero then things blow up... if zero c then set them to something very small! if (fl1.eq.0.d0) fl1=1.d-99 if (fl0.eq.0.d0) fl0=1.d-99 cc Note that the H, COMP, SR arrays include the center cc as Zone 1, but the Y array for numerics doesn't. cc Hence the value of k is one greater when dealing cc with H, COMP, SR than with Y sp1=h(k,19) sp0=h(k-1,19) st1=h(k,20) st0=h(k-1,20) x1=comp(k,1)+comp(k,2) x0=comp(k-1,1)+comp(k-1,2) y1=comp(k,4)+comp(k,3) y0=comp(k-1,4)+comp(k-1,3) XC121=comp(k,8) XC120=comp(k-1,8) XO161=comp(k,12) XO160=comp(k-1,12) c c Find right hand sides of the differential equations c c load common block TORHS at zone k, then call RHS do 10 i=1,15 cmp(i)=comp(k,i) 10 continue YY(1)=y(1,k) YY(2)=y(2,k) YY(3)=y(3,k) YY(4)=y(4,k) sp=sp1 st=st1 c set shell number for RHS call so that you don't need to c call EOS within RHS ishell=k call RHS(sr1,yy,dydx1) c constant for figuring radiative gradient cgr=3.d0*cls/(16.d0*cp4*cg*csb*cms) c Now extract some useful things for zone K (subscript 1) q1=h(k,9) ga1=h(k,6) gr1=h(k,7) g1=h(k,8) d1=dexp(h(k,5)) xd1=h(k,11) xt1=h(k,10) fk1=h(k,12) dlklt1=h(k,13) dlkld1=h(k,14) eps1=h(k,15) dlelt1=h(k,16) dleld1=h(k,17) cp1=p1*q1/d1/t1/ga1 c gr1=cgr*fl1*fk1*p1/(fmr1*t1**4) c find the partials of kappa w.r.t. P at constant T, and T at constant P dlklp1=dlkld1/xd1 dlklt1=dlklt1-q1*dlkld1 c c load common block TORHS at zone k-1, then call RHS c do 11 i=1,15 cmp(i)=comp(k-1,i) 11 continue YY(1)=y(1,k-1) YY(2)=y(2,k-1) YY(3)=y(3,k-1) YY(4)=y(4,k-1) sp=sp0 st=st0 c set shell number for RHS call so that you don't need to c call EOS within RHS c ishell=k ishell=k-1 call RHS(sr0,yy,dydx0) c Now extract some useful things for zone K-1 (subscript 0) q0=h(k-1,9) ga0=h(k-1,6) gr0=h(k-1,7) g0=h(k-1,8) d0=dexp(h(k-1,5)) xd0=h(k-1,11) xt0=h(k-1,10) fk0=h(k-1,12) dlklt0=h(k-1,13) dlkld0=h(k-1,14) eps0=h(k-1,15) dlelt0=h(k-1,16) dleld0=h(k-1,17) cp0=p0*q0/d0/t0/ga0 c gr0=cgr*fl0*fk0*p0/(fmr0*t0**4) c find the partials of kappa w.r.t. P at constant T, and T at constant P dlklp0=dlkld0/xd0 dlklt0=dlklt0-q0*dlkld0 c c Pressure and temperature derivatives of Xt and Ga at zone k c also derivatives of gradient when convective... c p1del=dexp(dlog(P1)+0.001d0) c tcut=5.5d0 call EOS(p1del,t1,x1,y1,XC121,tcut,qd,dum,gad,dd,xdd,xtd,ud, 1 cvd,etad,fmud) dlgalp1=(gad-ga1)/0.001d0/ga1 dlqlp1=(qd-q1)/0.001d0/q1 c if convective, compute d(ln del)/d(ln P) if (ga1.lt.gr1) then cpd=p1del*qd/dd/t1/gad fkd=dexp(dlog(fk1)+0.001d0*dlklp1) grd=cgr*fl1*fkd*p1del/(fmr1*t1**4) if (gad.lt.grd) then call TGRAD(fmr1,r1,t1,p1del,dd,alfa,0,qd,cpd,fkd,gad,grd,gd) else gd=grd endif dlglp1=(gd-g1)/0.001d0/g1 endif c t1del=dexp(dlog(T1)+0.001d0) call EOS(p1,t1del,x1,y1,XC121,tcut,qd,dum,gad,dd,xdd,xtd,ud, 1 cvd,etad,fmud) dlgalt1=(gad-ga1)/0.001d0/ga1 dlqlt1=(qd-q1)/0.001d0/q1 c if convective, compute d(ln del)/d(ln T) if (ga1.lt.gr1) then cpd=p1*qd/dd/t1del/gad fkd=dexp(dlog(fk1)+0.001d0*dlklt1) grd=cgr*fl1*fkd*p1/(fmr1*t1del**4) if (gad.lt.grd) then call TGRAD(fmr1,r1,t1del,p1,dd,alfa,0,qd,cpd,fkd,gad,grd,gd) else gd=grd endif dlglt1=(gd-g1)/0.001d0/g1 c also fix up dlgfl1 (would be 0 if adiabatic, 1 if radiative) fac=(g1-ga1)/(gr1-ga1) dlgfl1=fac c write (31,101) k, dlglp1,1.d0+dlklp1,dlglt1,-4+dlklt1,dlgfl1 c write (31,102) k, dlgalp1, dlgalt1 c write (31,101) k, ga1,gr1,g1,del1 101 format(i4,1p2e11.3,3x,2e11.3,3x,2e11.3) 102 format(i4,11x,1pe11.3,14x,e11.3) else c not convective, use derivitaves of delrad for derivatives of del dlglp1=1.d0+dlklp1 dlglt1=-4.d0+dlklt1 dlgfl1=1.d0 endif c c and at zone k-1 p0del=dexp(dlog(P0)+0.001d0) call EOS(p0del,t0,x0,y0,XC120,tcut,qd,dum,gad,dd,xdd,xtd,ud, 1 cvd,etad,fmud) dlgalp0=(gad-ga0)/0.001d0/ga0 dlqlp0=(qd-q0)/0.001d0/q0 c if convective, compute d(ln del)/d(ln P) at zone k-1 if (ga0.lt.gr0) then c use a fresh value for the gradient cpd=p0del*qd/dd/t0/gad fkd=dexp(dlog(fk0)+0.001d0*dlklp0) grd=cgr*fl0*fkd*p0del/(fmr0*t0**4) if (gad.lt.grd) then call TGRAD(fmr0,r0,t0,p0del,dd,alfa,0,qd,cpd,fkd,gad,grd,gd) else gd=grd endif dlglp0=(gd-g0)/0.001d0/g0 endif c t0del=dexp(dlog(T0)+0.001d0) call EOS(p0,t0del,x0,y0,XC120,tcut,qd,dum,gad,dd,xdd,xtd,ud, 1 cvd,etad,fmud) dlgalt0=(gad-ga0)/0.001d0/ga0 dlqlt0=(qd-q0)/0.001d0/q0 c if convective, compute d(ln del)/d(ln T) at zone k-1 if (ga0.lt.gr0) then cpd=p0*qd/dd/t0del/gad fkd=dexp(dlog(fk0)+0.001d0*dlklt0) grd=cgr*fl0*fkd*p0/(fmr0*t0del**4) if (gad.lt.grd) then call TGRAD(fmr0,r0,t0del,p0,dd,alfa,0,qd,cpd,fkd,gad,grd,gd) else gd=grd endif dlglt0=(gd-g0)/0.001d0/g0 c also fix up dlgfl0 (would be 0 if adiabatic, 1 if radiative) fac=(g0-ga0)/(gr0-ga0) dlgfl0=fac c write (31,101) k-1, dlglp0,1.d0+dlklp0,dlglt0,-4+dlklt0,dlgfl0 c write (31,102) k-1, dlgalp0, dlgalt0 c write (31,101) k-1, ga0,gr0,g0,del0 c write (31,101) else c not convective, use derivitaves of delrad for derivatives of del dlglp0=1.d0+dlklp0 dlglt0=-4.d0+dlklt0 dlgfl0=1.d0 endif c ================================================================== c And now for the difference equations... c the step in surface mass fraction: c ds=sr(k+1)-sr(k) ds=sr(k)-sr(k-1) c c Continuity equation (equation for lnR): c drnds1=dydx1(1) drnds0=dydx0(1) s(1,4+indexv(1))=1.d0+3.d0*ds*br*drnds1 s(1,indexv(1))=-1.d0+3.d0*ds*(1.d0-br)*drnds0 s(1,4+indexv(2))=ds*br*drnds1/Xd1 s(1,indexv(2))=ds*(1.d0-br)*drnds0/Xd0 s(1,4+indexv(3))=-ds*br*drnds1*Q1 s(1,indexv(3))=-ds*(1.d0-br)*drnds0*Q0 s(1,4+indexv(4))=0.d0 s(1,indexv(4))=0.d0 s(1,jsf)=dlog(r1)-dlog(r0)-ds*(br*drnds1+(1.d0-br)*drnds0) c c Hydrostatic equilibrium (equation for lnP) c dpnds1=dydx1(2) dpnds0=dydx0(2) s(2,4+indexv(1))=4.d0*ds*bp*dpnds1 s(2,indexv(1))=4.d0*ds*(1.d0-bp)*dpnds0 s(2,4+indexv(2))=1.d0+ds*bp*dpnds1 s(2,indexv(2))=-1.d0+ds*(1.d0-bp)*dpnds0 s(2,4+indexv(3))=0.d0 s(2,indexv(3))=0.d0 s(2,4+indexv(4))=0.d0 s(2,indexv(4))=0.d0 s(2,jsf)=dlog(p1)-dlog(p0)-ds*(bp*dpnds1+(1.d0-bp)*dpnds0) c c Thermal equilibrium (equation for lnT) c c gradients of gradients computed earlier whether convective or radiative... c dTnds1=dydx1(3) dTnds0=dydx0(3) s(3,4+indexv(1))=4.d0*ds*bt*dTnds1 s(3,indexv(1))=4.d0*ds*(1.d0-bt)*dTnds0 s(3,4+indexv(2))=ds*bt*dTnds1*(1.d0-dlglp1) s(3,indexv(2))=ds*(1.d0-bt)*dTnds0*(1.d0-dlglp0) s(3,4+indexv(3))=1.d0-ds*bt*dTnds1*dlglt1 s(3,indexv(3))=-1.d0-ds*(1.d0-bt)*dTnds0*dlglt0 s(3,4+indexv(4))=-ds*bt*dTnds1/fl1*dlgfl1 s(3,indexv(4))=-ds*(1.d0-bt)*dTnds0/fl0*dlgfl0 s(3,jsf)=dlog(T1)-dlog(T0)-ds*(bt*dTnds1+(1.d0-bt)*dTnds0) c c Energy conservation c c first some preliminaries, a-la Prather's thesis, for zone K. cl=-cms/cls sbar1=P1*q1/d1/ga1*(st1-ga1*sp1) dFlds1=dydx1(4) c Note that the derivatives of Sbar are flaky if you assume that c d/dlnP(dlnP/dt) is 1/dtime... the older (Prather) derivatives are c commented out c dsbTn1=sbar1*(q1+dlqlt1)+q1*p1/d1/ga1* c 1 (-dlgalt1*st1) c dsbPn1=sbar1*(1.d0+dlqlp1-1.d0/Xd1)- c 1 q1*p1/d1/ga1*(st1*dlgalp1) dsbTn1=sbar1*(q1+dlqlt1)+q1*p1/d1/ga1* 1 (1.d0-dlgalt1*st1) dsbPn1=sbar1*(1.d0+dlqlp1-1.d0/Xd1)- 1 q1*p1/d1/ga1*(ga1+st1*dlgalp1) c and some preliminaries, a-la Prather's thesis, for zone K-1. sbar0=P0*q0/d0/ga0*(st0-ga0*sp0) dFlds0=dydx0(4) c Note that the derivatives of Sbar are flaky if you assume that c d/dlnP(dlnP/dt) is 1/dtime... the older (Prather) derivatives are c commented out c dsbTn0=sbar0*(q0+dlqlt0)+q0*p0/d0/ga0* c 1 (-dlgalt0*st0) c dsbPn0=sbar0*(1.d0+dlqlp0-1.d0/Xd0)- c 1 q0*p0/d0/ga0*(st0*dlgalp0) dsbTn0=sbar0*(q0+dlqlt0)+q0*p0/d0/ga0* 1 (1.d0-dlgalt0*st0) dsbPn0=sbar0*(1.d0+dlqlp0-1.d0/Xd0)- 1 q0*p0/d0/ga0*(ga0+st0*dlgalp0) c s(4,4+indexv(1))=0.d0 s(4,indexv(1))=0.d0 s(4,4+indexv(2))=-cl*ds*bl*(eps1*(dleld1/Xd1)-dsbPn1/dtime) s(4,indexv(2))=-cl*ds*(1.d0-bl)* 1 (eps0*(dleld0/Xd0)-dsbPn0/dtime) s(4,4+indexv(3))=-cl*ds*bl*(eps1*(dlelt1-dleld1*q1) 1 -dsbTn1/dtime) s(4,indexv(3))=-cl*ds*(1.d0-bl)*(eps0*(dlelt0-dleld0*q0) 1 -dsbTn0/dtime) s(4,4+indexv(4))=1.d0 s(4,indexv(4))=-1.d0 s(4,jsf)=fl1-fl0-ds*(bl*dFlds1+(1.d0-bl)*dFlds0) c c Done with equations for the interior! c c print 100,k,(s(i,jsf),i=1,4) 100 format(i5,1p4e13.5) return end c subroutine DIFCBC(nyj,nyk,y,nsi,nsj,jsf,indexv,s) implicit double precision (a-h,o-z) parameter(jmax=1999) dimension y(nyj,nyk),indexv(nyj),s(nsi,nsj) c Common block CTRL contains most of the control parameters for passing common/CTRL/dllenv,dltenv,tcut,dtfac,sstep(7),tstep(7), 1 nmod,nsav,nprint,nzp,itstep c Common block STRUC contains model structure common/STRUC/h(jmax,30),comp(jmax,15),sr(jmax),nj,modno c Common block GLOB contains global model parameters: xat0 is the initial X common/GLOB/fms,xat0,z,alfa,time,dtime,pcl,tcl,fls,teffl 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 Common block EQRATE contains info for equilibrium hydrogen burning common/EQRATE/ieqep,ieqecn dimension cmp(15),dcdt(15) c c Central Boundary conditions c c some preliminaries: note that CBC's applied at zone 2. Thermo c quantities derived from pcl,tcl c fac=1.d0/dlog(10.d0) fmr=fms-sr(1) r=dexp(y(1,1)) fl=y(4,1) p=dexp(y(2,1)) t=dexp(y(3,1)) pc=dexp(pcl/fac) tc=dexp(tcl/fac) c p=pc t=tc c sp=h(1,19) st=h(1,20) x=comp(1,1)+comp(1,2) YY=comp(1,3)+comp(1,4) XC12=comp(1,8) XO16=comp(1,12) do 10 i=1,15 cmp(i)=comp(1,i) 10 continue c tcut=5.5d0 c tcut=6.d0 call EOS(P,T,x,YY,XC12,tcut,q,cp,ga,d,xd,xt,u,cv,eta,fmu) call BNUKE(dlog10(d),dlog10(T),cmp,dcdt,eps,ieqep,epptot, 1 ieqecn,ecno,e3a,eac,enu,1,dlelt,dleld) c Derivatives of del-ad and q w.r.t. p and t Pdel=dexp(dlog(p)+0.001d0) call EOS(Pdel,T,x,YY,XC12,tcut,q1,cp1,ga1,d1,xd1,xt1,u1, 1 cv1,eta1,fmu1) dlgalp=(ga1-ga)/0.001d0/ga dlqlp=(q1-q)/0.001d0/q Tdel=dexp(dlog(t)+0.001d0) call EOS(P,Tdel,x,YY,XC12,tcut,q1,cp1,ga1,d1,xd1,xt1,u1, 1 cv1,eta1,fmu1) dlgalt=(ga1-ga)/0.001d0/ga dlqlt=(q1-q)/0.001d0/q c c Continuity boundary condition: Relates R with Rho(P,T) through Mr c ccont=3.d0/cp4*cms/crs**3 s(3,4+indexv(1))=1.d0 s(3,indexv(1))=0.d0 s(3,4+indexv(2))=1.d0/3.d0/xd s(3,indexv(2))=0.d0 s(3,4+indexv(3))=-1.d0/3.d0*q s(3,indexv(3))=0.d0 s(3,4+indexv(4))=0.d0 s(3,indexv(4))=0.d0 s(3,jsf)=dlog(r)-1.d0/3.d0*(dlog(ccont)+dlog(fmr)-dlog(d)) c c Luminosity boundary condition: c clum=cms/cls dtnds=0.d0 dpnds=0.d0 c sbar=P*q/d/ga*(st-ga*sp) sbar=0.d0 c Note that the derivatives of Sbar are flaky if you assume that c d/dlnP(dlnP/dt) is 1/dtime... the older (Prather) derivatives are c commented out c dsblt=sbar*(q+dlqlt)+q*P/d/ga*(-dlgalt*st) c dsblp=sbar*(1.d0+dlqlp-1.d0/Xd)-q*P/d*(st/ga*dlgalp) c* dsblt=sbar*(q+dlqlt)+q*P/d/ga*(1.d0-dlgalt*st) c* dsblp=sbar*(1.d0+dlqlp-1.d0/Xd)-q*P/d*(1.d0+st/ga*dlgalp) dsblt=0.d0 dsblp=0.d0 s(4,4+indexv(1))=0.d0 s(4,indexv(1))=0.d0 s(4,4+indexv(2))=-clum*fmr*(eps*(dleld/xd)-dsblp/dtime) s(4,indexv(2))=0.d0 s(4,4+indexv(3))=-clum*fmr*(eps*(dlelt-dleld*q)-dsblt/dtime) s(4,indexv(3))=0.d0 s(4,4+indexv(4))=1.d0 s(4,indexv(4))=0.d0 s(4,jsf)=fl-clum*fmr*(eps-sbar/dtime) return end c subroutine DIFOBC(nyj,nyk,y,nsi,nsj,jsf,indexv,s) implicit double precision (a-h,o-z) parameter(jmax=1999) dimension indexv(nyj),s(nsi,nsj),cmp(15),y(nyj,nyk) c Common block STRUC contains model structure common/STRUC/h(jmax,30),comp(jmax,15),sr(jmax),nj,modno c Common block GLOB contains global model parameters: xat0 is the initial X common/GLOB/fms,xat0,z,alfa,time,dtime,pcl,tcl,fls,teffl 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 Common block KIPP contains positions for the corners of the OBC triangle COMMON/KIPP/g1,g2,g3,cl(3),ct(3),brl(3),bpl(3),btl(3),krot c c Surface boundary conditions for relaxation scheme c c rn=y(1,nj-1) c pn=y(2,nj-1) c tn=y(3,nj-1) c fl=y(4,nj-1) rn=y(1,nj) pn=y(2,nj) tn=y(3,nj) fl=y(4,nj) fac=1.d0/dlog(10.d0) c call OBCS(dlog10(fl),teffl,sr(nj),rl,pl,tl, c 1 dlRlL,dlRlTe,dlPlL,dlPlTe,dlTlL,dlTlTe) call OBCS(dlog10(fls),teffl,sr(nj),rl,pl,tl, 1 dlRlL,dlRlTe,dlPlL,dlPlTe,dlTlL,dlTlTe) c Derivatves of Te,L w.r.t. P,T at base dlTelP=1.d0/(dlPlTe-dlPlL*(dlTlTe/dlTlL)) dlTelT=1.d0/(dlTlTe-dlTlL*(dlPlTe/dlPlL)) dlLlP=1.d0/(dlPlL-dlPlTe*(dlTlL/dlTlTe)) dlLlT=1.d0/(dlTlL-dlTlTe*(dlPlL/dlPlTe)) c Derivative of R w.r.t. P,T at base dlRlP=dlRlTe*dlTelP+dlRlL*dlLlP dlRlT=dlRlTe*dlTelT+dlRlL*dlLlT c c Radius boundary condition c s(1,4+indexv(1))=fac s(1,indexv(1))=0.d0 s(1,4+indexv(2))=-fac*dlRlP s(1,indexv(2))=0.d0 s(1,4+indexv(3))=-fac*dlRlT s(1,indexv(3))=0.d0 s(1,4+indexv(4))=0.d0 s(1,indexv(4))=0.d0 rbasel=brl(1)+(pn*fac-bpl(1))*dlRlP+(tn*fac-btl(1))*dlRlT s(1,jsf)=rn*fac-rbasel c c Luminosity boundary condition c s(2,4+indexv(1))=0.d0 s(2,indexv(1))=0.d0 s(2,4+indexv(2))=-fac*dlLlP s(2,indexv(2))=0.d0 s(2,4+indexv(3))=-fac*dlLlT s(2,indexv(3))=0.d0 s(2,4+indexv(4))=fac/fl s(2,indexv(4))=0.d0 fll=cl(1)+(pn*fac-bpl(1))*dlLlP+(tn*fac-btl(1))*dlLlT s(2,jsf)=dlog10(fl)-fll c return end c *********************************************** SUBROUTINE GETEOS (TIN,VIN,tcut,X,Y,Z,P,PE,PV,PT,ET,E,TOTLN) C *********************************************** c This routine and the other EOS routines are due c to W. Dean Pesnell. C C EQUATION OF STATE C ARGUMENTS...TIN (DEGREES K), VIN=1/RHO (CM**3/GM) C METALS... C NA,AL ALWAYS IONIZED C MG,SI,FE INCLUDED AS SINGLE ELEMENT C ALL OTHERS IGNORED c c set tcut to be log of temperature above which full c ionization is enforced C C TABLE OF RETURNED QUANTITIES. C C FUNCTION NAME DERIVATIVE WITH RESPECT TO C TEMP. SP. VOL. C------------------------------------------------------- C PRESSURE I P I PT I PV I C INT. ENERGY I E I ET I EV I C ELEC. PRS. I PE I PET I PEV I C MOLAR DENSITYI BP I BPT I BPV I C------------------------------------------------------- C C C X = HYDROGEN MASS FRACTION C Y = HELIUM MASS FRACTION C Z = METALLIC MASS FRACTION C PARGRM = MEAN MOLECULAR MOLAR DENSITY WITHOUT C ELECTRONS C C R = GAS CONSTANT 8.31434E7 C A = STEFAN-BOLTZMAN CONSTANT 7.56471E-15 C BK = BOLTZMAN S CONSTANT 8.6170837E-5 C AVAGD = AVAGADRO S NUMBER 6.02217E23 C AD3 = A/3 C IMPLICIT DOUBLE PRECISION (A-H, O-Z) PARAMETER (ZERO = 0.D0,ONE = 1.D0, TWO=2.D0, * THRE=3.D0, FOR=4.D0, TEN=10.D0, AHF=0.5D0, * QRT=0.25D0 ) PARAMETER ( R = 8.31434D7, A = 7.56471D-15, * BK = 8.6170837D-5, AVAGD = 6.02217D23, * AD3 = A/3.D0 ) DATA T3OUT,T4OUT/1.665795163D-25,3.802592017D-28/ DATA T2OUT,T5OUT/ 5.347896D-35,6.614536D-34/ C C IONIZATION POTENTIALS FOR HYDROGEN AND HELIUM C DATA XH,XHE,XHE2/13.595D0,24.581D0,54.403D0/ DATA C1,C2,C3/4.0092926D-9,1.00797D0,4.0026D0/ DATA XM,CM,ZPZP/7.9D0,0.7D0,0.12014D0/ DATA PREC / 1.D-10 / DATA ONHLF/1.5D0/ tl=dlog10(tin) V = VIN T = TIN IF( V .LE. ZERO ) GO TO 11 IF( T .LE. ZERO ) GO TO 10 FRE = ZERO ENT = ZERO PARGRM = X/C2 + Y/C3 RMUC = ONE/PARGRM RT = R*T TT4 = T**4 TK = ONE/(T*BK) SQT = DSQRT(T) C C1=ORIGINAL(C1(0.33334622))/R T1 = V*SQT**3*C1 T2 = T2OUT IF( T .GT. 2.D3 ) T2 = DEXP(-XH*TK) T3 = T3OUT IF( T .GT. 5.D3 ) T3 = DEXP(-XHE*TK) T4 = T4OUT IF( T .GT. 1.D4 ) T4 = DEXP(-XHE2*TK) T5 = T5OUT IF( T .GT. 1.2D3 ) T5 = DEXP(-XM*TK) c c check for tcut and therefore enforce complete ionization c smoothly interpolate between full ionization at tcut+0.2 c and unionized at tcut and below. Force full ionization c by setting t1 equal to about 1000. c t1max=10000d0 dlt1max=4.d0 c if (tl.lt.tcut+0.2d0.and.tl.gt.tcut) then factor=(tl-tcut)/0.2d0 t1=10.d0**(dlog10(t1)+factor*dlt1max) elseif (tl.ge.tcut+0.2d0) then t1=t1max endif c D = T1*T2 B = FOR*T1*T3 C = B*T1*T4 DD = TWO*CM*T1*T5 ZNA = Z*2.48D-3/24.969D0 ZMG = Z*ZPZP/45.807D0 C C CONVERGE ON ELECTRON DENSITY USING THE SAHA C EQUATION. C C GES IS THE MOLAR DENSITY OF ELECTRONS. C GES = (X+Y*AHF)/(ONE+Y/(FOR*C)) IF( GES .LT. X ) GES = AHF*(DSQRT(D*(D+FOR*X))-D) IF( GES .LT. 1.D-6*Z ) GES = 1.D-6*Z XC2 = X/C2 YC3 = Y/C3 C C NEWTON METHOD FOR ELECTRON DENSITY. C DO 1 I=1,25 T2 = C/GES+GES+B GEP = XC2*D/(GES+D)+YC3*(B+TWO*C/GES)/T2 @ + ZMG*DD/(GES+DD) + ZNA T1 = ONE+XC2*D/(D+GES)**2+YC3/T2* @ (TWO*C/GES**2+(B+TWO*C/GES)*(ONE-C/GES**2)/T2) @ + ZMG*DD/(GES+DD)**2 DGES = (GEP-GES)/T1 GES = GES+DGES IF( DABS(DGES)/GES .LT. PREC ) GOTO 3 1 CONTINUE GOTO 12 3 CONTINUE C C ELECTRON PRESSURE C PE = RT*GES/V C C TOTLN = 1/MU = X/C2+Y/C4+Z/C3+GES C TOTLN = PARGRM+GES XX = D/(GES+D) T2 = GES+B+C/GES YY = B/T2 ZZ = C/(GES*T2) WW = DD/(GES+DD) C C DERIVATIVES OF THE SAHA EQUATION FOR THE C PRESSURE AND INTERNAL ENERGY TEMPERATURE AND C DENSITY DERIVATIVES. C T1 = YC3*(B+TWO*C/GES) QC0 = ONE+XC2*XX/(GES+D)+ZMG*WW/(GES+DD)+YC3/T2* @ (TWO*C/GES**2+(B+TWO*C/GES)*(ONE-C/GES**2)/T2) QC1 = XC2*(ONE-XX)/(GES+D) QC4 = ZMG*(ONE-WW)/(GES+DD) QC2 = (YC3-T1/T2)/T2 QC3 = (YC3*TWO-T1/T2)/(GES*T2) QGV = (QC1*D+QC2*B+QC3*TWO*C+QC4*DD)/(QC0*V) QP1 = D*(ONHLF+XH*TK)/T QP2 = B*(ONHLF+XHE*TK)/T QP3 = C*(THRE+(XHE+XHE2)*TK)/T QP4 = DD*(ONHLF+XM*TK)/T QGT = (QC1*QP1+QC2*QP2+QC3*QP3+QC4* QP4)/QC0 C C ELECTRON PRESSURE DERIVATIVES. C PET = ONE + QGT/GES PEV =-ONE + QGV/GES C C PRESSURE DUE TO THE IDEAL GAS C P = RT*TOTLN/V PT = P/T+RT*QGT/V PV = RT*QGV/V-P/V C C BP IS R/MU C BP = R*TOTLN BPV = R*QGV BPT = R*QGT C C ADD THE RADIATION PRESSURE C P = P+AD3*TT4 PT = PT+FOR*AD3*TT4/T C C IONIZATION ENERGY C EI = (R/BK)*(XH*XX*XC2+YC3*(XHE*YY+(XHE+XHE2)*ZZ) @ +ZMG*XM*WW + ZNA*5.524D0 ) C TOTAL INTERNAL ENERGY E = ONHLF*RT*TOTLN + A*V*TT4 + EI EV = T*PT-P QXT = ((ONE-XX)*QP1-XX*QGT)/(GES+D) DT2 = QGT*(ONE-C/GES**2)+QP2+QP3/GES QYT = (QP2-B*DT2/T2)/T2 QZT = (QP3-C*QGT/GES-C*DT2/T2)/(T2*GES) QWT = ((ONE-WW)*QP4-WW*QGT)/(GES+DD) EIT = (R/BK)*(XH*QXT*XC2+YC3*(XHE*QYT+(XHE+XHE2)*QZT)+ @ ZMG*XM*QWT) ET = ONHLF*R*(TOTLN+T*QGT)+FOR*A*V*TT4/T+EIT RETURN 10 WRITE(6,100) T,V 100 FORMAT(1H0,27HNEGATIVE TEMP IN EOS T=,1PE10.3,2X, * 2HV=,E10.3) STOP 11 WRITE(6,101) T,V 101 FORMAT(1H0,27HNEGATIVE DENSITY IN EOS T=,1PE10.3,2X, * 2HV=,E10.3) STOP 12 WRITE(6,102) T,V 102 FORMAT(1H0,27HNO CONVERGENCE IN EOS T=,1PE10.3,2X, * 2HV=,E10.3) STOP END c subroutine EOS(P,t,xin,yin,zin,tcut,q,cp,dela, 1 d,xd,xt,u,cv,eta,fmu) implicit double precision (a-h,o-z) c c C EQUATION OF STATE FOR A MIXTURE OF X (HYDROGEN) Y (HELIUM) C and metals. Pressure and Temperature c are independent variables. Saha equation solver assumes c non-degenerate material, Eggleton Faulkner and Flannery c technique for Fermi-Dirac statistics... c c INPUT: c P = pressure (cgs) c T = temperature (K) c XIN = hydrogen abundance by mass c YIN = helium abundance by mass c Zin = "metal" abundance by mass c Tcut = assumed fully ionized above tcut+0.2; between Tcut c and Tcut+0.2, interpolates between full ionization c and Saha c c OUTPUT: c q = - dln(rho)/dlnT at constant pressure c Cp = dU/dT at constant pressure c dela = adiabatic temperature gradient (=dlnT/dlnP at constant S) c Ro = density (g/cm3) c Xd = dlnP/dlnrho at constant T c Xt = dlnP/dlnT at constant rho c U = internal energy per gram c Cv = dU/dT at constant volume c eta = degeneracy parameter c fmu = mean atomic weight per particle C DOUBLE PRECISION NE,dne1 fln10=2.30258509299404d0 C x=xin y=yin c For starters, set metals as the rest - i.e. ignore input z. z=1.d0-xin-yin c ptl=dlog(p)/fln10 tl=dlog(t)/fln10 iful=0 ne=0.d0 c c First compute Saha-style E.O.S. if TL tcut+0.2, or if Saha says so c or if pressure greater than 1e20 c if (tl.gt.tcut+0.2d0.or.iful.ge.1.or.ptl.gt.20.d0) then q=qi cp=cpi dudtp=dudtpi dela=delai d=di dl=dli xd=xdi xt=xti u=ui ul=uli cv=cvi eta=etai fmu=fmui c c If log(T) between Tcut and Tcut+0.2d0, interpolate between full c ionization and partial ionization for the quantities. c elseif (Tl.gt.tcut.and.tl.lt.Tcut+0.2d0) then fac=(tl-tcut)/0.2d0 q=q+fac*(qi-q) Cp=Cp+fac*(Cpi-Cp) Dela=Dela+fac*(Delai-Dela) dl=dl+fac*(dlog10(di)-dl) c d=10.d0**dl d=dexp(2.302585093d0*dl) Xd=Xd+fac*(Xdi-Xd) Xt=Xt+fac*(Xti-Xt) Ul=Ul+fac*(dlog10(Ui)-Ul) c u=10.d0**Ul u=dexp(2.302585093d0*Ul) Cv=Cv+fac*(Cvi-Cv) eta=eta+fac*(etai-eta) fmu=fmu+fac*(fmui-fmu) endif c c Trouble? Report! c if (xt.lt.0.d0.or.DELA.lt.0.d0.or.cv.lt.0.d0.or.xd.lt.0.d0) then print 131, ptl,tl,x,y,c,o,dlog10(d),xd,q,xt,DELA,cv,cp 1 ,d1,d2,d2-d1,u1,u2,u2-u1,eta1,eta write (31,131) ptl,tl,x,y,c,o,dlog10(d),xd,q,xt,DELA,cv,cp 1 ,d1,d2,d2-d1,u1,u2,u2-u1,eta1,eta 131 format(' Trouble in EOS land!'/' log P=',f9.5,'log T=',f9.5, 1 ' X,Y,C,0 :',4f7.4,/,' log ro =',f9.5/ 1 ' Xd=',f9.5,' Q=',f9.5,' Xt=',f9.5/ 2 ' DelAd=',f9.5,' Cv=',1pe10.3,' Cp=',1pe10.3/ 3 ' d1 =',e10.3,' d2=',e10.3,' d2-d1=',e10.3/ 4 ' U1 =',e10.3,' U2=',e10.3,' U2-U1=',e10.3/ 5 ' ETA1=',e10.3,' ETA=',e10.3) endif 100 format(3f8.4,1p3e12.4,0p2f9.5) c c Okay, all done! c RETURN END C SUBROUTINE STATE(PTL,TL,AX1,AX2,AX3,IFUL,D,UT,ETA,MU,NE,ITNE) IMPLICIT DOUBLE PRECISION (A-H,O-Z) c c Calculates equation of state using Saha equation a la Mihalas for c partial ionization. c C Does not include pressure ionization (except via temperature cutoffs). c C INPUTS: PTL = LOG TOTAL PRESSURE C TL = LOG TEMPERATURE C AX1 = MASS FRACTION OF HYDROGEN C AX2 = MASS FRACTION OF HELIUM C AX3 = MASS FRACTION OF 'METAL' c NE = FIRST GUESS AT ELECRON DENSITY c C OUTPUTS: C D = DENSITY C UT = INTERNAL ENERGY PER GRAM C ETA = DEGENERACY PARAMETER c Mu = MEAN ATOMIC WEIGHT PER PARTICLE C NE = NUMBER DENSITY OF ELECTRONS C ITNE = NUMBER OF ITERATIONS ON NE C IFUL = SET TO 1 IF FINDS COMPLETE IONIZATION C DOUBLE PRECISION Y(0:10,5),NE,MUI,MUE,X(3),MU DOUBLE PRECISION Netot COMMON/ATDAT/G(0:2,3),XI(0:2,3),A(3),IZ(3) 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 Y/55*0.D0/ c CSAHA=(2.D0*CPI*CK*CME)**(3.D0/2.D0) c CSAHA=CSAHA/CH/CH/CH csaha=2.414703187d15 C NSPEC=3 X(1)=AX1 X(2)=AX2 X(3)=AX3 c c T=10.D0**(TL) T=dexp(2.302585093d0*TL) c Ptot=10.D0**(PTL) Ptot=dexp(2.302585093d0*PTL) Prad=2.5219D-15*(T**4) Pgas=PTOT-PRAD if (pgas.le.0.d0) then print *,' Total pressure less than gas pressure?!' print '(A,f10.5,a,f10.5)','log P=',ptl,'log T=',tl print *,' Setting it equal to 0.01*P, but BEWARE' pgas=0.01d0*ptot endif PgasL=dlog10(Pgas) THRES=1.d-5 c If X+Y+Z doesn't add up to one, fudge with a heavy, atomic weight c of 20 and charge of 1. sumx=x(1)+x(2)+x(3) if (sumx.lt.1.d0) then X(3)=x(3)+1.d0-sumx endif C C COMPUTE . . . C AVZOA, the mass-weighted c MUI, the mean atomic weight per nucleus C MUE, the mean atomic weight per electron, free or bound OOMUI=0.D0 AVZOA=0.D0 DO 5 I=1,NSPEC OOMUI=OOMUI+X(I)/A(I) AVZOA=AVZOA+X(I)*IZ(I)/A(I) 5 CONTINUE c MUI=1.D0/OOMUI MUE=1.d0/avzoa C C COMPUTE NE C Use the Saha equation to find the new NE without pressure ionization c CALL SAHA(PGASL,DL,TL,NSPEC,X,NE,Y,ITNE) D=MUI*(Pgas/cna/ck/t-Ne/cna) dl=dlog10(d) c Now is a good time to compute Mu, the mean atomic weight per free particle fmuefr=D*CNA/NE mu=1.d0/(1.d0/mui+1.d0/fmuefr) c c Check for "complete" ionization c NETOT=D/MUe*CNA c print '(10hNe/Netot ,1pe12.5)',ne/netot IF (NE/NETOT.gt.0.999D0) IFUL=1 C C CACULATE INTERNAL ENERGY PER GRAM C C Calculate properties of e- gas. PPGE=NE*CK*T UPGE=1.5d0*ppge/d C RADIATION URAD=3.D0*PRAD/D C TRANSLATIONAL PPGI=D*CNA*CK*T/MUI UGT=UPGE+1.5D0*PPGI/D C IONIZATION ENERGY FOR EACH SPECIES UGI=0.D0 DO 50 K=1,NSPEC IF (X(K).GE.1.D-10) THEN UGII=0.D0 DO 60 IJ=1,IZ(K) DO 70 JJ=0,IJ-1 UGII=UGII+X(K)*Y(IJ,K)*XI(JJ,K)*CEV 70 CONTINUE 60 CONTINUE UGII=UGII*CNA/A(K) UGI=UGI+UGII ENDIF 50 CONTINUE c c Add in dissosiation energy of H2, assume always dissociated XIH2=4.477d0 UH2=0.5d0*X(1)*CEV*XIH2*CNA UGI=UGI+UH2 C TOTAL INTERNAL ENERGY PER GRAM c print '(1pe12.5/e12.5/e12.5/e12.5)',upge,1.5d0*ppgi/d,ugt,ugi UT=URAD+UGT+UGI c print '(1p5e14.6)',upge,1.5d0*ppgi/d,ugt,ugi,urad c Degeneracy parameter in nondegenerate limit ETA=-36.113d0+dlog(Ne)-3.4539d0*TL c RETURN END C subroutine SAHA(PGL,DL,TL,NSPEC,X,NE,Y,ITNE) implicit double precision (a-h,o-z) c c Subroutine to compute Ne for an ensemble of NSPEC atomic species c NSTAGE in the parameter statement is the maximum number of ionization c stages expected. Can be called with either (Pgas,T) or (Rho,T) c as the independent variables. c c If the sum of the abundances is less than 1, it adds a final species c that is a one-electron metal c c Follows algorithm in Mihalas (around p. 118) but with (Pg,T) as the c known quantities. Thus the linearization is slightly different, c but all notation is the same as Mihalas. c c Input: PGL = log gas pressure c note: if PGL is input as greater than 90 c then this routine c uses DL as the first independent variable c DL = log density c TL = log temperature c NSPEC = number of species to consider c X = array of mass fractions c c Output: Ne = number density of free electrons c Y(Nstage,NSPEC)= fraction of given stage that is ionized c ITNE = number of iterations performed c parameter(NSTAGE=10) double precision Ne,Ne0,Ntot,Netot,NNuc dimension Y(0:10,5),X(5) dimension phi(0:10,5),prophi(5) dimension P(0:nstage,5),S(5),alfa(5) COMMON/ATDAT/G(0:2,3),XI(0:2,3),A(3),IZ(3) 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 first index in, i.e. XI(J,K) is ionization stage, second is species c c read in atomic and other parameters if this is the first call c CSAHA=CH**2/2.d0/CPI/CME/CK CSAHA=0.5d0*dexp(1.5d0*dlog(CSAHA)) CKEV=CK/CEV THRESH=1.d-5 C c now to work! c c T=10.d0**TL T=dexp(2.302585093d0*TL) FKT=CKEV*T c TM32=10.d0**(-3.d0*TL/2.d0) TM32=dexp(-3.d0*(TL*2.302585093d0)/2.d0) CSTM32=CSAHA*TM32 c compute total number fraction of each species UAm1=0.d0 nK=nspec do 5 k=1,nspec UAm1=UAm1+x(k)/A(k) 5 continue do 7 k=1,nK alfa(k)=X(k)/A(k)/UAm1 7 continue c print '(1p5e12.5)',(alfa(k),k=1,nK) c c if Pgas is an input variable, compute Ntot, the total c number density of all paritcles, from the Pressure c if (PGL.lt.90.d0) then c Pgas=10.d0**PGL Pgas=dexp(2.302585093d0*PGL) Ntot=Pgas/CK/T if (ntot.le.0.d0) ntot=1.d5 c and use this to make a first guess at Ne0 Ne0=Ntot/100.0 else c c otherwise, use the input density to determine c the number density of nuclei c rho=10.d0**DL rho=dexp(2.302585093d0*DL) Nnuc=rho*uam1*cna c and the total number of electrons, free or bound fnetot=0.d0 do 9 k=1,nk fnetot=fnetot+alfa(k)*iz(k) c print '(i4,2f10.4,1pe12.4,0pi4,1pe12.4)', c 1 k,pgl,tl,alfa(k),iz(k),fnetot 9 continue netot=nnuc*fnetot Ne0=Netot/100.d0 endif c c More preliminaries: phi's and stuff do 115 k=1,nk if (x(k).le.1.d-10) go to 115 prophi(k)=1.d0 do 125 j=iz(k)-1,0,-1 flpjk=dlog(g(j,k)/g(j+1,k)*csaha*tm32)+(xi(j,k)/fkt) c don't let the exponential factor get too big... causes trouble in unionized c material... if (flpjk.gt.28.d0) flpjk=28 phi(j,k)=dexp(flpjk) prophi(k)=prophi(k)*phi(j,k) 125 continue 115 continue C C ITERATE ON NE: c compute the P's and S's. Note Mihalas's convention that the c product terms involving the last ionization stage ie. P(iz(k),k) c are unity. c itne=0 verg=1.d0 do while (itne.le.40.and.verg.ge.thresh) itne=itne+1 sumbar=0.d0 c dsbdNe is the derivative of sumbar w.r.t. Ne0 dsbdNe=0.d0 do 15 k=1,nK S(k)=1.d0 if (x(k).le.1.d-10) go to 15 P(iz(k),k)=1.d0 sumjP=iz(k)*1.d0 sumjdP=0.d0 dSkdNe=0.d0 do 25 j=iz(k)-1,0,-1 P(j,k)=P(j+1,k)*ne0*phi(j,k) sumjP=sumjP+j*P(j,k) c dPjkdNe is the derivative of P(J,K) w.r.t. Ne0 dPjkdNe=(iz(k)-j)*ne0**(iz(k)-j-1)*prophi(k) sumjdP=sumjdP+j*dPjkdNe S(k)=S(k)+P(j,k) dSkdNe=dSkdNe+dPjkdNe 25 continue sumbar=sumbar+alfa(k)*sumjP/S(k) dsbdNe=dsbdNe+alfa(k)*(sumjdP/S(k)-sumjP/S(k)*dSkdNe/S(k)) c 15 continue c c the new computed value of Ne c if (PGL.lt.90.d0) then ne=(Ntot-ne0)*sumbar else ne=nnuc*sumbar endif c c print *,iter,' HII', P(1,1)/S(1),' HI',P(0,1)/S(1) c print *,iter,' HeII', P(1,2)/S(2),' HeI',P(0,2)/S(2) c print *,iter,' HeIII', P(2,2)/S(2) c print *,iter,' XMII', P(1,3)/S(3),'XMI',P(0,3)/S(3) c print * c c Compute a new (hopefully better) value for Ne: c Linear correction to ne0 c if (PGL.lt.90.d0) then dne=(ne-ne0)/(1.d0+sumbar-(Ntot-ne0)*dsbdNe) else dne=(ne-ne0)/(1.d0-nnuc*dsbdNe) endif c print '(I3,1p3e13.5,3x,2e12.5)',iter,ne0,ne,dne/ne0,netot c c CONVERGENCE CRITERIA: c difference between computed and guessed value of ne verg=dabs((ne-ne0)/ne0) c electron density/total particle density if (ntot.le.0.d0) then frace=0.5d0 else frace=dabs(ne0/ntot) endif verg=dmin1(frace*100.d0,verg) c Add the correction. c if the correction is too large, reduce it a tad. This is essential c when ne/netot is tiny and the correction tries to subtract the c value of ne0 from itself (to within machine accuracy) fac=dabs(dne/ne0) if (fac.ge.0.1d0) dne=(1.d0-1.d-6)*dne c and make the correction ne0=ne0+dne if (ne0.le.0.d0) ne0=1.d0 c If iteration counter is getting too large, then reduce convergence c threshold. if (itne.eq.20) then write (31,131) PGL,TL,X(1),X(2),X(3) write (31,'(I3,1p3e13.5,3x,e12.5)') iter,ne0,ne,dne/ne0 131 format(' Saha routine in trouble : PGL=',f9.4,' TL=',f9.4, 1 ' X=',f7.4,' Y=',f7.4,' C=',f7.4/ 'Reducing thresh') thresh=thresh*10.d0 endif if (itne.eq.30) then write (31,132) PGL,TL write (31,'(I3,1p3e13.5,3x,e12.5)') iter,ne0,ne,dne/ne0 132 format(' Saha routine in trouble again : PGL=',f9.4, 1 ' TL=',f9.4) thresh=thresh*10.d0 endif c ******** End of iteration loop ***************** enddo c c If electron density small, return with a guess at it based on slight c ionization if (frace.le.thresh) ne0=dsqrt(Pgas/ck/t) c store the new value of ne in the correct variable ne=ne0 c c compute fractional ionization do 28 k=1,nspec do 27 j=0,iz(k) y(j,k)=p(j,k)/s(k) 27 continue 28 continue c return end C block data c SUBROUTINE RATDAT C C ATOMIC DATA FOR THE CHOSEN SPECIES C IMPLICIT DOUBLE PRECISION (A-H,O-Z) COMMON/ATDAT/G(0:2,3),XI(0:2,3),A(3),IZ(3) DATA A/1.00782505d0,4.00260308d0,14.d0/ DATA IZ/1,2,1/ DATA G/2.d0,1.d0,0.d0, 1 1.d0,2.d0,1.d0, 2 2.d0,1.d0,0.d0/ DATA XI/13.598d0,0.d0,0.d0, 1 24.587d0,54.416d0,0.d0, 2 5.000d0,0.0d0,0.0d0/ c a stub c temp=1234.56d0 c return c END c subroutine IONIZED(PL,TL,X,Y,z,d,U,eta,fmu,Pe,Pi,ue, 1 Xd,q,dudtP) IMPLICIT DOUBLE PRECISION (A-H,O-Z) C C FULLY IONIZED EQUATION OF STATE: ITERATE ON DENSITY TO GIVE PROPER c TOTAL PRESSURE GIVEN THAT ALL ELECTRONS ARE FREE c C INPUT: PL = log(total pressure) c TL = log(temperature) c X = hydrogen mass fraction c Y = helium mass fraction c C = carbon mass fraction c c OUTPUT: d = density c U = total energy density c eta = degeneracy parameter c fmu = mean atomic weight per free particle c Pe = electron pressure c Pi = ion pressure c Ue = energy density in electron gas c Xd = dlnP/dlnrho at constant T c q = -dlnrho/dlnT at constant P c dudtP = dU/dT at constant P C c Common block CONST contains fundamental physical constants COMMON/CONST/CMS,CRS,CLS,CG,CSB,CPI,CP4,CCL,CME,CH,CNA,CEV,CK,CSYR DOUBLE PRECISION mue,mui c c Convergence threshold for fully ionized iterations c THRESH=1.d-5 c c Some preliminaries: subtract off radiation pressure to obtain c gas pressure alone c Ptot=10.d0**pl Ptot=dexp(2.302585093d0*pl) c T=10.d0**tl T=dexp(2.302585093d0*tl) Prad=2.5219d-15*t**4 Pgas=Ptot-Prad c Assume that anything left after x,y,c is o Z=1.d0-X-Y c Mean atomic weight per ion and per electron (free or bound) oomui=X/1.007825d0 + Y/4.002603d0 + Z/14.0d0 oomue=X/1.007825d0 + 2.d0*Y/4.002603d0 + 0.5d0*Z mue=1.d0/oomue mui=1.d0/oomui fmu=1.d0/(oomui+oomue) c Sum of ionization potentials in eV, for energy computation sumxix=13.598d0 sumxiy=79.003d0 sumxiz=1500.d0 C First guess at D0 and NE: assume nondegenerate gas with no Coulomb term D0=Pgas/CNA/CK/T/(1.D0/mui+1.D0/mue) c Don't allow the first guess of D0 to exceed 1e6, otherwise relativistic c corrections go haywire if (d0.gt.1.d6) d0=1.d6 c ************************************************************************ c C Iterations on D to get Pe+Pi=Pgas; continue iterating until density c changes by less than THRESH and derived Pgas matches true value c drho=d0 itfi=0 DO while ((dabs(drho/d0).gt.THRESH 1 .OR.dabs(dpgas/pgas).gt.THRESH) 1 .AND.itfi.lt.30) itfi=itfi+1 c Pressure of ions assumed to be a perfect gas Pi=d0/mui*CNA*CK*T C Compute pressure of electron gas given guess at domue via EFF DOmue0=D0/mue call EFF(domue0,T,eta,Pe,Ued,dpedom,dpedt,duedtr) c Compute new value of total gas pressure with this guess c at the density. Pgas1=Pe+Pi dPgas=Pgas-Pgas1 c Compute the corresponding linearized correction to the density c dPe/drho: dpedro=dpedom/mue c dPi/drho: dpidro=CNA*CK*T/mui c Use these derivatives to compute the linearized density correction fac=dpidro+dpedro drho=dpgas/fac c reduce the density correction if it is very large if (dabs(drho/d0).ge.1.d0) then c write (31,*)' Reducing drho by 50%, T=',t,' Rho=',d0,' Iter=' c 1 ,itfi drho=dsign(d0/2.d0,drho) elseif (dabs(drho/d0).gt.0.1d0) then c write (31,*)' Reducing drho by 10%, T=',t,' Rho=',d0,' Iter=' c 1 ,itfi drho=drho*0.9d0 endif c Apply the density correction d1=d0+drho c print *,itfi,d1,drho/d1 d0=d1 c enddo c C Done iterating on D for full ionization c ******************************************************************** c All done; set up for return d=d0 IF (ITFI.GT.29) then PRINT *, ' Fully Ionized Iterations failed to converge' print *, ' T=',dlog10(t),' P=',dlog10(PL) endif c c Internal energy c c Radiation pressure energy density Urad=3.d0*Prad/d c Energy DENSITY of electron gas Ue=Ued/d c Translational energy density of ions Ui=1.5d0*Pi/d c Excitation energy density Uex=cna*(X*sumxix+Y*sumxiy/4.d0+z*sumxiz/14.d0) Uex=Uex*cev c Total energy U=Urad+Ue+Ui+Uex c c Compute derivatives analytically c pressure derivatives dpidt=d*cna*ck/mui dpidro=cna*ck*t/mui dpedro=dpedom/mue dprdro=0.d0 dprdt=4.d0*prad/t c dpdro=dpidro+dpedro+dprdro dpdt=dpidt+dpedt+dprdt drodt=-dpdt/dpdro c c Energy derivative w.r.t. T at constant rho duidt=1.5d0*cna*ck/mui c c print *,"Gam=",gam c duedt=duedtr/d durdt=4.d0*3.d0*prad/t/d c dudtr=duidt+duedt+durdt dudtp=dudtr-t/d**2*dpdt*drodt+ptot/d**2*drodt c ********************************************************************** c Derivatives for output xd=dpdro/ptot*d xt=dpdt/ptot*T q=xt/xd return end c subroutine EFF(domue,T,psi,Pe,Ue,dPeDoM,dPedT,dUedTr) implicit double precision (a-h,o-z) c dimension rhat(4,4),phat(4,4),uhat(4,4) data rhat/2.315472d0,7.128660d0,7.504998d0,2.665350d0, 1 7.837752d0,23.507934d0,23.311317d0,7.987465d0, 2 9.215560d0,26.834068d0,25.082745d0,8.020509d0, 3 3.693280d0,10.333176d0,9.168960d0,2.668248d0/ data phat/2.315472d0,6.748104d0,6.564912d0,2.132280d0, 1 7.837752d0,21.439740d0,19.080088d0,5.478100d0, 2 9.215560d0,23.551504d0,19.015888d0,4.679944d0, 3 3.693280d0,8.859868d0,6.500712d0,1.334124d0/ data uhat/3.473208d0,10.122156d0,9.847368d0,3.198420d0, 1 16.121172d0,43.477194d0,37.852852d0,10.496830d0, 2 23.971040d0,60.392810d0,47.782844d0,11.361074d0, 3 11.079840d0,26.579604d0,19.502136d0,4.002372d0/ M=4 N=4 c c given Domue and T, compute rho/mue, Pe, Ue using the thermodynamically c consistent interpolation formulae of Eggleton Faulkner and Flannery (1973) c also returns analytic values of d(Pe)/d(rho/mue) at constant T, c d(Pe)/dT at constant rho, and d(Ue)/dT at constant rho c c note that Ue is internal energy per unit VOLUME c flcomp is the compton wavelength of the electron c flcomp=2.42580d-10 flcomp=2.42631d-10 fac=flcomp**3/8.d0/3.1415927d0/9.1094d-28/2.998d10**2 c beta=t/9.1094d-28/2.998d10**2*1.3807d-16 rhoh0=domue*flcomp**3/8.d0/3.141592654d0/1.66054d-24 c c Obtain f from domue via Newton's method using EFF eq. 15 c c f0=domue/beta**1.5*1.693d-7 f0=domue/dexp(1.5d0*dlog(beta))*1.693d-7 do 10 it=1,20 f=f0 g=beta*dsqrt(1.d0+f) dbdfg=-beta/2.d0/(1.d0+f) dbdgf=beta/g dgdfb=beta/2.d0/dsqrt(1.d0+f) c Now evaluate the polynomial for the rhohat corresponding to c the guessed f (and g) sumr=0.d0 sumdrf=0.d0 sumdrg=0.d0 do 20 im=1,M do 30 in=1,N fmgn=f**(im-1)*g**(in-1) sumr=sumr+rhat(im,in)*fmgn sumdrf=sumdrf+rhat(im,in)*fmgn*(im-1.d0)/f sumdrg=sumdrg+rhat(im,in)*fmgn*(in-1.d0)/g 30 continue 20 continue fnorm=(1.d0+f)**(M-1)*(1.d0+g)**(N-1) c g32=dexp(1.5d0*dlog(g)) ggp132=dexp(1.5d0*dlog(g+g*g)) c rhohat=f/(1.d0+f)*g**1.5d0*(1.d0+g)**1.5d0*sumr/fnorm c rhohat=f/(1.d0+f)*g32*gp132*sumr/fnorm rhohat=f/(1.d0+f)*ggp132*sumr/fnorm drdfg=rhohat*(1.d0/f/(1.d0+f)+sumdrf/sumr-(m-1.d0)/(1.d0+f)) drdgf=rhohat* 1 (1.5d0/g+1.5d0/(1.d0+g)+sumdrg/sumr-(n-1.d0)/(1.d0+g)) drdfb=(drdfg*dbdgf-drdgf*dbdfg)/dbdgf c compute the correction to f delta=(rhohat-rhoh0)/drdfb c print *,it,f,delta/f f0=f-delta c if the correction is small enough, then u.r. done if (dabs(delta/f).lt.1.d-7) go to 11 10 continue c 11 if (it.ge.20) then print *, ' E.F.F. iterations failed to converge' write (31,*) ' E.F.F. iterations failed to converge' write (31,*) ' log rho/mue=',dlog10(domue),' log T=',dlog10(t) endif c f=f0 c if (f.lt.0.d0) print *,' EFF: f=',f0 c if (f.lt.0.d0) f=1.d-10 g=beta*dsqrt(1.d0+f) ggp132=dexp(1.5d0*dlog(g+g*g)) dbdfg=-beta/(2.d0*(1.d0+f)) dbdgf=beta/g c c now use the value of f to compute Pe and Ue and Psi c sump=0.d0 sumu=0.d0 sumdpf=0.d0 sumdpg=0.d0 sumduf=0.d0 sumdug=0.d0 do 40 im=1,M do 50 in=1,N fmgn=f**(im-1)*g**(in-1) sump=sump+phat(im,in)*fmgn sumdpf=sumdpf+phat(im,in)*fmgn*(im-1.d0)/f sumdpg=sumdpg+phat(im,in)*fmgn*(in-1.d0)/g sumu=sumu+uhat(im,in)*fmgn sumduf=sumduf+uhat(im,in)*fmgn*(im-1.d0)/f sumdug=sumdug+uhat(im,in)*fmgn*(in-1.d0)/g 50 continue 40 continue fnorm=(1.d0+f)**(M-1)*(1.d0+g)**(N-1) c x=dsqrt(1.d0+f) psi=2.d0*x+dlog((x-1.d0)/(x+1.d0)) c g52=dexp(2.5d0*dlog(g)) c gp132=dexp(1.5d0*dlog(g+1.d0)) g52p132=g*ggp132 c PeHat=f/(1.d0+f)*g**2.5d0*(1.d0+g)**1.5d0*sump/fnorm c UeHat=f/(1.d0+f)*g**2.5d0*(1.d0+g)**1.5d0*sumu/fnorm c PeHat=f/(1.d0+f)*g52*gp132*sump/fnorm PeHat=f/(1.d0+f)*g52p132*sump/fnorm c UeHat=f/(1.d0+f)*g52*gp132*sumu/fnorm UeHat=f/(1.d0+f)*g52p132*sumu/fnorm c c compute derivatives w.r.t. f at constant g dpdfg=PeHat*(1.d0/f/(1.d0+f)+sumdpf/sump-(m-1)/(1.d0+f)) dudfg=UeHat*(1.d0/f/(1.d0+f)+sumduf/sumu-(m-1)/(1.d0+f)) c and derivatives w.r.t. g at constant f dpdgf=PeHat*(2.5d0/g+1.5d0/(1.d0+g)+sumdpg/sump-(n-1)/(1.d0+g)) dudgf=UeHat*(2.5d0/g+1.5d0/(1.d0+g)+sumdug/sumu-(n-1)/(1.d0+g)) c compute due/dT at constant rho = duhhat/dbeta*8*pi*k/flam**3 dudbr=(dudfg*drdgf-dudgf*drdfg)/(dbdfg*drdgf-dbdgf*drdfg) dUedTr=dudbr*8.d0*3.14159254d0*1.38d-16/flcomp**3 c now dPe/d(rho/me) at constant beta dPhDrh=(dpdgf*dbdfg-dpdfg*dbdgf)/(dbdfg*drdgf-drdfg*dbdgf) dPedOm=dphdrh*9.109d-28*2.998d10**2/1.6724d-24 c and dPe/dT:rho=(dPh/dbeta)*8*pi*k/flam**3 dphdb=(dpdfg*drdgf-dpdgf*drdfg)/(dbdfg*drdgf-dbdgf*drdfg) dPedT=dPhdb*8.d0*3.141592654d0*1.38d-16/flcomp**3 c Pe=PeHat/fac Ue=UeHat/fac c return end c subroutine iocond(dlin,tlin,xin,yin,ochl) implicit double precision (a-h,o-z) c c Returns the conductive opacity for a mixture (X,Y,Z=1-X-Y) at c a given value of log density and log temperature. Uses the Iben c fit to the Hubbard and Lampe opacity tables (Ap. J. 196, 545). c Implementation by S.Kawaler, August 1990. Checked against c actual H&L tables for helium and carbon. Well within 20% c agreement for log rho>0, can be way off at lower densities, so c caveat emptor... c dl=dlin tl=tlin x=xin y=yin cln=dlog(10.d0) c c If log rho is greater than 6, use Lamb(1974)'s extrapolation c of the tables for the conductive opacity when relativistically c degenerate... Though I suppose you could use the Itoh et al. c opacities in this regime... c if (dlin.gt.6.d0) then ochl1=-1.d0 dl=6.d0 endif c coarse approximation to and ue: ue=1.d0/(0.5d0*(1.d0+X)) avz2oa=x+y+(1.d0-x-y)*3.d0 c (the statement number here is for branching to when dlin > 6) 10 fldelt=dl-dlog10(ue)-1.5d0*(tl-6.d0) c delta=10.d0**fldelt delta=dexp(cln*fldelt) fleta=-0.55255+(2.d0/3.d0)*fldelt c eta=10.d0**fleta eta=dexp(cln*fleta) c Iben equations (A2)-(A5) if (fldelt.lt.0.645d0) then fltokt=-3.2862d0+dlog10(delta*(1.d0+0.024417*delta)) else a1=-3.29243d0+dlog10(delta*(1.d0+0.02804*delta)) if (fldelt.lt.2.0d0) then fltokt=a1 else b1=-4.80946d0+dlog10(delta*delta*(1.d0+9.376d0/eta/eta)) if (fldelt.gt.2.5d0) then fltokt=b1 else fltokt=2.d0*a1*(2.5d0-fldelt)+2.d0*b1*(fldelt-2.d0) endif endif endif c (A8) thru (A10) for Pe/NkT a2=dlog10(1.d0+0.021876*delta) if (fldelt.lt.1.5d0) then peonkt=a2 else b2=dlog10(0.4d0*eta*(1.d0+4.1124d0/eta/eta)) if (fldelt.gt.2.0d0) then peonkt=b2 else peonkt=2.d0*a2*(2.d0-fldelt)+2.0d0*b2*(fldelt-1.5d0) endif endif c (A11) and (A12) for nednedeta if (delta.lt.40.d0) then dlnede=1.d0-0.01d0*delta*(2.8966d0-0.034838d0*delta) else dlnede=(1.5d0/eta)*(1.d0-0.8225/eta/eta) endif c and (A7) for alfa temp=dlog10(ue*avz2oa+dlnede) alfa=-2.033983d0+fldelt-0.5d0*(tl-6.d0)-peonkt+temp c c Now calculate the value of Theta c c for Hydrogen if (x.le.0.d0) then flthx=0.d0 else if (alfa.le.-3.d0) then flthx=1.048d0-0.124d0*alfa elseif (alfa.gt.-1.d0) then flthx=0.185d0-0.585d0*alfa c note that .585 in above is from Iben's code; paper says 0.558 else flthx=0.13d0-alfa*(0.745d0+0.105d0*alfa) endif endif c for Helium if (y.le.0.d0) then flthy=0.d0 else if (alfa.le.-3.d0) then flthy=0.937d0-0.111d0*alfa elseif (alfa.gt.0.d0) then flthy=0.24d0-0.6d0*alfa else flthy=0.24d0-alfa*(0.55d0+0.0689d0*alfa) endif endif c for Carbon if (alfa.lt.-2.5d0) then flthc=1.27d0-0.1d0*alfa elseif (alfa.gt.0.5d0) then flthc=0.843d0-0.785d0*alfa else flthc=0.727d0-alfa*(0.511d0+0.0778d0*alfa) endif c Zhere=1.d0-x-y c ochl=x*10.d0**flthx+y*10.d0**flthy+zhere*10.d0**flthc ochl= x*dexp(cln*flthx)+y*dexp(cln*flthy)+zhere*dexp(cln*flthc) c ochl=ochl/10.d0**(tl-6.d0)/10.d0**fltokt c ochl=ochl/10.d0**(tl-6.d0+fltokt) ochl=ochl/dexp(cln*(tl-6.d0+fltokt)) c c If input density greater than 10^6, then extrapolate H&L opacities c linearly using the opacity at 10^5 and 10^6 c if (dlin.gt.6.d0.and.ochl1.le.0.d0) then ochl1=ochl dl=5.d0 go to 10 elseif (dlin.gt.6.d0.and.ochl1.gt.0.d0) then flochl=dlog10(ochl1)+(dlog10(ochl1)-dlog10(ochl))*(dlin-6.d0) ochl=10.d0**flochl return endif c return end c SUBROUTINE iorad(DL,TL,X,Y,Z,C12,X16,FKRAD) IMPLICIT DOUBLE PRECISION (A-H,O-Z) C C Compute radiative opacity using Iben's 1975 fits to C the Cox+Stewart opacities, in conjunction with the Christy opacity c fits C C INPUT: DL = log density (g/cm**3) c TL = log temperature (K) c X = hydrogen mass fraction c Y = He4 mass fraction c Z = mass fraction of all "metals" c C12 = C12 mass fraction c X14 = N14 mass fraction c X16 = O16 mass fraction c c OUTPUT: fkrad = radiative opacity c dlkrro = d(log kappa rad)/d(log rho) at constant T c dlkrt = d(log kappa rad)/d(log T) at constant rho c dimension Xsaha(5),Ysaha(11,5) FLT6=TL-6.d0 T6=10.d0**(FLT6) T=T6*1.d6 RHO=10.d0**(DL) Zresid=1.d0-X-Y-C12-X16 FMU = X+Y+3.d0*C12+4.d0*X16+(4.528d0)*Zresid - 1.d0 ZB=C12+X16 c Electron scattering opacity D=0.05d0*(FLT6-1.7d0) FKE=(0.2d0-D-DSQRT(D*D+0.0004d0))*(1.d0+X) c c Iben Equations (B3)-(B9) A=(1.25d0+0.488d0*dsqrt(FMU)+0.092d0*FMU)/0.67d0 B=3.86d0+0.252d0*dsqrt(FMU)+0.018d0*FMU C=(2.019d-4*RHO/(T6**1.7d0)) C=C**2.425d0 AP=1.d0+C*(1.d0+C/24.55d0) FLDB=-A+B*FLT6 c Iben Equation (B2) FLK1=0.67d0*(DL-FLDB)+dlog10(AP) FK1=10.d0**FLK1 c c Iben Equations (B12)-(B20) c c note that FLKRHO below is from Iben's code directly and doesn't c appear in the paper DLT=0.15d0*(0.13d0+0.212d0*X+FLT6) FLKRHO=-0.174d0-DLT+DSQRT(DLT*DLT+0.002809d0) FKRHO=10.d0**FLKRHO FLD1=-2.3d0+4.833d0*FLT6 FLD2=-(1.86d0+0.3125d0*X)+3.86d0*FLT6 FLDB=DMIN1(FLD1,FLD2) IF (T6.GT.1.d0) THEN FLD0=-(1.68d0+0.35d0*X)+1.8d0*FLT6 ELSE FLD0=-(1.68d0+0.35d0*X)+(3.42d0-0.52d0*X)*FLT6 ENDIF FAC=FLT6-0.22d0+0.1375d0*X A1=1.22d0*dexp(-(1.74d0-0.755d0*X)*(FAC*FAC)) FAC=FLT6-0.18d0+0.1625d0*X W=4.05d0*dexp(-(0.306d0-0.04125d0*X)*(FAC*FAC)) AZ0=A1*dexp(-2.76d0*(DL-FLD0)*(DL-FLD0)/(W*W)) AZ=AZ0*(Z+ZB)/0.02d0 c Iben Equation (B12) FLK2=FKRHO*(DL-FLDB)+AZ FK2=10.d0**FLK2 c c Christy radiative opacity. Note: need PE for this part. Here we c use the GETEOS routine to get it... c The Christie opacity begins to suck at densities of c about 1.d-3 and above when T is of order 1e5. c Can be off by a factor of two before that spike... c IF (T6.le.1.5d0) then Xsaha(1)=X Xsaha(2)=Y V=1.d0/RHO call GETEOS(T,V,5.5d0,x,y,z,p,pe,pv,pt,et,e,totln) v12=dsqrt(v) v14=dsqrt(v12) c-- The commented lines here are Icko's original version------- c t1=t/1.d4 c ses=(5.4d-13)*v/t1 c trt=sqrt(t1) c t2=t1*t1 c t3=t2*t1 c t4=t3*t1 c d1=20.+5.*t3+t4 c sz=z/(trt*d1) c t46=t4*t2 c d1=(1.4d3)*t1 c d2=d1+t46 c d3=1.d6+0.1d0*t46 c sy=1./d2+1.5/d3 c sy=y*sy c t10=t46*t4 c d1=2.d6+2.1*t10 c s1x=trt*t4/d1 c t4v14=t4*v14 c d1=1.d0+0.05d0*t4v14 c d2=4.5*d1*t3 c d3=d2+250.d0 c s2x=d1/(t3*d3) c sx=(s1x+s2x)*x c fkch=ses+sx+sy+sz c fkch=pe*fkch c ---------------------------------------------------------------- T4=T/1.d4 FLNT4=dlog(T4) T44=dexp(4.d0*FLNT4) T45=dexp(5.d0*FLNT4) T46=dexp(6.d0*FLNT4) c TEMPX=4.0d-3/T44+2.0d-4*v14 TEMPX=1.d0/(4.5d0*T46+1.d0/T4/TEMPX) TEMPX=X*(TEMPX+dsqrt(T4)/(2.0d6/T44+2.1d0*T46)) TEMPY=1.d0/(1.4d3*T4+T46)+1.5d0/(1.0d6+0.1d0*T46) TEMPY=Y*TEMPY TEMPZ=Z*dsqrt(T4)/(20.d0*T4+5.d0*T44+T45) FKCH=PE*(5.4d-13*V/T4+TEMPX+TEMPY+TEMPZ) ELSE FKCH=0.d0 ENDIF c IF (X.GT.0.D0.AND.T6.GE.1.5D0) THEN c print *," Case 1" FKRAD=FKE+FK2 ELSEIF (X.GT.0.d0.AND.T6.LE.1.d0) THEN c print *," Case 2" FKRAD=FKCH ELSEIF (X.GT.0.d0.AND.T6.LT.1.5d0.AND.T6.GT.1.d0) THEN c print *," Case 3" FKRAD=2.d0*FKCH*(1.5d0-T6)+2.d0*(FK2+FKE)*(T6-1.d0) ELSEIF (X.EQ.0.d0.AND.ZB.GT.Z) THEN c print *," Case 4" FKRAD=FKE+FK1 ELSEIF (X.EQ.0.d0.AND.ZB.LE.Z.AND.T6.GE.1.d0) THEN c print *," Case 5" FKRAD=(FKE+FK1)*ZB/Z+(FKE+FK2)*(1.d0-ZB/Z) ELSE print *,' Tl=',TL,' Dl=',dl print *,' X=',X,'Y=',Y,'C=',C12,'O=',O16 PAUSE "Impossible logic in Iben-style Opacity" ENDIF RETURN END c subroutine LUDCMP(a,n,np,indx,d) implicit double precision (a-h,o-z) c c Given an NxN matrix A, with dimension NP, this routine replaces it by c the LU decomposition ofa rowwise permutation of itself. A and N are c input. A is the output, arranged in equation 2.3.14 on p. 34. INDX is c an output vector which records the row permutation resulting from c partial pivoting. D is output as +/- 1 depending on whether the number c of row interchanges was even or odd, respectively. Used in combination c with LUBKSB to solve linear equations or invert a matrix. c c Numerical Recipes, p. 38-9. c INTEGER n,np,indx(n),NMAX,i,imax,j,k DOUBLE PRECISION d,a(np,np),TINY PARAMETER(nmax=500, tiny=1.d-40) DOUBLE PRECISION aamax,dum,sum,vv(NMAX) d=1.0 do 12 i=1,n aamax=0.d0 do 11 j=1,n if (abs(a(i,j)).gt.aamax) aamax=abs(a(i,j)) 11 continue if (aamax.eq.0.d0) PAUSE ' Singular matrix in LUDCMP' vv(i)=1.d0/aamax 12 continue do 19 j=1,n do 14 i=1,j-1 sum=a(i,j) do 13 k=1,i-1 sum=sum-a(i,k)*a(k,j) 13 continue a(i,j)=sum 14 continue aamax=0.d0 do 16 i=j,n sum=a(i,j) do 15 k=1,j-1 sum=sum-a(i,k)*a(k,j) 15 continue a(i,j)=sum dum=vv(i)*abs(sum) if (dum.ge.aamax) then imax=i aamax=dum endif 16 continue if (j.ne.imax) then do 17 k=1,n dum=a(imax,k) a(imax,k)=a(j,k) a(j,k)=dum 17 continue d=-d vv(imax)=vv(j) endif indx(j)=imax if (a(j,j).eq.0.d0) a(j,j)=tiny if (j.ne.n) then dum=1.d0/a(j,j) do 18 i=j+1,n a(i,j)=a(i,j)*dum 18 continue endif 19 continue return end c subroutine LUBKSB(a,n,np,indx,b) implicit double precision (a-h,o-z) c c Solves the set of N linear equations A.X=B. Here A is input, not as c the matrix A but rather its LU decomposition, determined by the routine c LUDCMP. INDX is input as the permutation vector returned by LUDCMP. B c is input as teh right hand side vector B and returns with the solution c vector X. A, N, NP, and INDX are not modified by this routine and can c be left in place for successive calls with different right-hand sides B. c This routie takes into account the possibility that B will begin with c many zero elements, so it is efficient for use in matrix inversion. c c Numerical Recipes, p. 39 c INTEGER n,np,indx(n) DOUBLE PRECISION a(np,np),b(n) INTEGER i,ii,j,ll DOUBLE PRECISION sum ii=0 do 12 i=1,n ll=indx(i) sum=b(ll) b(ll)=b(i) if (ii.ne.0) then do 11 j=ii,i-1 sum=sum-a(i,j)*b(j) 11 continue else if (sum.ne.0.d0) then ii=i endif b(i)=sum 12 continue do 14 i=n,1,-1 sum=b(i) do 13 j=i+1,n sum=sum-a(i,j)*b(j) 13 continue b(i)=sum/a(i,i) 14 continue return end c 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 c subroutine OBCS(fll,Teffl,sstop,rl,pl,tl, 1 dlRlL,dlRlTe,dlPlL,dlPlTe,dlTlL,dlTlTe) implicit double precision (a-h,o-z) c c Checks to see if the input values of Fll, Te are within the current c OBC interpolation triangle. If not, it grinds away to find a triangle c that does surround the point. It then computes the conditions at the c base of the envelopes of the triangle, and computes the numerical c logarithmic derivatives of R, P, and T w.r.t. log Te and log L. c Finally, computes the value of logR, logP, and logT for the specified c log(L) and log(Te). c c The triangle has three corners at (cl(i),ct(i)) in the H-R diagram. c At the base of each of the three envelopes, the quanities logP, logT, c and log R are stored in bpl(i), btl(i), and brl(i) respectively c c Common block CTRL contains most of the control parameters for ISUEVO c its used to get dll and dlt common/CTRL/dll,dlt,tcut,dtfac,sstep(7),tstep(7), 1 nmod,nsav,nprint,nzp,itstep c common/KIPP/g1,g2,g3,cl(3),ct(3),brl(3),bpl(3),btl(3),krot common/KFIRST/newrun c if first step of run, compute three new corner positions if (newrun.ne.2) then cl(1)=fll-0.2d0/2.2415927d0 ct(1)=teffl-0.2d0/2.2415927d0 cl(2)=cl(1) ct(2)=ct(1)+dlt cl(3)=cl(1)+dll ct(3)=ct(1) newrun=2 krot=1 endif c icorn1=0 icorn2=0 icorn3=0 c call SETOBC(ICORN1,ICORN2,ICORN3,teffl,fll) c if (icorn1.eq.1) then call ENVELOP(cl(1),ct(1),sstop,rbotl,pbotl,tbotl) write (4,102) 1,rbotl,pbotl,tbotl,ct(1),cl(1) 102 format(' Env #',i2,'; Rbot= ',f7.4,' Pbot= ',f7.4, 1 ' Tbot=',f6.4,' Te=',f6.4,' L= ', f7.4) brl(1)=rbotl bpl(1)=pbotl btl(1)=tbotl endif if (icorn2.eq.1) then call ENVELOP(cl(2),ct(2),sstop,rbotl,pbotl,tbotl) write (4,102) 2,rbotl,pbotl,tbotl,ct(2),cl(2) brl(2)=rbotl bpl(2)=pbotl btl(2)=tbotl endif if (icorn3.eq.1) then call ENVELOP(cl(3),ct(3),sstop,rbotl,pbotl,tbotl) write (4,102) 3,rbotl,pbotl,tbotl,ct(3),cl(3) brl(3)=rbotl bpl(3)=pbotl btl(3)=tbotl endif c Care taken to ensure that the correct differences are taken w.r.t. c triangle parity and orientation. c logarithmic Derivatives with respect to log(L) at constant temperature if (dabs(cl(3)-cl(1)).gt.1d-8) then dlRlL=(brl(3)-brl(1))/(cl(3)-cl(1)) dlPlL=(bpl(3)-bpl(1))/(cl(3)-cl(1)) dlTlL=(btl(3)-btl(1))/(cl(3)-cl(1)) else dlRlL=(brl(2)-brl(1))/(cl(2)-cl(1)) dlPlL=(bpl(2)-bpl(1))/(cl(2)-cl(1)) dlTlL=(btl(2)-btl(1))/(cl(2)-cl(1)) endif c logarithmic Derivatives with respect to log(Te) at consant L if (dabs(ct(3)-ct(1)).gt.1d-8) then dlRlTe=(brl(3)-brl(1))/(ct(3)-ct(1)) dlPlTe=(bpl(3)-bpl(1))/(ct(3)-ct(1)) dlTlTe=(btl(3)-btl(1))/(ct(3)-ct(1)) else dlRlTe=(brl(2)-brl(1))/(ct(2)-ct(1)) dlPlTe=(bpl(2)-bpl(1))/(ct(2)-ct(1)) dlTlTe=(btl(2)-btl(1))/(ct(2)-ct(1)) endif c Interpolation calculation of quantities at base at the point dll=fll-cl(1) dlte=teffl-ct(1) rl=brl(1)+dll*dlRlL+dlte*dlRlTe pl=bpl(1)+dll*dlPlL+dlte*dlPlTe tl=btl(1)+dll*dlTlL+dlte*dlTlTe if (icorn1.eq.1.or.icorn2.eq.1.or.icorn3.eq.1) then write (4,104) dlRlL,dlPlL,dlTlL,dlRlTe,dlPlTe,dlTlTe 104 format (" Derivs: L:(",2(f8.4,","),f8.4,")", 1 " Te: (",2(f8.4,","),f8.4,")") endif c return end c subroutine ENVELOP(fllin,tefflin,sstop,rbotl,pbotl,tbotl) implicit double precision (a-h,o-z) parameter(jmax=1999) c c Envelope integration down to specified mass sstop. Integrates an c envelope at log(L/Lo)=fllin, log(Te)=tefflin c c Common block STRUC contains model structure common/STRUC/h(jmax,30),comp(jmax,15),sr(jmax),nj,modno c Common block GLOB contains global model parameters: xat0 is the initial X common/GLOB/fms,xat0,z,alfa,time,dtime,pcl,tcl,fls,teffl c Common block CONST contains fundamental physical constants COMMON/CONST/CMS,CRS,CLS,CG,CSB,CPI,CP4,CCL,CME,CH,CNA,CEV,CK,CSYR common/TORHS/cmp(15),sp,st,ishell dimension y1(4) COMMON/PATH/KMAX,KOUNT,DXSAV,XP(2000),YP(50,2000) external RHS,rkqs kmax=2000 dxsav=0.d0 c fl=10.d0**fllin teff=10.d0**tefflin x=comp(nj,1)+comp(nj,2) y=comp(nj,4)+comp(nj,3) c12=comp(nj,8) o16=comp(nj,12) c c First the atmospheric integration to obtain the pressure and mass c at the photosphere. call ATMINT(fl,teff,alfa,x,y,z,c12,o16,fms,fmenv,Ro,peff) reff=dsqrt(cls*fl/cp4/csb/teff**4)/crs y1(1)=dlog(reff) y1(2)=dlog(peff) y1(3)=dlog(Teff) y1(4)=fl do 10 i=1,15 cmp(i)=comp(nj,i) 10 continue c c Now integrate downwards to interior mass point with Sr=Sstop c using the RHS subroutine for the derivatives of the physical c quantities. c itmp=ishell ishell=0 sfirst=fmenv/10.d0 eps=3.d-6 call ODEINT(y1,4,fmenv,sstop,eps,sfirst,0.d0,nok,nbad, 1 RHS,rkqs) ishell=itmp rbotl=y1(1)/dlog(10.d0) pbotl=y1(2)/dlog(10.d0) tbotl=y1(3)/dlog(10.d0) c return end c subroutine SETOBC(Icorn1,Icorn2,Icorn3,teffl,fll) implicit double precision (a-h,o-z) c c Determine the new points in the OBC grid within the H-R diagram, and c set the Icorn flag to signal the integrator. Uses values of g1,g2,g3 c computed by IZITOUT to flip triangle once to try and enclose the point c in the triangle. Does NOT reset the Icorns so that multiple calls can c result in multiple resets of corner positions without forgetting what was c done before. c common/KIPP/g1,g2,g3,cl(3),ct(3),brl(3),bpl(3),btl(3),krot c c Calculate new postions in L,Te space for corners of triangle c 100 continue g1=krot*((cl(2)-cl(3))*(Teffl-ct(2)) 1 +(ct(3)-ct(2))*(fll-cl(2))) g2=krot*((cl(3)-cl(1))*(Teffl-ct(3)) 2 +(ct(1)-ct(3))*(fll-cl(3))) g3=krot*((cl(1)-cl(2))*(Teffl-ct(1)) 3 +(ct(2)-ct(1))*(fll-cl(1))) c G1 is less than zero: Compute new position for point 1, change c sense of rotation, and set flag for new calculation of point 1 if (g1.lt.0.d0) then cl(1)=cl(3)+cl(2)-cl(1) ct(1)=ct(3)+ct(2)-ct(1) icorn1=1 krot=-krot go to 100 endif c G2 is less than zero: Compute new position for point 2, change c sense of rotation, and set flag for new calculation of point 2 if (g2.lt.0.d0) then cl(2)=2.d0*cl(1)-cl(2) ct(2)=2.d0*ct(1)-ct(2) icorn2=1 krot=-krot go to 100 endif c G3 is less than zero: Compute new position for point 3, change c sense of rotation, and set flag for new calculation of point 3 if (g3.lt.0.d0) then cl(3)=2.d0*cl(1)-cl(3) ct(3)=2.d0*ct(1)-ct(3) icorn3=1 krot=-krot go to 100 endif c return end c subroutine PTLTE(pl,tl,teobcl,flobcl) implicit double precision (a-h,o-z) parameter(jmax=1999) c c Subroutine to compute the value of log(Te) for a model c from the values of ln(P),ln(T) at the last interior zone c c Common block STRUC contains model structure common/STRUC/h(jmax,30),comp(jmax,15),sr(jmax),nj,modno c Common block GLOB contains global model parameters: xat0 is the initial X common/GLOB/fms,xat0,z,alfa,time,dtime,pcl,tcl,fls,teffl c Common block KIPP contains info about the OBC envelope triangle common/KIPP/g1,g2,g3,cl(3),ct(3),brl(3),bpl(3),btl(3),krot 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 Compute derivatives of P,T w.r.t. L,Te. Make sure that the correct c triangle corners are used in derivatives! c logarithmic Derivatives with respect to log(L) at constant temperature c if (dabs(cl(3)-cl(1)).gt.1d-8) then dlRlL=(brl(3)-brl(1))/(cl(3)-cl(1)) dlPlL=(bpl(3)-bpl(1))/(cl(3)-cl(1)) dlTlL=(btl(3)-btl(1))/(cl(3)-cl(1)) else dlRlL=(brl(2)-brl(1))/(cl(2)-cl(1)) dlPlL=(bpl(2)-bpl(1))/(cl(2)-cl(1)) dlTlL=(btl(2)-btl(1))/(cl(2)-cl(1)) endif c logarithmic Derivatives with respect to log(Te) at consant L if (dabs(ct(3)-ct(1)).gt.1d-8) then dlRlTe=(brl(3)-brl(1))/(ct(3)-ct(1)) dlPlTe=(bpl(3)-bpl(1))/(ct(3)-ct(1)) dlTlTe=(btl(3)-btl(1))/(ct(3)-ct(1)) else dlRlTe=(brl(2)-brl(1))/(ct(2)-ct(1)) dlPlTe=(bpl(2)-bpl(1))/(ct(2)-ct(1)) dlTlTe=(btl(2)-btl(1))/(ct(2)-ct(1)) endif c Derivatve of Te w.r.t. P,T at base dlTelP=1.d0/(dlPlTe-dlPlL*(dlTlTe/dlTlL)) dlTelT=1.d0/(dlTlTe-dlTlL*(dlPlTe/dlPlL)) c Derivatve of L w.r.t. P,T at base dlLlP=1.d0/(dlPlL-dlPlTe*(dlTlL/dlTlTe)) dlLlT=1.d0/(dlTlL-dlTlTe*(dlPlL/dlPlTe)) c Effective temperature teobcl=ct(1)+(tl-btl(1))*dlTelT+(pl-bpl(1))*dlTelP c Luminosity flobcl=cl(1)+(tl-btl(1))*dlLlT+(pl-bpl(1))*dlLlP c return end c c subroutine neutrino(dl,tl,x,y,c,o,z, 1 eplas,ephot,epair,erec,ebrem,etot) implicit double precision (a-h,o-z) double precision mue c c Neutrino rates from Itoh et al. (1996), Ap.J. Supp 102, 411 c (plama, photo, pair, bremmstrahlung) c and Beaudet, Petrosian, and Salpeter (1967), Ap. J. 150, 979. c (recomb.) c c Standardized for ISUEVO by SDK: December 1990 c original BPS updated to Itoh: September 1996 c cln=dlog(10.d0) c t=10.d0**tl t=dexp(cln*tl) c rho=10.d0**dl rho=dexp(cln*dl) mue=2.d0/(1.d0+X) zbar=x+2.d0*y+6.d0*c+8.d0*o if (t.lt.1.d7.or.rho.lt.1.d3) then eplas=0.d0 ephot=0.d0 epair=0.d0 erec=0.d0 ebrem=0.d0 etot=0.d0 return endif c flam=t/5.9302d9 cv=0.5d0+2.0d0*0.232d0 cvp=1.d0-cv ca=0.5d0 cap=1.d0-ca temp=1.018d0*dexp(0.6666667*dlog(rho/mue)) tf=5.9302d9*(dexp(0.5*dlog(1.d0+temp))-1) call PLASMA(flam,rho,mue,t,cv,cvp,ca,cap,qplas) call BREMM(t,rho,mue,x,y,c,o,zbar,cv,cvp,ca,cap,qbrem) call PAIR(flam,rho,mue,cv,cvp,ca,cap,qpair) call PHOTO(flam,rho,t,mue,cv,cvp,ca,cap,qphot) call RECOMB(flam,t,rho,zbar,qrec) eplas=qplas/rho ephot=qphot/rho epair=qpair/rho erec=qrec/rho ebrem=qbrem/rho etot=eplas+ephot+epair+erec+ebrem return end c subroutine PLASMA(flam,rho,mue,t,cv,cvp,ca,cap,qplas) c from Itoh et al. 1996, checks based on plots in that paper and tables c in the CDROM version implicit double precision (a-h,o-z) double precision mue top=1.1095d11*rho/mue temp=dexp(0.6666667d0*dlog(1.019d-6*rho/mue)) bot=t*t*dexp(0.5d0*dlog(1+temp)) gam=dexp(0.5d0*dlog(top/bot)) ft=2.4d0+0.6d0*dexp(0.5d0*dlog(gam))+0.51d0*gam 1 +1.25d0*dexp(1.5d0*dlog(gam)) fl=(8.6d0*gam*gam+1.35d0*dexp(3.5d0*dlog(gam)))/ 1 (225d0-17*gam+gam*gam) x=0.16666667*(17.5+dlog10(2*rho/mue)-3*dlog10(t)) y=0.16666667*(-24.5+dlog10(2*rho/mue)+3*dlog10(t)) n=2 if (x*x.gt.0.49d0.or.y.lt.0) then fxy=1.d0 else fxy=1.05d0+(0.39-1.25*x-0.35d0*sin(4.5d0*x) 1 -0.3d0*dexp(-(4.5d0*x+0.9d0)**2)) 2 *exp(-(min(0.d0,y-1.6d0+1.25*x)/(0.57d0-0.25d0*x))**2) endif qv=3.00d21*flam**9*gam**6*dexp(-gam)*(ft+fl)*fxy qplas=(cv**2+n*cvp**2)*qv return end c subroutine BREMM(t,rho,mue,x,y,c,o,z,cv,cvp,ca,cap,qbrem) c Itoh et al 1996; checked with their figures 6 and 7. For carbon, c doesn't do well at high density regime, but okay elsewhere c requires separate computation routines for He, C,and O implicit double precision (a-h,o-z) c double precision mue c cln=2.302585093d0 t8=t/1.d8 rome=rho/mue temp=1.018d0*dexp(0.6666667*dlog(rho/1d6/mue)) tf=5.9302d9*(dexp(0.5*dlog(1.d0+temp))-1.d0) c if ((y.gt.0.5d0.and.t.gt.0.01d0*tf).or.(t.gt.0.3d0*tf)) then c weakly degenerate electrons a0=23.5d0 a1=6.83d4 a2=7.81d8 a3=230.d0 a4=6.75d5 a5=7.66d9 b1=1.47d0 b2=0.0329 b3=7.75d5*dexp(1.5d0*dlog(t8)) 1 +247d0*dexp(3.85d0*dlog(t8)) b4=4.07d0+0.0240d0*dexp(1.4d0*dlog(t8)) b5=4.57d-5*dexp(-0.11d0*dlog(t8)) c bot1=a3+a4/(t8*t8)+a5/(dexp(5.d0*dlog(t8))) bot1=bot1*(1.d0+1.d-9*rome) bot2=b3/rome+b4+b5*dexp(0.656d0*dlog(rome)) ggas=1.d0/bot1 + 1.d0/bot2 eta=rome/(7.05d6*dexp(1.5d0*dlog(t8))+5.12d4*t8*t8*t8) c bot1=a0+a1/t8/t8+a2*dexp(-5d0*dlog(t8)) bot2=1.d0+b1/eta+b2/eta/eta fgas=1.d0/bot1+1.26d0*(1.d0+1.d0/eta)/bot2 c n=2 c zbar is the mean value (mass-weighted) of the nuclear charge Z zbar=x+2.d0*y+6.d0*c+8.d0*o qgas=0.5738d0*(2.d0*zbar)*dexp(6.d0*dlog(t8))*rho temp1=0.5*((cv*cv+ca*ca)+n*(cvp*cvp+cap*cap))*fgas temp2=0.5*((cv*cv-ca*ca)+n*(cvp*cvp-cap*cap))*ggas qgas=qgas*(temp1-temp2) qbrem=qgas return else c c strongly degenerate electrons (only valid for 1