program teos implicit none save include 'vector_eos.dek' c.. c..tests the iben eos routine c.. c..ionmax = number of isotopes in the network c..xmass = mass fractions c..ymass = molar fractions c..aion = number of nucleons c..zion = number of protons integer ionmax parameter (ionmax=3) double precision xmass(ionmax),ymass(ionmax), 1 aion(ionmax),zion(ionmax),temp,den,abar,zbar c..set the mass fractions, z's and a's of the composition c..hydrogen xmass(1) = 0.75d0 aion(1) = 1.0d0 zion(1) = 1.0d0 c..helium xmass(2) = 0.23d0 aion(2) = 4.0d0 zion(2) = 2.0d0 c..carbon 12 xmass(3) = 0.02d0 aion(3) = 12.0d0 zion(3) = 6.0d0 c..get abar, zbar and a few other composition variables call azbar(xmass,aion,zion,ionmax, 1 ymass,abar,zbar) c..set the input vector temp_row(1) = 1.0d8 den_row(1) = 1.0d6 abar_row(1) = abar zbar_row(1) = zbar c..here the pipeline is only 1 element long jlo_eos = 1 jhi_eos = 1 c..call the eos call eosiben(xmass,aion,zion,ionmax) c..write out the results call pretty_eos_out('iben: ') stop 'normal termination' end subroutine azbar(xmass,aion,zion,ionmax, 1 ymass,abar,zbar) implicit none save c..this routine calculates composition variables for an eos routine c..input: c..mass fractions = xmass(1:ionmax) c..number of nucleons = aion(1:ionmax) c..charge of nucleus = zion(1:ionmax) c..number of isotopes = ionmax c..output: c..molar abundances = ymass(1:ionmax), c..mean number of nucleons = abar c..mean nucleon charge = zbar c..declare integer i,ionmax double precision xmass(ionmax),aion(ionmax),zion(ionmax), 1 ymass(ionmax),abar,zbar,zbarxx,ytot1 zbarxx = 0.0d0 ytot1 = 0.0d0 do i=1,ionmax ymass(i) = xmass(i)/aion(i) ytot1 = ytot1 + ymass(i) zbarxx = zbarxx + zion(i) * ymass(i) enddo abar = 1.0d0/ytot1 zbar = zbarxx * abar return end subroutine pretty_eos_out(whose) implicit none save include 'vector_eos.dek' c.. c..writes a pretty output for the eos tester c.. c..declare integer j character*7 whose c..popular formats 01 format(1x,t2,a,t11,'total',t24,'ion',t34,'electron', 1 t46,'positron',t58,'radiation',t70,'coulomb') 02 format(1x,t2,a,1p6e12.4) 03 format(1x,t2,a6,1pe12.4,t22,a6,1pe12.4, 1 t42,a6,1pe12.4,t62,a6,1pe12.4) c..loop over all elements in the pipe do j=jlo_eos,jhi_eos c..the input write(6,03) 'temp =',temp_row(1),'den =',den_row(1), 1 'abar =',abar_row(1),'zbar =',zbar_row(1) write(6,*) ' ' c..and the output c..first the totals from each of the components write(6,01) whose write(6,02) 'pres =', 1 ptot_row(j),pion_row(j),pele_row(j), 2 ppos_row(j),prad_row(j),pcou_row(j) write(6,02) 'ener =', 1 etot_row(j),eion_row(j),eele_row(j), 2 epos_row(j),erad_row(j),ecou_row(j) write(6,02) 'entr =', 1 stot_row(j),sion_row(j),sele_row(j), 2 spos_row(j),srad_row(j),scou_row(j) c..derivatives of the totals with respect to the input variables write(6,*) ' ' write(6,03) 'dp/dd=',dpd_row(j),'dp/dt=',dpt_row(j), 1 'dp/da=',dpa_row(j),'dp/dz=',dpz_row(j) write(6,03) 'de/dd=',ded_row(j),'de/dt=',det_row(j), 1 'de/da=',dea_row(j),'de/dz=',dez_row(j) write(6,03) 'ds/dd=',dsd_row(j),'ds/dt=',dst_row(j), 1 'ds/da=',dsa_row(j),'ds/dz=',dsz_row(j) c..derivatives of the electron-positron compoenets with c..respect to the input variables write(6,*) ' ' write(6,03) 'dpepd=',dpepd_row(j),'dpept=',dpept_row(j), 1 'dpepa=',dpepa_row(j),'dpepz=',dpepz_row(j) write(6,03) 'deepd=',deepd_row(j),'deept=',deept_row(j), 1 'deepa=',deepa_row(j),'deepz=',deepz_row(j) write(6,03) 'dsepd=',dsepd_row(j),'dsept=',dsept_row(j), 1 'dsepa=',dsepa_row(j),'dsepz=',dsepz_row(j) c..the thermodynamic consistency relations, these should all be c..at the floating poiint limit of zero write(6,*) ' ' write(6,03) 'maxw1=',dse_row(j),'maxw2=',dpe_row(j), 1 'maxw3=',dsp_row(j) c..number density of electrons, poistrons, matter electrons, and ions write(6,03) 'xne =',xne_row(j),'xnp =',xnp_row(j), 1 'xnem =',xnem_row(j),'xni =',xni_row(j) c..derivatibves of the electron number density with c..respect to the input variables write(6,03) 'dxned=',dxned_row(j),'dxnet=',dxnet_row(j), 1 'dxnea=',dxnea_row(j),'dxnez=',dxnez_row(j) c..electron chemical potential, positron chemical potential c..and derivatives of electron chemical potential with respect c..to the input variables write(6,03) 'eta =',etaele_row(j),'etap =',etapos_row(j) write(6,03) 'detad=',detad_row(j),'detat=',detat_row(j), 1 'detaa=',detaa_row(j),'detaz=',detaz_row(j) c..specific heats, and ratio of electostatic to thermal energy write(6,03) 'cp =',cp_row(j),'cv =',cv_row(j), 1 'plasg=',plasg_row(j) c..the 3 gammas and the sound speed write(6,03) 'gam1 =',gam1_row(j),'gam2 =',gam2_row(j), 1 'gam3 =',gam3_row(j),'csond=',cs_row(j) write(6,*) ' ' enddo return end subroutine eosiben(xmass,aion,zion,ionmax) implicit double precision (a-h,o-z) include 'vector_eos.dek' c.. c..front end for the iben eos c.. c..declare integer ionmax double precision xmass(ionmax),aion(ionmax),zion(ionmax), 1 aonly(9),zonly(9) double precision clight,nabad,dse,dpe,dsp,avo,kerg parameter (clight = 2.99792458e10, 1 kerg = 1.380658d-16, 2 avo = 6.0221367d23) common /xyx/xavi(9),dmy(5),dmx(9) common /eosp/x,hell,hell1,zl,zls,xyznuc,xyzele,xyexyt,wtnorm, 1z2oa,z2oah,z2oa13,gliq0,gsol0,deblam,deg1,sphtc1 common /evoc1/xc12,z,z2,z12,z4,xhsurf,xcsurf, 1t,rho,p,pgas,blt,blp,blr,v,vrad,vad,opacr,opacc,entrop,entrpn, 2dly,dlys,pgt1,deg,gamd,eppg,sphtcp,sphtq,vaovap, 3vaovat,cpocpp,cpocpt,epoepp,epoept,drorp,drort,e11,e33,e34, 4e17,ee7,e112,e114,e116,e444,e412,e414,e416,e418,enuc,enucre,ble, 5dledlt,dledlp,dledlr,egrav,epp,ecn,opac,ok,blk,blok,dltp,qstop, 6dlkdlt,dlkdlp,sl,t0,qtop,dqtop,h0,ha0,pgpr,dqbyq,dybyy1, 7dybyy2,dybyy3,twkz,twkz1,rwkz,rwkz1,th2,ta,tb,trapid, 8convc1,convc2,rmixs,alpha1,dtv,dtmpmx,enuch,enuche, 9dtc,dtime,rtime,dlymax,dxmax,enoeg,gamd1,eppg1,cm0,cm1,cp0,cpm, 1 prad,pnuc,erad,eion,eele,entrpe,entrpr,psilnr,psilnt common /evoc2/tpr75,tpr750,cpm1,cpm0,vaovar,cpocpr,epoepr, 1dludlr,dlpdlt,dlpdlr,depdlr,dlkdlr,uint,eatbd,ekt, 2duidlt,duidlr,dlyl,ebind,ekin1,gravsm,ekin2,eneutb,eneut,eneufr, 3epopg,propg,enucr,c12c12,screen,dwork,duint,gamcl,psi, 4dqgdr,dqgdt,dqrdlr,dqrdlt,dqtdlr,dqtdlt,egls,deglst,deglsr,qsdot common /evoc3/cp1,ce0,tpr0,tpr1,dalsl,dalts,blts,blsl,sl0, 1sl1,ts0,ts1,blsl0,blsl1,blts0,blts1,sigho,sighm,sighb,dsigh, 2delt,time,dtmaxi,dtmini,dhecma,dhecmi,dx4dtc,q1min,dmlmi,dmpc, 3dmlc,xcntr,ycntr,xh1,xh2,xhe1,xhe2,dlyms,stime,aclock, 4ep,depdlt,depdlp,elhyd,elhyd0,elhe,elhe0, 5e1212,e1216,convc3,dludlt,dludlp,zneut common /ievoc1/i3,n,kstati,kstate,lskip,kskip,kion,jloss,isc,ienv, 6inuc,np1,ixtype,lrecy,lxtype,jrecy,jt,jl,itmax,it,ihalve,iconv, 7icycm,icycp,icycpp,icyc,iquit,jcon70,irecy,ne,itstp,icom,nfreez, 8lyo,lxo,istart,iguess,nmax,ipassi,ijktp,kx,iclk,iaccrt,iadrad data aonly/1.0d0, 3.0d0, 4.0d0, 12.0d0, 14.0d0, 1 16.0d0,18.0d0, 22.0d0, 25.0d0/ data zonly/1.0d0, 2.0d0, 2.0d0, 6.0d0, 7.0d0, 8.0d0, 8.0d0, 1 10.0d0, 12.0d0/ c..order of the xavi array is h1, he3, he4, c12, n14, o16, o18, ne22, mg25 c..without extensive work to generalize, this is the only composition allowed. c..see the above data statement do i=1,9 xavi(i) = 1.0e-20 enddo do i=1,ionmax if (aion(i) .eq. aonly(1) .and. zion(i) .eq. zonly(1)) 1 xavi(1) = xmass(i) if (aion(i) .eq. aonly(2) .and. zion(i) .eq. zonly(2)) 1 xavi(2) = xmass(i) if (aion(i) .eq. aonly(3) .and. zion(i) .eq. zonly(3)) 1 xavi(3) = xmass(i) if (aion(i) .eq. aonly(4) .and. zion(i) .eq. zonly(4)) 1 xavi(4) = xmass(i) if (aion(i) .eq. aonly(5) .and. zion(i) .eq. zonly(5)) 1 xavi(5) = xmass(i) if (aion(i) .eq. aonly(6) .and. zion(i) .eq. zonly(6)) 1 xavi(6) = xmass(i) if (aion(i) .eq. aonly(7) .and. zion(i) .eq. zonly(7)) 1 xavi(7) = xmass(i) if (aion(i) .eq. aonly(8) .and. zion(i) .eq. zonly(8)) 1 xavi(8) = xmass(i) if (aion(i) .eq. aonly(9) .and. zion(i) .eq. zonly(9)) 1 xavi(9) = xmass(i) enddo c..the least we can do is normalize the composition sum = 0.0d0 do i=1,9 sum = sum + xavi(i) enddo sum = 1.0d0/sum do i=1,9 xavi(i) = xavi(i) * sum enddo c..mystery parameter dltp dltp = 1.0d-3 kstate = 3 c..start the vectorization loop eosfail = .false. do j = jlo_eos,jhi_eos tin = temp_row(j) t = tin * 1.0d-6 rho = den_row(j) c..call the eos call eosprp call dens c..convert the derivatives dpt = dlpdlt/tin*p dpd = dlpdlr*p/rho det = duidlt/tin ded = duidlr/rho dpet = depdlt/tin*ep dped = depdlr*ep/rho c..the temperature and density exponents (c&g 9.81 9.82) c..the specific heat at constant volume (c&g 9.92) c..the third adiabatic exponent (c&g 9.93) c..the first adiabatic exponent (c&g 9.97) c..the second adiabatic exponent (c&g 9.105) c..the specific heat at constant pressure (c&g 9.98) c..and relativistic formula for the sound speed (c&g 14.29) zz = p/rho chit = tin/p * dpt chid = dpd/zz cv = det x = zz * chit/(tin * cv) gam3 = x + 1.0d0 gam1 = chit*x + chid nabad = x/gam1 gam2 = 1.0d0/(1.0d0 - nabad) cp = cv * gam1/chid c cp = sphtcp z = 1.0d0 + (uint + clight*clight)/zz sound = clight * sqrt(gam1/z) c..these are the 3 thermodynamic consistency checks c..each is zero if the consistency is perfect dse = 1.0d20 c dse = tin*dst/det - 1.0d0 dpe = (ded*rho**2 + tin*dpt)/p - 1.0d0 c dsp = -dsd*(rho**2/dpt) - 1.0d0 dsp = 1.0d20 c..load the output ptot_row(j) = p dpt_row(j) = dpt dpd_row(j) = dpd dpa_row(j) = 0.0d0 dpz_row(j) = 0.0d0 etot_row(j) = uint det_row(j) = det ded_row(j) = ded dea_row(j) = 0.0d0 dez_row(j) = 0.0d0 stot_row(j) = entrop * kerg * avo dst_row(j) = 0.0d0 dsd_row(j) = 0.0d0 dsa_row(j) = 0.0d0 dsz_row(j) = 0.0d0 prad_row(j) = prad erad_row(j) = erad srad_row(j) = entrpr * kerg * avo pion_row(j) = pnuc eion_row(j) = eion sion_row(j) = entrpn * kerg * avo xni_row(j) = avo/abar_row(j) * den_row(j) pele_row(j) = ep ppos_row(j) = 0.0d0 dpept_row(j) = dpet dpepd_row(j) = dped dpepa_row(j) = 0.0d0 dpepz_row(j) = 0.0d0 eele_row(j) = eele epos_row(j) = 0.0d0 deept_row(j) = 0.0d0 deepd_row(j) = 0.0d0 deepa_row(j) = 0.0d0 deepz_row(j) = 0.0d0 sele_row(j) = entrpe * kerg * avo spos_row(j) = 0.0d0 dsept_row(j) = 0.0d0 dsepd_row(j) = 0.0d0 dsepa_row(j) = 0.0d0 dsepz_row(j) = 0.0d0 xnem_row(j) = zbar_row(j) * xni_row(j) xne_row(j) = xnem_row(j) dxnet_row(j) = 0.0d0 dxned_row(j) = zbar_row(j)/abar_row(j) * avo dxnea_row(j) = 0.0d0 dxnez_row(j) = 0.0d0 xnp_row(j) = 0.0d0 etaele_row(j) = psi detat_row(j) = 0.0d0 detad_row(j) = 0.0d0 detaa_row(j) = 0.0d0 detaz_row(j) = 0.0d0 etapos_row(j) = 0.0d0 pcou_row(j) = 0.0d0 ecou_row(j) = 0.0d0 scou_row(j) = 0.0d0 dse_row(j) = dse dpe_row(j) = dpe dsp_row(j) = dsp cv_row(j) = cv cp_row(j) = cp gam1_row(j) = gam1 gam2_row(j) = gam2 gam3_row(j) = gam3 cs_row(j) = sound enddo return end c.. c.. c.. c.. c.. subroutine eosprp implicit double precision (a-h,o-z) common /xyx/xavi(9),dmy(5),dmx(9) common /eosp/x,hell,hell1,zl,zls,xyznuc,xyzele,xyexyt,wtnorm, 1z2oa,z2oah,z2oa13,gliq0,gsol0,deblam,deg1,sphtc1 x=xavi(1) hell=xavi(2)+xavi(3) hell1=xavi(2)/2.992622+xavi(3)/3.971539 z2oah=3.*xavi(4)+3.5*xavi(5)+4.*xavi(6)+3.56*xavi(7)+4.54*xavi(8) 1+5.76*xavi(9) z2oa=z2oah+xavi(1)+1.3333333*xavi(2)+xavi(3) z2oah=z2oah/3. z2oa13=x+hell*2.52+15.724*xavi(4)+20.328*xavi(5)+25.4*xavi(6) 1 +24.4*xavi(7)+35.7*xavi(8)+49.3*xavi(9) zoa13=xavi(1)+0.66666667*xavi(2)+0.5*(xavi(3)+xavi(4)+xavi(5) 1 +xavi(6))+0.44444444*xavi(7)+0.454545*xavi(8)+0.48*xavi(9) zoa13=zoa13**0.33333333 z53oa=xavi(1)+1.058267*xavi(2)+0.793700*xavi(3)+1.650964*xavi(4) 1 +1.829653*xavi(5)+2.0*xavi(6)+1.777778*xavi(7)+2.109813*xavi(8) 2 +2.515912*xavi(9) z2oa2=xavi(1)+0.4444444*xavi(2)+0.25*(xavi(3)+xavi(4)+xavi(5) 1 +xavi(6))+0.197531*xavi(7)+0.206612*xavi(8)+0.230400*xavi(9) oneoa=xavi(1)+xavi(2)/3.+xavi(3)/4.+xavi(4)/12.+xavi(5)/14. 1 +xavi(6)/16.+xavi(7)/18.+xavi(8)/22.+xavi(9)/26. gliq0=z53oa*zoa13/oneoa z2oa43=xavi(1)+0.924482*xavi(2)+0.629961*xavi(3)+1.310371*xavi(4) 1 +1.452196*xavi(5)+1.587401*xavi(6)+1.356698*xavi(7) 2 +1.622195*xavi(8)+1.969892*xavi(9) gsol0=z2oa43/oneoa deblam=sqrt(z2oa2) zl=xavi(4)+xavi(5)+xavi(6)+xavi(7)+xavi(8)+xavi(9) xyznuc=xavi(1)+hell1+xavi(4)/11.90686+xavi(5)/13.8944 1 +xavi(6)/15.87079+xavi(7)/17.85945+xavi(8)/21.8207+xavi(9)/26. xyzele=(1.+xavi(1))/2.+xavi(2)/6. xyztot=xyznuc+xyzele wtnorm=1./xyztot xyexyt=xyzele/xyztot deg1=.109754*xyzele sphtc1=(8.31433e+13)/wtnorm return end c.. c.. c.. c.. c.. c.. subroutine dens c******************************************************************** c ***** p=p(rho,t,composition) ***** c******************************************************************** implicit double precision (a-h,o-z) common /junk/gurk(100),iiaa(100) common /eosp/x,hell,hell1,zl,zls,xyznuc,xyzele,xyexyt,wtnorm, 1z2oa,z2oah,z2oa13,gliq0,gsol0,deblam,deg1,sphtc1 common /evoc1/xc12,z,z2,z12,z4,xhsurf,xcsurf, 1t,rho,p,pgas,blt,blp,blr,v,vrad,vad,opacr,opacc,entrop,entrpn, 2dly,dlys,pgt1,deg,gamd,eppg,sphtcp,sphtq,vaovap, 3vaovat,cpocpp,cpocpt,epoepp,epoept,drorp,drort,e11,e33,e34, 4e17,ee7,e112,e114,e116,e444,e412,e414,e416,e418,enuc,enucre,ble, 5dledlt,dledlp,dledlr,egrav,epp,ecn,opac,ok,blk,blok,dltp,qstop, 6dlkdlt,dlkdlp,sl,t0,qtop,dqtop,h0,ha0,pgpr,dqbyq,dybyy1, 7dybyy2,dybyy3,twkz,twkz1,rwkz,rwkz1,th2,ta,tb,trapid, 8convc1,convc2,rmixs,alpha1,dtv,dtmpmx,enuch,enuche, 9dtc,dtime,rtime,dlymax,dxmax,enoeg,gamd1,eppg1,cm0,cm1,cp0,cpm, 1 prad,pnuc,erad,eion,eele,entrpe,entrpr,psilnr,psilnt common /evoc2/tpr75,tpr750,cpm1,cpm0,vaovar,cpocpr,epoepr, 1dludlr,dlpdlt,dlpdlr,depdlr,dlkdlr,uint,eatbd,ekt, 2duidlt,duidlr,dlyl,ebind,ekin1,gravsm,ekin2,eneutb,eneut,eneufr, 3epopg,propg,enucr,c12c12,screen,dwork,duint,gamcl,psi, 4dqgdr,dqgdt,dqrdlr,dqrdlt,dqtdlr,dqtdlt,egls,deglst,deglsr,qsdot common /evoc3/cp1,ce0,tpr0,tpr1,dalsl,dalts,blts,blsl,sl0, 1sl1,ts0,ts1,blsl0,blsl1,blts0,blts1,sigho,sighm,sighb,dsigh, 2delt,time,dtmaxi,dtmini,dhecma,dhecmi,dx4dtc,q1min,dmlmi,dmpc, 3dmlc,xcntr,ycntr,xh1,xh2,xhe1,xhe2,dlyms,stime,aclock, 4ep,depdlt,depdlp,elhyd,elhyd0,elhe,elhe0, 5e1212,e1216,convc3,dludlt,dludlp,zneut common /ievoc1/i3,n,kstati,kstate,lskip,kskip,kion,jloss,isc,ienv, 6inuc,np1,ixtype,lrecy,lxtype,jrecy,jt,jl,itmax,it,ihalve,iconv, 7icycm,icycp,icycpp,icyc,iquit,jcon70,irecy,ne,itstp,icom,nfreez, 8lyo,lxo,istart,iguess,nmax,ipassi,ijktp,kx,iclk,iaccrt,iadrad c******************************************************************** c kstate=1 ,perfect gas partial ionization rad pressure c kstate=2 ,perfect gas completely ionized rad pressure c kstate=3 ,relativistic electron degeneracy allowed, nucleons c are gaseous, and Coulomb interactions are neglected c kstate=4 ,electrons as in kstate=3, but nucleons may be c gas, liquid, or solid c kstate=5 ,similar to kstate=4, but better eos which is c thermodynamically self consistent and continuous c enter with rho,t,x,hell1,deg1,xyexyt,gamd,eppg,sphtc1,kstate, c leave with p,blr,,sphtcp,sphtq,vaovar,vaovat,cpocpr, c cpocpt,epoepr,epoept,drorp,drort,vad,etcetera c******************************************************************** i=icom psi=-100. t2=t*t t4=t2*t2 prad=(2.52116e+9)*t4 pperg=(8.31433e+13)*rho*t entrpr=4.*prad/pperg t32or=(t**1.5)/rho egls=1.5 deglst=0. deglsr=0. gamcl=0. go to(50, 400,1400,1400,1400),kstate c******************************************************************** c *** equation of state in regions of partial ionization c includes h2 molecules and a representative c electron donor of low ionization potential *** c******************************************************************** 50 jwkz=1 eatbd=0. alpha2=.00199*alpha1*z wtpg=(8.31433e+13)*rho*t if(t .le. th2) go to 53 jwkz=2 caa=0. h2bind=0. delta=0. vardyk=0. vardya=0. vardyb=0. alpha3=0. alpha4=0. alpha5=0. cm=1. dvarat=0. dvarbt=0. go to 55 53 th1=(5.04039e-3)/t vardyk=exp(2.302585*(th1*(th1*(-th1*(3.268766e-3)+5.619127e-2)- 14.925164)+12.53351)) vardya=2.302585*th1*(th1*((9.806298e-3)*th1-.1123825)+4.925164) vardyb=2.5-2.302585*th1*th1*(.1123825-(1.96126e-2)*th1) h2bind=vardya-2.5 smxp=wtpg*x/vardyk c******************************************************************** c *** vardya=dln(vardyk)/dln(t)=2.5+b/kt, vardyb=-(d/dt)(b/k), c where b=binding energy of h2 and k=boltzman's const *** c******************************************************************** 55 some=t2*sqrt(t) chi1=.157776/t chi2=.285261/t chi3=.631104/t chim=(8.704082e-2)/t xsa1=2.5+chi1 xsa2=2.5+chi2 xsa3=2.5+chi3 xsam=2.5+chim a1=some*exp(35.94006-xsa1) a2=some*exp(37.32632-xsa2) a3=some*exp(35.94006-xsa3) am=some*exp(36.85632-xsam) alamb=10000000. if(x .gt. 1.e-10)alamb=a1/(wtpg*x) c1=.5*alamb*(sqrt(1.+4./alamb)-1.) if(alamb .gt. 10000.) c1=1. cm=1. if(alpha2 .lt. 1.e-15) go to 57 alamb=am/(wtpg*alpha2) cm=.5*alamb*(sqrt(1.+4./alamb)-1.) if(alamb .gt. 10000.) cm=1. 57 ep1=wtpg*(x*c1+alpha2*cm) ep=ep1 trial=0. kdens=1 iit=1 60 trial=trial+1. if(trial .gt. 100.) go to 130 if(jwkz .eq. 2) go to 75 62 a1ep=1.+a1/ep assist=8.*smxp/(a1ep*a1ep) if(assist .le. .00001) go to 65 ca=a1ep*(sqrt(1.+assist)-1.)/(4.*smxp) go to 70 65 ca=1./a1ep 70 caa=ca*ca*smxp c1=ca*a1/ep go to 80 75 c1=a1/(a1+ep) ca=abs(1.-c1) 80 cm=am/(am+ep) com=ep*ep+a2*(a3+ep) if(com .ge. (1.e-36)) go to 85 c2=0. c3=0. go to 100 85 c2=a2*ep/com c3=a2*a3/com 100 c1c2c3=c1*x+hell1*(c2+c3+c3) alpha3=cm*alpha2 ep=wtpg*(c1c2c3+alpha3) go to(105,110),iit 105 iit=2 epa=ep seea=ep1-epa if(seea .eq. 0.) go to 130 go to(62,75),jwkz 110 seeb=epa-ep if(seeb .eq. 0.) go to 130 if(abs(epa-ep1) .le. .00001*epa) go to 130 if(abs(ep-epa) .le. .00001*epa) go to 130 seea=seea/abs(seea) seeb=seeb/abs(seeb) if(abs(seeb-seea) .le. 1.) go to 125 up=(ep-epa)/(epa-ep1) if(up .gt. .9999) go to 115 epnew=(ep-up*epa)/(1.-up) if(epnew .gt. 0.) go to 120 115 epnew=ep/2. ep1=ep go to 60 120 tiny1=min(epnew,ep1,epa) tiny2=min(epnew,epa,ep) c eptry1=amax1(abs(epnew-ep1),abs(epnew-ep)) c eptry2=amax1(abs(epnew-epa),abs(epnew-ep)) eptry1=max(abs(epnew-ep1),abs(epnew-ep)) eptry2=max(abs(epnew-epa),abs(epnew-ep)) epa=epnew ep=epa if(eptry1 .le. .00001*tiny1 .or. eptry2 .le. .00001*tiny2) 1go to 130 125 ep1=ep iit=1 go to 60 130 c1c2c3=c1*x+hell1*(c2+c3+c3) alpha3=cm*alpha2 wte=c1c2c3+alpha3 wtn=x*(1.-caa)+hell1+alpha2 ep=wtpg*wte pnuc=wtpg*wtn wt=1./(wte+wtn) wte=1./wte wtn=1./wtn pgas=wtpg/wt p1=pgas+prad pgbyp=pgas/p1 prbypg=4.*prad/pgas c******************************************************************** c *****program for specific heats,vad,v,etc***** c******************************************************************** 150 xaa=vardya u22=c2*(1.-c2-c3-c3) u21=c2*(xsa2*(1.-c2-c3)-xsa3*c3) u32=c3*(2.-c2-c3-c3) u31=c3*(xsa2*(1.-c2-c3)+xsa3*(1.-c3)) um2=cm*(1.-cm) um1=xsam*um2 alpha4=alpha2*um1 alpha5=alpha2*um2 pgep=(pgas/ep)-1. caad=caa+caa wtpgep=wte ru111=hell1*(u21+u31+u31) ru222=hell1*(u22+u32+u32) cid=1.+caad*(1.-c1*wtpgep*x) u11c=((1.-c1)*xsa1+caad*(xsa1+vardya+wtpgep* 1(ru111+alpha4)))/cid u12c=((1.-c1)+caad*(2.+wtpgep*(ru222+alpha5)))/cid u11=c1*u11c u12=c1*u12c uaa1=(ca*xsa1-(ca+c1)*u11c)/2. uaa2=(ca-(ca+c1)*u12c)/2. u11x=u11*x u12x=u12*x uaa1x=uaa1*x uaa2x=uaa2*x uep=1./(1.+wt*(pgep*(u12x+ru222+alpha5)+uaa2x)) uet=uep*wt*(pgep*(u11x+ru111+alpha4)+uaa1x) u1t=u11-u12*uet u2t=u21-u22*uet u3t=u31-u32*uet umt=um1-um2*uet uaat=uaa1-uaa2*uet u1tx=u1t*x uaatx=uaat*x sphtcp=(-xaa*uaatx+caa*vardyb*x+xsa1*u1tx+hell1*(xsa2*u2t 1+(xsa2+xsa3)*u3t)+alpha2*xsam*umt)*wt+2.5 uepwt=uep*wt xup=(-xaa*uaa2x+xsa1*u12x+hell1*(xsa2*u22+(xsa2+xsa3)*u32) 1+alpha5*xsam)*uepwt+1. sxut=(-uaatx+u1tx+hell1*(u2t+u3t+u3t)+alpha2*umt)*wt+1. sxup=(-uaa2x+u12x+ru222+alpha5)*uepwt+1. sxut=sxut+prbypg*sxup xupn=xup+prbypg*sxup sphtcp=sphtcp+prbypg*(3.+sxut+xup) sphtq=xupn/rho vad=xupn/(pgbyp*sphtcp) sphtcp=sphtcp/wt caadp1=caad+1. sigma=1.-c1/caadp1 sigma1=c1*sigma if(jwkz .eq. 1) delta=caad/(caadp1-c1) sigmap=sigma1*x epmult=1.+wte*(sigmap+ru222+alpha5) emvard=1.-vardya delem=delta*emvard sigmat=sigmap*(xsa1-delem) rmult=1.-wte*sigmap*delta tmult=1.+wte*(sigmat+ru111+alpha4) depdlr=rmult/epmult depdlt=tmult/epmult u1tt=sigma*(xsa1-delem-depdlt) u2tt=u21-u22*depdlt u3tt=u31-u32*depdlt umtt=um1-um2*depdlt uatt=ca*(u1tt+depdlt-xsa1) u1tt=c1*u1tt uaatt=0. u1r=-(delta+depdlr)*sigma u2r=-u22*depdlr u3r=-u32*depdlr umr=-um2*depdlr uar=ca*(u1r+depdlr) u1r=c1*u1r uaar=0. if(jwkz .eq. 2) go to 155 uaatt=-.5*(u1tt+uatt) uaar=-.5*(u1r+uar) 155 continue erad = 3.*prad/rho eion = 1.5*pgas/rho ekt=(1.5*pgas+3.*prad)/rho eatbd=x*((1.-c1)*chi1+h2bind*caa)+hell1*((1.-c3)*(chi2+chi3) 1-c2*chi2)+alpha2*(1.-cm)*chim pgasn=(8.31433e+13)*t eatbd=eatbd*pgasn uint=ekt-eatbd duidlt=(12.*prad+1.5*(pnuc+ep*depdlt))/rho 1+pgasn*(x*(chi1*u1tt-(1.+h2bind)*uaatt+caa*vardyb) 2+hell1*(chi2*u2tt+(chi2+chi3)*u3tt) 3+alpha2*chim*umtt) duidlr=(-3.*prad+1.5*ep*(depdlr-1.))/rho 1+pgasn*(x*(chi1*u1r-(1.+h2bind)*uaar) 2+hell1*(chi2*u2r+(chi2+chi3)*u3r)+alpha2*chim*umr) dqgdr=(duidlr-p1/rho)/rho dqgdt=duidlt/t go to(160,170,180),kdens c******************************************************************** c ***** compute derivatives ***** c******************************************************************** 160 vadi=vad cpi=sphtcp epi=sphtq uinti=uint duidri=duidlr duidti=duidlt pgasi=pgas rhoi=rho p1i=p1 pgbypi=pgbyp prbypi=prbypg pradi=prad ti=t wti=wt wtpgi=wtpg epl=ep vrdyai=vardya vrdybi=vardyb c10=c1 c20=c2 c30=c3 cmet0=cm caa0=caa ca0=ca xsa10=xsa1 xsa20=xsa2 xsa30=xsa3 xsam0=xsam xaa0=xaa chi10=chi1 chi20=chi2 chi30=chi3 chim0=chim h2bin0=h2bind dqdri=dqgdr dqdti=dqgdt depdri=depdlr depdti=depdlt kdens=2 dc1t=u1tt*dltp dc2t=u2tt*dltp dc3t=u3tt*dltp dcmt=umtt*dltp dcat=uatt*dltp dc1r=u1r*dltp dc2r=u2r*dltp dc3r=u3r*dltp dcmr=umr*dltp dcar=uar*dltp drorp=sxup/pgbyp drort=-sxut if(jwkz .eq. 2) go to 165 dvarat=-2.302585*th1*(th1*((2.9418894e-2)*th1-.224765)+4.925164) dvarbt=2.302585*th1*th1*(.224765-(5.88378e-2)*th1) 165 dlpdlt=(4.*prad+pnuc+ep*depdlt-wtpg*x*uaatt)/p1 dlpdlr=(ep*depdlr+pnuc-wtpg*x*uaar)/p1 c******************************************************************** c *** derivatives wrt density *** c******************************************************************** rho=rho+rho*dltp wtpg=wtpg*(1.+dltp) c1=c10+dc1r c2=c20+dc2r c3=c30+dc3r cm=cmet0+dcmr ca=ca+dcar if(jwkz .eq. 2) ca=1.-c1 caa=(1.-c1-ca)/2. go to 130 170 vaovar=(vad-vadi)/(dltp*vadi) cpocpr=(sphtcp-cpi)/(dltp*cpi) epoepr=(sphtq-epi)/(dltp*epi) sxutr=-(sxut+drort)/(dltp*drort) dludlr=cpocpr+.5*(dlpdlr+1.+sxutr) dqrdlr=(dqgdr-dqdri)/dltp dqtdlr=(dqgdt-dqdti)/dltp c******************************************************************** c *** derivatives wrt temperature *** c******************************************************************** kdens=3 t=t+t*dltp prad=prad*(1.+4.*dltp) rho=rhoi chi1=chi10*(1.-dltp) chi2=chi20*(1.-dltp) chi3=chi30*(1.-dltp) chim=chim0*(1.-dltp) xsa1=2.5+chi1 xsa2=2.5+chi2 xsa3=2.5+chi3 xsam=2.5+chim vardya=vrdyai+dvarat*dltp vardyb=vrdybi+dvarbt*dltp if(jwkz .eq. 1) h2bind=vardya-2.5 c1=c10+dc1t c2=c20+dc2t c3=c30+dc3t cm=cmet0+dcmt ca=ca0+dcat if(jwkz .eq. 2) ca=1.-c1 caa=(1.-c1-ca)/2. go to 130 180 vaovat=(vad-vadi)/(dltp*vadi) cpocpt=(sphtcp-cpi)/(dltp*cpi) epoept=(sphtq-epi)/(dltp*epi) sxutt=-(sxut+drort)/(dltp*drort) dludlt=cpocpt+1.+.5*(dlpdlt+sxutt) dqrdlt=(dqgdr-dqdri)/dltp dqtdlt=(dqgdt-dqdti)/dltp sphtcp=cpi*(8.31433e+13) sphtq=epi uint=uinti duidlr=duidri duidlt=duidti dqgdr=dqdri dqgdt=dqdti depdlt=depdti depdlr=depdri rho=rhoi wt=wti pgas=pgasi ep=epl entrpn=xyznuc*(log(t32or/xyznuc)+15.150977) entrop=entrpn+xyzele*(log(t32or/xyzele)+4.5818588)+entrpr vad=vadi p1=p1i pgbyp=pgbypi prbypg=prbypi prad=pradi t=ti c1=c10 c2=c20 c3=c30 cm=cmet0 ca=ca0 caa=caa0 h2bind=h2bin0 gamcl=(0.227/t)*gsol0*(rho**.33333333) go to 2000 c******************************************************************** c******************************************************************** c *** perfect gas equation of state *** c******************************************************************** c******************************************************************** 400 pgas=rho*t*sphtc1 p1=pgas+prad pgbyp=pgas/p1 prbypg=4.*prad/pgas epn=4.-3.*pgbyp cpn=16.-pgbyp*(12.+1.5*pgbyp) sphtcp=sphtc1*cpn/(pgbyp*pgbyp) drorp=1./pgbyp drort=-epn/pgbyp sphtq=-drort/rho vad=epn/cpn cpd=(32.-12.*pgbyp)/cpn epd=4./epn ep=pgas*xyzele*wtnorm erad = 3.*prad/rho eion = 1.5*pgas/rho ekt=(1.5*pgas+3.*prad)/rho dlpdlr=pgbyp dlpdlt=epn uint=ekt duidlt=(1.5*pgas+12.*prad)/rho duidlr=-3.*prad/rho depdlr=1. depdlt=1. vaovar=pgbyp*pgbyp*(1.-pgbyp)*(12.-4.5*pgbyp)/(epn*cpn) vaovat=-vaovar*3. cpocpr=-(1.-pgbyp)*cpd cpocpt=-3.*cpocpr epoepr=(7.-8./pgbyp)/(rho*sphtq) epoept=12.*(1./pgbyp-1.)/(rho*sphtq) sxutr=-prbypg/(1.+prbypg) dludlr=cpocpr+.5*(dlpdlr+1.+sxutr) dludlt=cpocpt+1.+.5*(dlpdlt-3.*sxutr) rhosq=rho*rho rhot=rho*t dqgdr=-(pgas+4.*prad)/rhosq dqgdt=(1.5*pgas+12.*prad)/rhot dqrdlr=(pgas+8.*prad)/rhosq dqrdlt=-(pgas+16.*prad)/rhosq dqtdlr=-12.*prad/rhot dqtdlt=-3.*dqtdlr gamcl=(0.227/t)*gsol0*(rho**.33333333) entrpn=xyznuc*(log(t32or/xyznuc)+15.150977) entrpe=xyzele*(log(t32or/xyzele)+4.5818588) entrop=entrpn+entrpe+entrpr go to 2000 c******************************************************************** c******************************************************************** c ***** relativistically degenerate electrons ***** c******************************************************************** c******************************************************************** 1400 dega=rho*xyzele/(t*sqrt(t)) ccxx if(i .ge. lskip) go to 400 t1=t*1000000. sphtci=(8.31433e+7)*xyznuc c *** need guess for psi *** eta=log(.12375*dega*(1.+.0437*dega)) if(dega .lt. 10.) go to 1430 degap=(1.e-4)*((rho*xyzele)**.66666667) if(degap .gt. .0001) go to 1410 eta=.30022*(dega**.66666667)*(1.-.25*degap) go to 1420 1410 eta=(6004.4/t)*(sqrt(1.+degap)-1.) 1420 eta=eta*(1.-.825/(eta*eta)) 1430 psi=eta kdens=1 if(psi .lt. 0.) go to 20 if(psi .lt. 3.2) go to 10 ef=1.+psi*psi/4. go to 30 10 ef=.43+psi*(.35+.17*psi) go to 30 20 ahelp=exp(psi-2.) ahelp=ahelp+ahelp ef=(ahelp+ahelp)*(1.-ahelp) 30 psi=ef 1450 pgasn=sphtci*t1 rhoi=rho rheft=rho*xyzele rhook=(1.e-6)*rheft m1=1 trial=1. ptry1=1. delp1=1. psi1=psi 1460 call degen(t1,psi,rhe,rhes,rhet,ep,pes,pet,ue,ues,uet,se) ptry=rhe delp=rheft-rhe delpa=abs(delp) c write(6,111) psi,delpa,rhook if(delpa .lt. rhook) go to 1470 if(trial .gt. 1.5) go to 1461 slope=rhes go to 1464 1461 if(trial .lt. 2.5) go to 1463 1462 if(delp1 .gt. 0. .and. delp2 .gt. 0.) go to 1463 if(delp1 .lt. 0. .and. delp2 .lt. 0.) go to 1463 if(delp .gt. 0. .and. delp2 .gt. 0.) go to 1463 if(delp .lt. 0. .and. delp2 .lt. 0.) go to 1463 psi1=psi2 delp1=delp2 ptry1=ptry2 delta=psi-psi1 1463 slope=(ptry-ptry1)/delta 1464 ptry2=ptry1 psi2=psi1 delp2=delp1 ptry1=ptry psi1=psi delp1=delp delta=delp/slope psi=psi+delta trial=trial+1. c write(6,111) trial,psi,delta,rhe,ep,ue c 111 format(1x,1p6e12.4) if(trial .lt. 51.) go to 1460 1470 m1=2 call degen(t1,psi,rhe,rhes,rhet,ep,pes,pet,ue,ues,uet,se) rhon=rhe/xyzele pnuc=pgasn*rhon if(kstate .ge. 4) go to 1472 pgas=ep+pnuc pd=prad+pgas pvkt=-(xyzele*pes/rhes+pgasn)*rhon*rhon ptkt=(pet-(pes*rhet/rhes))+sphtci*rhon+4.*prad/t1 evkt=-(ues/rhes)*rhe*rhe +3.*prad etkt=(uet-(ues*rhet/rhes))*xyzele+1.5*sphtci+12.*prad/(rhon*t1) erad = 3.*prad/rhon eion = 1.5*pgasn eele = ue * xyzele ekt=ue*xyzele+1.5*pgasn+3.*prad/rhon entrpn=xyznuc*(log(t32or/xyznuc)+15.150977) entrpe=xyzele*se entrop=entrpe+entrpn+entrpr pgls=1.0 dpglsr=0.0 dpglst=0.0 c write(6,*) psi go to 1600 c******************************************************************** c******************************************************************** c *** gas,liquid,solid equation of state for nuclei *** c a version of this was used by iben and tutukov in 1984 c******************************************************************** c******************************************************************** 1472 gamcl=(2.27e+5/t1)*gsol0*(rhon**.33333333) debye=(3.47e+3)*deblam*sqrt(rhon)/t1 plam=2.24*debye plamsq=plam*plam egls=1.5 deglsr=0. deglst=0. pgls=1. dpglsr=0. dpglst=0. if(kstate .eq. 5) go to 1500 twoth=2./3. if(gamcl .le. 1.) go to 1480 if(gamcl .gt. 160.) go to 1475 egl=1.5+1.1*log10(gamcl)+plamsq/12. egls=egl angl=1.1/2.30258509 deglr=(angl/3.+plamsq/12.) deglt=-angl-plamsq/6. deglst=deglt deglsr=deglr go to 1480 1475 as=(0.43-0.15084*plamsq)*plam bs=1.+(0.43-0.05028*plamsq)*plam aobs=3.*as/(bs*bs) dest=aobs desr=-0.5*aobs deglst=dest deglsr=desr es=3./bs egls=es 1480 gamkt=gamcl*pgasn*gurk(16) ucoul=-gamkt pcoul=-gamkt*rhon/3. pgasnn=(pgasn*(egls+deglsr)+(pcoul+pcoul)/rhon)*twoth sphcii=sphtci*(egls+deglst) pgasmm=sphcii*twoth pvkt=-(xyzele*pes/rhes+pgasnn)*rhon*rhon ptkt=(pet-(pes*rhet/rhes))+pgasmm*rhon+4.*prad/t1 evkt=-(ues/rhes)*rhe*rhe+3.*prad-pnuc*deglsr-pcoul etkt=(uet-(ues*rhet/rhes))*xyzele+sphcii+12.*prad/(rhon*t1) erad = 3.*prad/rhon eion = egls*pgasn eele = ue * xyzele ekt=ue*xyzele+egls*pgasn+3.*prad/rhon+ucoul pnuc=twoth*egls*pnuc pgas=ep+pnuc+pcoul pd=prad+pgas entrpn=xyznuc*(log(t32or/xyznuc)+15.150977) entrpe=xyzele*se entrop=entrpn+entrpe+entrpr go to 1600 c********************************************************************* c********************************************************************* c This is an equation of state constructed by c Jim MacDonald in 1987, using: c c Abe (Prog. Theor. Phys., 22, 213, 1959) c Bowers and Salpeter(Phys. Rev., 119, 1180, 1960) c --- Debye-Huckel theory c Slattery et al (Phys. Rev. A, 21, 2087, 1980) c Slattery et al (Phys. Rev. A, 26, 2255, 1982) c --- Monte-Carlo simulation c Iyetomi and Ichimaru (Phys. Rev., 25, 2434, 1982) c --- quantum corrections c Hansen (Phys. Rev. A., 8, 3096, 1973) c Hansen and Mazighi(Phys. Rev. A, 18, 1282, 1978) c --- modern Debye theory c Cohen and Heffer (Phys. Rev., 99, 1128, 1955) c Carr (Phys. Rev., 122, 1437, 1961) c********************************************************************* c********************************************************************* c********************************************************************* c********************************************************************* 1500 gamh=(2.27e+5/t1)*(rhon**0.33333333) gamliq=gliq0*gamh gamsol=gsol0*gamh alnbl=log(plam) cee=1./sqrt(1.+(2.2704177e-3)*plamsq) ajay=(plamsq/36.)*cee*(2.+cee)/(1.+cee) aljp=cee*(2.+cee*cee)*plamsq/36. al2jpp=(cee**5.)*plamsq/12. c********************************************************************* c ***** gas or liquid phase ***** c********************************************************************* gamp5=sqrt(gamliq) gamp25=sqrt(gamp5) gam1p25=gamp25*gamliq gam1p5=gamliq*gamp5 gam2=gamliq*gamliq gam3=gamliq*gam2 alng=log(gamliq) if(gamliq .gt. 1.) go to 1510 hg=0.57735027*gam1p5+gam3*(-0.104584+0.172110*alng 1 -0.033724*gam1p5) ghp=0.8660254*gam1p5+gam3*(-0.141641+0.516333*alng 1 -0.151758*gam1p5) g2hpp=0.4330126*gam1p5+gam3*(0.233051+1.032666*alng 1 -0.531153*gam1p5) go to 1530 1510 if(gamliq .gt. 200.) go to 1520 hg=0.897744*gamliq-3.801720*gamp25+0.758240/gamp25+0.81487*alng 1 +2.584778 ghp=0.897744*gamliq-0.95043*gamp25-0.18956/gamp25+0.81487 g2hpp=0.712823*gamp25+0.23695/gamp25-0.81487 go to 1530 1520 hg=0.895929*gamliq-2.96467*alng-103.5892/gamliq+9.39581 ghp=0.895929*gamliq-2.964673+103.5892/gamliq g2hpp=2.964673-207.1784/gamliq 1530 efliq=3.*alnbl-1.5*alng-1.32351-hg+ajay if(gamsol .gt. 160.) go to 1540 1535 ghp=ghp*gurk(16) aljp=aljp*gurk(16) g2hpp=g2hpp*gurk(16) al2jpp=al2jpp*gurk(16) pgls=1.-ghp/3.+aljp/2. egls=1.5-ghp+aljp sgls=-3.*alnbl+1.5*alng+2.82351+hg-ghp-ajay+aljp gamhlp=ghp+g2hpp almhlp=aljp+al2jpp dpglsr=-gamhlp/9.+almhlp/4. dpglst=gamhlp/3.-almhlp/2. deglsr=-gamhlp/3.+almhlp/2. deglst=gamhlp-almhlp sghlp1=1.5+ghp+ghp+g2hpp sghlp2=3.+aljp+aljp+al2jpp dsglsr=sghlp1/3.-sghlp2/2. dsglst=-sghlp1+sghlp2 gamcl=gamliq eflgsl=efliq go to 1550 c******************************************************************** c ***** solid phase ***** c******************************************************************** 1540 alph=0.5711 alpha=1.-alph alam1=1.0643 alam2=2.9438 ex1=plam/alam1 ex2=plam/alam2 d1=debfn(ex1) d2=debfn(ex2) amc1=exp(ex1) amc2=exp(ex2) uamc1=1./amc1 uamc2=1./amc2 alx1=1.125*ex1+3.*log(1.-uamc1)-d1 alx2=1.125*ex2+3.*log(1.-uamc2)-d2 gam2=gamsol*gamsol efsol=-0.895929*gamsol-1612.5/gam2+alph*alx1+alpha*alx2 c *** the smallest free energy wins *** if(efsol .gt. efliq) go to 1535 eflqsl=efsol xalpx1=1.125*ex1+3.*d1 xalpx2=1.125*ex2+3.*d2 x2lpp1=3.*d1-9.*ex1/(amc1-1.) x2lpp2=3.*d2-9.*ex2/(amc2-1.) uhlp1=(-0.895929*gamsol+3225./gam2)*gurk(16) uhlp2=(alph*xalpx1+alpha*xalpx2)*gurk(16) akx1=xalpx1-alx1 akx2=xalpx2-alx2 pgls=uhlp1/3.+uhlp2/2. egls=uhlp1+uhlp2 shlp=(9675./gam2)*gurk(16) sgls=shlp/2.+alph*akx1+alpha*akx2 uhlp=(0.895929*gamsol+6450./gam2)*gurk(16) alhlp=alph*(xalpx1+x2lpp1)+alpha*(xalpx1+x2lpp2) deglsr=-uhlp/3.+alhlp/2. deglst=uhlp-alhlp dpglsr=-uhlp/9.+alhlp/4. dpglst=uhlp/3.-alhlp/2. axkhlp=alph*(-12.*d1+9.*ex1/(amc1-1.)) 1 +alpha*(-12.*d2+9.*ex2/(amc2-1.)) dsglsr=-shlp/3.+axkhlp/2. dsglst=shlp-axkhlp gamcl=gamsol 1550 pgasnn= pgasn*(pgls+dpglsr) pgasmm=sphtci*(pgls+dpglst) sphcii=sphtci*(egls+deglst) pvkt=-(xyzele*pes/rhes+pgasnn)*rhon*rhon ptkt=(pet-(pes*rhet/rhes))+pgasmm*rhon+4.*prad/t1 evkt=-(ues/rhes)*rhe*rhe+3.*prad-pnuc*deglsr etkt=(uet-(ues*rhet/rhes))*xyzele+sphcii+12.*prad/(rhon*t1) erad = 3.*prad/rhon eion = egls*pgasn eele = ue * xyzele ekt=ue*xyzele+egls*pgasn+3.*prad/rhon pnuc=pgls*pnuc pgas=ep+pnuc pd=prad+pgas entrpe=xyzele*se entrpn=xyznuc*sgls entrop=entrpn+entrpe+entrpr c******************************************************************** c******************************************************************** c egrav=sphtq*(dp/dtime)-sphtcp*(dt/dtime) c******************************************************************** c******************************************************************** c******************************************************************** 1600 sphtq=-(pd+evkt)/pvkt sphtcp=etkt+sphtq*ptkt duidlr=-evkt/rhon duidlt=t1*etkt dqgdr=-(evkt+pd)/(rhon*rhon) dqgdt=etkt*1000000. vad=(pd/t1)*(sphtq/sphtcp) go to(1610,1620,1630),kdens 1610 depdlr=pes*rhe/rhes depdlt=pet-pes*rhet/rhes dlpdlr=(depdlr+pnuc*(1.+dpglsr/pgls))/pd dlpdlt=(t1*depdlt+pnuc*(1.+dpglst/pgls)+4.*prad)/pd depdlr=depdlr/ep depdlt=depdlt*t1/ep vadi=vad epi=sphtq cpi=sphtcp pgasi=pgas pradi=prad ekti=ekt uint=ekt psii=psi psilnr=rhe/rhes psilnt=-t1*(rhet/rhes) duidri=duidlr duidti=duidlt dqdri=dqgdr dqdti=dqgdt ugls=egls duglst=deglst duglsr=deglsr gamsv=gamcl drort=rhon*t1*(ptkt/pvkt) drorp=-rhon*(pd/pvkt) pi=pd ti=t1 pe1=ep entrpi=entrop entrni=entrpn kdens=2 c *** derivatives wrt "psi" *** dltpn=dltp*psi psi=psi+dltpn go to 1470 1620 vaovar=(vad/vadi-1.)/dltpn cpocpr=(sphtcp/cpi-1.)/dltpn epoepr=(sphtq/epi-1.)/dltpn sxutr=(rhon*t1*(ptkt/pvkt)-drort)/(dltpn*drort) dqrdlr=(dqgdr-dqdri)/dltpn dqtdlr=(dqgdt-dqdti)/dltpn c *** derivatives wrt temperature *** kdens=3 pray=1.+dltp t1=t1*pray pgasn=pgasn*pray pray=pray*pray pray=pray*pray prad=prad*pray psi=psii go to 1470 1630 vaovat=(vad/vadi-1.)/dltp cpocpt=(sphtcp/cpi-1.)/dltp epoept=(sphtq/epi-1.)/dltp c write(6,*) rhon,t1,ptkt,pvkt c write(6,*) drort,dltp,drort sxutt=(rhon*t1*(ptkt/pvkt)-drort)/(dltp*drort) dqrdlt=(dqgdr-dqdri)/dltp dqtdlt=(dqgdt-dqdti)/dltp c *** final derivatives wrt temperature and density *** vaovat=vaovat+vaovar*psilnt cpocpt=cpocpt+cpocpr*psilnt epoept=epoept+epoepr*psilnt dqrdlt=dqrdlt+dqrdlr*psilnt dqtdlt=dqtdlt+dqtdlr*psilnt vaovar=vaovar*psilnr cpocpr=cpocpr*psilnr epoepr=epoepr*psilnr dqrdlr=dqrdlr*psilnr dqtdlr=dqtdlr*psilnr sxutt=sxutt+sxutr*psilnt sxutr=sxutr*psilnr dludlr=cpocpr+.5*(dlpdlr+1.+sxutr) dludlt=cpocpt+1.+.5*(dlpdlt+sxutt) sphtcp=cpi*1000000. sphtq=epi t1=ti pgas=pgasi prad=pradi vad=vadi ep=pe1 p1=pi prbypg=4.*prad/pgas xyztot=xyznuc+xyzele duidlr=duidri duidlt=duidti dqgdr=dqdri dqgdt=dqdti egls=ugls deglst=duglst deglsr=duglsr gamcl=gamsv entrop=entrpi entrpn=entrni 2000 continue blr=log(rho) c ep=ep*(1.e-17) c p=p1*(1.e-17) p=p1 propg=prad/pgas c pgas=pgas*(1.e-17) c sphtq=sphtq*(1.e+17) epopg=ep/pgas return end c.. c.. c.. c.. c.. function debfn(x) implicit double precision (a-h,o-z) c c calculates debye function for n=3. c parameter (c1=1.0/20.0 ,c2=-1.0/1680.0 , 1 c3=1.0/90720.0,c4=-1.0/4435200.0, 2 c5=1.0/207567360.0) if (x.le.2.75) then c c for small x use first seven terms of series expansion. c x2=x**2 debfn=1.0-0.375*x+x2*(c1+x2*(c2+x2*(c3+x2*(c4+x2*c5)))) else c c for large x use first four terms of asymptotic expansion. c e1=exp(-x) x2=x**2 x3=x*x2 t1=e1*(x3+3.0*x2+6.0*x+6.0) e2=e1**2 t2=e2*(0.5*x3+0.75*x2+0.75*x+0.375) e3=e1*e2 t3=e3*(x3/3.0+x2/3.0+x/4.5+1.0/13.5) debfn=3.0*(6.493939404-(t1+t2+t3))/x3 endif return end c.. c.. c.. c.. c.. subroutine degen(t,ef,rhe,rhes,rhet,ep,pes,pet,ue,ues,uet, 1se) c******************************************************************** c******************************************************************** c ***** private communication from flannery, faulkner,and c eggleton in the spring of 1972 ***** c******************************************************************** c******************************************************************** implicit double precision (a-h,o-z) dimension efd(4,3),efe(5),ge(5) dimension c(64) data n,nm3,nm,enm5h/4,1,16,1.5/ data a1,a2,ct,ge(1),efe(1)/2.92196e+6,1.44069e+24, 11.68629e-10,1.0,1.0/ data c/ 2.315456, 6.748456, 6.565124, 2.132124, 7.837772, * 21.439348, 19.079812, 5.478236, 9.215592, 23.550788, 19.015416, * 4.680220, 3.693276, 8.859896, 6.500728, 1.334108, 1.157728, * 3.770052, 4.015890, 1.403566, 8.283474, 26.184274, 28.210176, * 10.309376, 14.755506, 45.032608, 46.907996, 16.630894, 7.386552, * 22.159656, 22.438128, 7.665024, 2.315456, 7.129408, 7.504622, * 2.665155, 7.837772, 23.507117, 23.311140, 7.987820, 9.215592, * 26.832553, 25.083223, 8.020911, 3.693276, 10.333240, 9.168944, * 2.668216, 3.473184, 10.122684, 9.847686, 3.198186, 16.121246, * 43.475686, 37.851844, 10.497404, 23.971098, 60.391420, 47.781924, * 11.361602, 11.079828, 26.579688, 19.502184, 4.002324/ efp1=ef+1. efofp1=ef/efp1 rtfp1=sqrt(efp1) rtfp12=rtfp1+rtfp1 g=ct*rtfp1*t gplus1=g+1. gogp1=g/gplus1 efdk=sqrt(gogp1)*ef*g/((efp1**n)*(gplus1**nm3)) a2gfdk=a2*g*efdk do 1 m=2,n efe(m)=efe(m-1)*ef 1 ge(m)=ge(m-1)*g do 2 m=1,4 do 2 mm=1,3 2 efd(m,mm)=0. ls=0 do 3 j=1,n jm1=j-1 do 3 k=1,n gjfk=ge(j)*efe(k) km1=k-1 ls=ls+1 do 3 m=1,4 l=nm*(m-1)+ls add=c(l)*gjfk efd(m,1)=efd(m,1)+add efd(m,2)=efd(m,2)+jm1*add 3 efd(m,3)=efd(m,3)+km1*add b1=efd(2,1)/(efd(3,1)*rtfp1) c3=1.5-enm5h*gogp1 c14=c3+1. d1=.5*t/efp1 d2=1.-n*efofp1 d3=c14+efd(4,2)/efd(4,1) d4=efd(3,2)/efd(3,1) d5=efd(3,3)/efd(3,1) ep=a2gfdk*efd(1,1) pet=ep*(c14+efd(1,2)/efd(1,1))/t pes=(ep*(d2+efd(1,3)/efd(1,1))/ef)+pet*d1 rhe=a1*efdk*efd(3,1) rhet=rhe*(c3+d4)/t rhes=rhe*(d2+d5)/ef +rhet*d1 ue=a2gfdk*efd(4,1)/rhe uet=ue*(d3/t-rhet/rhe) ues=ue*((d3*.5/efp1)-(rhes/rhe)+(d2+efd(4,3)/efd(4,1))/ef) se=efd(2,1)/(efd(3,1)*rtfp1)-log(ef/(1.+rtfp1)**2) c *** se=mu(electron)*(m(hydrogen)/k(boltzman))*(electron entropy)*** return end