program teos implicit none save include 'vector_eos.dek' c.. c..tests the arnett 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 eosarnett c..write out the results call pretty_eos_out('arnett:') 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 eosarnett implicit none save include 'arnett_eos.dek' include 'vector_eos.dek' c.. c..top level interface routine for arnett's equation of state c.. c..declare integer kin,ktop,j double precision chit,chid,cv,cp,gam1,gam2,gam3,nabad,sound, 1 bdum,zz,z,dpd,dped,ded,x,afac,dse,dpe,dsp, 2 c,con2,me,mecc,kerg,avo parameter (c = 2.99792458d10, con2=c*c, 2 me = 9.1093897d-28, 3 mecc = me * con2, 4 kerg = 1.380658d-16, 5 avo = 6.0221367d23) c..set the vector input eosfail = .false. kin = jlo_eos ktop = jhi_eos do j=jlo_eos,jhi_eos tem(j) = temp_row(j) rho(j) = den_row(j) ye(j) = zbar_row(j)/abar_row(j) enddo c..do it call state(kin,ktop,abar_row) c..set the vector output do j=jlo_eos,jhi_eos afac = -1.0d0/den_row(j)**2 dpd = afac * pv(j) dped = afac * pev(j) ded = afac * ev(j) 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(j)/rho(j) chit = tem(j)/p(j) * pt(j) chid = dpd/zz cv = et(j) x = zz * chit/(tem(j) * cv) gam3 = x + 1.0d0 gam1 = chit*x + chid nabad = x/gam1 gam2 = 1.0d0/(1.0d0 - nabad) cp = cv * gam1/chid z = 1.0d0 + (e(j) + con2)/zz sound = c * 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 = tem(j)*st(j)/et(j) - 1.0d0 dpe = (ded*rho(j)**2 + tem(j)*pt(j))/p(j) - 1.0d0 c dsp = -sd(j)*(rho(j)**2/pt(j)) - 1.0d0 dsp = 1.0d20 c..load the output ptot_row(j) = p(j) dpt_row(j) = pt(j) dpd_row(j) = dpd dpa_row(j) = 0.0d0 dpz_row(j) = 0.0d0 etot_row(j) = e(j) det_row(j) = et(j) ded_row(j) = ded dea_row(j) = 0.0d0 dez_row(j) = 0.0d0 stot_row(j) = 0.0d0 dst_row(j) = 0.0d0 dsd_row(j) = 0.0d0 dsa_row(j) = 0.0d0 dsz_row(j) = 0.0d0 prad_row(j) = prad(j) erad_row(j) = erad(j) srad_row(j) = 0.0d0 pion_row(j) = pgas(j) eion_row(j) = egas(j) sion_row(j) = 0.0d0 xni_row(j) = avo/abar_row(j) * den_row(j) pele_row(j) = pe(j) ppos_row(j) = 0.0d0 dpept_row(j) = pet(j) dpepd_row(j) = dped dpepa_row(j) = 0.0d0 dpepz_row(j) = 0.0d0 eele_row(j) = ee(j) 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) = 0.0d0 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) = ye(j) * avo * den_row(j) xne_row(j) = xnem_row(j) dxnet_row(j) = 0.0d0 dxned_row(j) = ye(j) * avo dxnea_row(j) = 0.0d0 dxnez_row(j) = 0.0d0 xnp_row(j) = 0.0d0 etaele_row(j) = eta(j) detat_row(j) = etabt(j) detad_row(j) = etabd(j) * ye(j) detaa_row(j) = 0.0d0 detaz_row(j) = 0.0d0 etapos_row(j) = 0.0d0 bdum = kerg * temp_row(j)/mecc if (bdum .gt. 0.02) etapos_row(j) = -eta(j) - 2.0d0/bdum 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 end do return end subroutine state(kin,ktop,abar_row) implicit none save include 'arnett_eos.dek' c.. c..equation of state c..includes ideal ion gas, radiation pressure, and noninteracting c..ferm-dirac gas of electrons and positrons c.. c..declare integer kin,ktop,ieos,k double precision 1 abar_row(1),zero,half,one,three,four,third,cmev, 2 denye,t1,v1,pegas,pegst,pegsv,eegas,eegst,eegsv, 3 sen,si,fact1,fact2,ynuc,t2,t3,t4,pgst, 4 pgsv,prdt,prdv,egst,egsv,erdt,erdv, 5 arad,rgas,cergs,coul1,coul2,xbhc c c 1 rgas = 8.3145112d+7, c parameter (arad = 7.56590624067000303d-15, parameter (arad = 7.5641000E-15, 1 rgas = 8.3170000E+07, 2 cergs = 9.648d17) c.. data ieos/0/ data zero/0.0d0/,half/0.5d0/,one/1.0d0/,three/3.0d0/, 1 four/4.0d0/ c---------------------------------------------------------- c energy constant in cgs cmev = cergs third = one/three c..ieos flag reset in ndeo to 1 if ( ieos .eq. 0 ) call ndeo(ieos) c..electron component c..fills all k values, even if outside table; these will be written over call neos(kin,ktop) do k = kin, ktop denye = ye(k)*rho(k) if( denye .lt. dd(2) )then t1 = tem(k) v1 = 1.0d0/rho(k) call nrfd(denye,t1,pegas,pegst,pegsv,eegas,eegst,eegsv, 1 sen,si,k) pegsv = pegsv/v1 eegas = eegas*v1 eegst = eegst*v1 sen = sen*ye(k) c..other variables already converted if( denye .gt. dd(1) .and. denye .lt. dd(2) )then c..morph into nonrelativistic limit fact1 = ( denye - dd(1) )/( dd(2) - dd(1)) fact1 = max( fact1, 0.0d0 ) fact1 = min( fact1, 1.0d0 ) fact2 = 1.0d0 - fact1 pe(k) = pegas*fact2 + pe(k) *fact1 pet(k) = pegst*fact2 + pet(k)*fact1 pev(k) = pegsv*fact2 + pev(k)*fact1 ee(k) = eegas*fact2 + ee(k) *fact1 eet(k) = eegst*fact2 + eet(k)*fact1 eev(k) = eegsv*fact2 + eev(k)*fact1 elseif( denye .lt. dd(1) )then c..nonrelativistic limit fact1 = 0 pe(k) = pegas pet(k) = pegst pev(k) = pegsv ee(k) = eegas eet(k) = eegst eev(k) = eegsv endif endif enddo c------------------------------------------------------------------- c..ions and radiation component do k = kin, ktop ynuc = 1.0d0/abar_row(k) t1 = tem(k) t2 = t1*t1 t3 = t2*t1 t4 = t2*t2 v1 = 1.0d0/rho(k) denye = rho(k)*ye(k) pgas(k) = rgas*t1/v1 *ynuc pgst = rgas/v1 *ynuc pgsv = -pgas(k)/v1 prad(k) = arad*t4/three prdt = four*arad*t3/three prdv = zero egas(k) = rgas*t1*1.5d0 *ynuc egst = rgas* 1.5d0 *ynuc egsv = zero erad(k) = three*prad(k)*v1 erdt = three*prdt*v1 erdv = three*prad(k) c..add coulomb effects to ionic eos; factors are (z**(5/3) - z), c..which are zero for z=0 and 1 c..boo, this assumes a certain composition order, comment out c coul1 = denye/1.818d14 c coul2 = coul1**third c xbhc = coul2 * ( x( 4,k)* 1.175d0 c 1 + x( 5,k)* 13.81d0 + x( 6,k)* 24.00d0 c 2 + x( 7,k)* 36.41d0 + x( 8,k)* 50.90d0 c 3 + x( 9,k)* 77.32d0 + x(10,k)*230.18d0 c 4 + x(11,k)*216.0d0 + x(12,k)*202.19d0 ) c xbhc = coul2 * 13.81d0 c pcoul(k) = -cmev * xbhc * rho(k)/3.0d0 c ecoul(k) = -cmev * xbhc c ecv(k) = - ecoul(k)/3.00d0 * rho(k) c pcv(k) = - pcoul(k)/0.75d0 * rho(k) c write(6,111) pcoul(k),ecoul(k) c write(6,111) pcv(k),ecv(k) c 111 format(1x,1p2e14.6) pcoul(k) = 0.0d0 ecoul(k) = 0.0d0 ecv(k) = 0.0d0 pcv(k) = 0.0d0 c sum all terms c nuclear terms are not included here, only ideal gas of nuclei p(k) = pgas(k) + prad(k) + pe(k) + pcoul(k) + pdeg(k) pt(k) = pgst + prdt + pet(k) pv(k) = pgsv + prdv + pev(k) + pcv(k) + pdv(k) e(k) = egas(k) + erad(k) + ee(k) + ecoul(k) + edeg(k) et(k) = egst + erdt + eet(k) ev(k) = egsv + erdv + eev(k) + ecv(k) + edv(k) enddo return end subroutine ndeo(ieos) implicit none save include 'arnett_eos.dek' c.. c..reads Arnett tables for equation of state c.. c..declare integer ieos,k,n,n1,n2,n3,n4 double precision d12,d23,d34,d41,d24,d13,d1,d2,d3,d4 c..ieos positive prevents second call to deo ieos = 1 c write(*,*)'Tabular EOS for electrons, with degeneracy and pairs' c..uses unit 14 for extreme relativistic fermi-dirac data c write(*,*)'using eos.d ER EOS table: ',nerfd open(14,file='eos.d',status='old') rewind 14 read(14,11) ( asi(k),aj2(k),aj3(k),aphi(k), k=1, nerfd) close(14) c..uses unit 14 for nonrelativistic fermi-dirac data c write(*,*)'using eos1.d NR EOS table: ',nerfd open(14,file='eos1.d',status='old') rewind 14 read(14,11) ( bsi(k),bj2(k),bj3(k),bphi(k), k=1, nerfd) close(14) c..uses unit 14 for electron-pair data c write(*,*)'using eos7.d EOS table: ',ndd,ntt c open(14,file='eos7.d',status='old') open(14,file='eos8.d',status='old') do 102 n=1, ndd read(14,940) dd(n) read(14,940) ( tt(k), array(n,k,1), array(n,k,2), 1 array(n,k,3), k=1, ntt) 102 continue close(14) c write(*,*)'All EOS tables read' c..convert to cgs units (T9 --> K, P(d-less) --> dyne/cm**2) do 103 k = 1, ntt tt(k) = 1.0d9*tt(k) 103 continue c..pcons = 1.43991d24 in cgs (dynes/cm**2=erg/cm**3) do 105 k = 1, ntt do 104 n = 1, ndd array(n,k,1) = array(n,k,1)*1.43991d24 array(n,k,2) = array(n,k,2)*1.43991d24 104 continue 105 continue c..construct temperature differences do 200 k = 1, ntt if( k .le. 2 )then n1 = 1 n2 = 2 n3 = 3 n4 = 4 elseif( k .eq. ntt )then n4 = ntt n3 = n4-1 n2 = n4-2 n1 = n4-3 else n1 = k-2 n2 = k-1 n3 = k n4 = k+1 endif d12 = tt(n1)-tt(n2) d23 = tt(n2)-tt(n3) d34 = tt(n3)-tt(n4) d41 = tt(n4)-tt(n1) d24 = tt(n2)-tt(n4) d13 = tt(n1)-tt(n3) c d1 = - d12 * d13 * d41 d2 = - d12 * d23 * d24 d3 = d13 * d23 * d34 d4 = d41 * d24 * d34 ttd(1,k) = d1 ttd(2,k) = d2 ttd(3,k) = d3 ttd(4,k) = d4 200 continue c..construct density differences once and for all (used in eos) do 210 k = 1, ndd if( k .le. 2 )then n1 = 1 n2 = 2 n3 = 3 n4 = 4 elseif( k .eq. ndd )then n4 = ndd n3 = n4-1 n2 = n4-2 n1 = n4-3 else n1 = k-2 n2 = k-1 n3 = k n4 = k+1 endif d12 = dd(n1)-dd(n2) d23 = dd(n2)-dd(n3) d34 = dd(n3)-dd(n4) d41 = dd(n4)-dd(n1) d24 = dd(n2)-dd(n4) d13 = dd(n1)-dd(n3) d1 = - d12 * d13 * d41 d2 = - d12 * d23 * d24 d3 = d13 * d23 * d34 d4 = d41 * d24 * d34 ddd(1,k) = d1 ddd(2,k) = d2 ddd(3,k) = d3 ddd(4,k) = d4 210 continue 11 format(f5.1,f12.5,f12.5,f12.5) 940 format(1p4e17.8) return end subroutine neos(kin,kk) implicit none save include 'arnett_eos.dek' c.. c..equation of state: fermi-dirac gas of electrons and positrons c.. c..declare integer kin,kk,k,ii,kd,maxd,mind,kkt,maxt,mint,n1,n2,n3,n4,jd double precision ff(4),f(3),df(3),g(4,3),dg(4,3), 1 one,two,three,third,eight,pcdeg, 2 gamdeg,del1,del2,del3, 3 del4,del12,del23,del24,del34,del41,del13,b1,b2,b3, 4 b4,c1,c2,c3,c4,alf1,alf2,alf3,alf4,bet1,bet2,bet3, 5 bet4,pee,ped,peet,erho,erhod,erhot,etaele,detadd,detadt c.. c.. one = 1 two = 2 three = 3 third = one/three eight = 8 pcdeg = 1.43991d24 / 24.0d0 if( kin .lt. 1 .or. kk .ge. kdm-1 )then write(*,*)'NEOS error: kin,kk ',kin,kk stop'neos' endif do k=kin,kk if( tem(k) .le. 0.0d0 )then write(*,*)'NEOS error: tem ',k,tem(k) stop'neos' endif if( rho(k) .le. 0.0d0 )then write(*,*)'NEOS error: rho ',k,rho(k) stop'neos' endif if( ye(k) .le. 0.0d0 )then write(*,*)'NEOS error: ye ',k,ye(k) stop'neos' endif enddo c--Density and Temperature given do 50 ii = kin, kk d(ii) = rho(ii)*ye(ii) t(ii) = tem(ii) c..degenerate electrons, limiting case, analtyic c x = ( d(ii) / 0.9735d6 )**third c x2 = x * x c x3 = x2 * x c w = sqrt( one + x2 ) c..chandrasekhar functions c if( x .gt. 0.1d0 )then c fx = x*( two*x2 - three)*w c 1 + three * log( w + x ) c gx = eight * x3 *( w - one ) - fx c gamdeg = x *( three*w*( two * x2 - one ) c 1 + ( x2*(two*x2 - three) + three )/w c 2 ) / ( three * fx ) c else c..low x limit.................................... c x5 = x2*x3 c x7 = x5*x2 c x9 = x7*x2 c fx = 1.6d0*x5 - 4.0d0/7.0d0*x7 + 1.0d0/3.0d0*x9 c gx = 2.4d0*x5 - 3.0d0/7.0d0*x7 + 1.0d0/6.0d0*x9 c gamdeg = 5.0d0/3.0d0 c endif c pdeg(ii) = pcdeg * fx c edeg(ii) = pcdeg * gx / rho(ii) pdeg(ii) = 0.0d0 edeg(ii) = 0.0d0 c..thermal part: from tables..................................... c..get indices c..!!!!!! WARNING: silently extrapolates out of table !!!!!!!!!! c..density................................................... if( d(ii) .lt. 0.0d0 )then stop'neos: d < 0' endif kd = ( log10( d(ii) ) + 4.0 )*8.0 + 2.0 c..8.0=1/.125 is log10 increment in density per table element; c..4.0 is log10 of first d value in table; order is 3 (cubic fit) if( kd .ge. ndd )then maxd = ndd mind = maxd - 3 kd = maxd elseif( kd .le. 2 )then mind = 1 maxd = mind + 3 kd = mind else maxd = kd + 1 mind = maxd - 3 endif c..maxd = k + nd/2 and mind = maxd - nd for order nd c..temperature................................................ if( t(ii) .lt. 0.0d0 )then write(*,*)t(ii),ii,d(ii),kd stop'neos: t < 0' endif c..for eos7.d c kkt = ( log10( t(ii) ) - 7.0 )*20.0 + 2.0 c.. c..for eos8.d kkt = ( log10( t(ii) ) - 5.0 )*20.0 + 2.0 c..20.0=1/.05 is log10 increment in temperature per table element; c..7.0 is log10 of first t value in table; order is 3 (cubic fit) if( kkt .ge. ntt )then maxt = ntt mint = maxt - 3 kkt = maxt elseif( kkt .le. 2 )then mint = 1 maxt = mint + 3 kkt = mint else maxt = kkt + 1 mint = maxt - 3 endif c..maxt = k + nt/2 and mint = maxt - nt for order nt c-------------------------------------------------------------- c..interpolate in temperature n1 = mint n2 = mint + 1 n3 = mint + 2 n4 = maxt del1 = t(ii) - tt(n1) del2 = t(ii) - tt(n2) del3 = t(ii) - tt(n3) del4 = t(ii) - tt(n4) del1 = t(ii) - tt(n1) del2 = t(ii) - tt(n2) del3 = t(ii) - tt(n3) del4 = t(ii) - tt(n4) del12 = del1*del2 del23 = del2*del3 del24 = del2*del4 del34 = del3*del4 del41 = del4*del1 del13 = del1*del3 b1 = del2 * del34 b2 = del1 * del34 b3 = del12 * del4 b4 = del12 * del3 c1 = del23 + del24 + del34 c2 = del13 + del41 + del34 c3 = del12 + del41 + del24 c4 = del12 + del13 + del23 alf1 = b1/ttd(1,kkt) alf2 = b2/ttd(2,kkt) alf3 = b3/ttd(3,kkt) alf4 = b4/ttd(4,kkt) bet1 = c1/ttd(1,kkt) bet2 = c2/ttd(2,kkt) bet3 = c3/ttd(3,kkt) bet4 = c4/ttd(4,kkt) c..loop over densities...................................... do 300 jd = mind, maxd c..set up variables (P=1,rho*E=2,mu/kt=3) ff(1) = array(jd,n1,1) ff(2) = array(jd,n2,1) ff(3) = array(jd,n3,1) ff(4) = array(jd,n4,1) f(1) = alf1*ff(1) + alf2*ff(2) + alf3*ff(3) + alf4*ff(4) df(1) = bet1*ff(1) + bet2*ff(2) + bet3*ff(3) + bet4*ff(4) ff(1) = array(jd,n1,2) ff(2) = array(jd,n2,2) ff(3) = array(jd,n3,2) ff(4) = array(jd,n4,2) f(2) = alf1*ff(1) + alf2*ff(2) + alf3*ff(3) + alf4*ff(4) df(2) = bet1*ff(1) + bet2*ff(2) + bet3*ff(3) + bet4*ff(4) ff(1) = array(jd,n1,3) ff(2) = array(jd,n2,3) ff(3) = array(jd,n3,3) ff(4) = array(jd,n4,3) f(3) = alf1*ff(1) + alf2*ff(2) + alf3*ff(3) + alf4*ff(4) df(3) = bet1*ff(1) + bet2*ff(2) + bet3*ff(3) + bet4*ff(4) c..save values at density points for density interpolation g (jd-mind+1,1) = f (1) dg(jd-mind+1,1) = df(1) g (jd-mind+1,2) = f (2) dg(jd-mind+1,2) = df(2) g (jd-mind+1,3) = f (3) dg(jd-mind+1,3) = df(3) 300 continue c............................................................. c..interpolate in density n1 = mind n2 = mind + 1 n3 = mind + 2 n4 = maxd del1 = d(ii) - dd(n1) del2 = d(ii) - dd(n2) del3 = d(ii) - dd(n3) del4 = d(ii) - dd(n4) del12 = del1*del2 del23 = del2*del3 del24 = del2*del4 del34 = del3*del4 del41 = del4*del1 del13 = del1*del3 b1 = del2 * del34 b2 = del1 * del34 b3 = del12 * del4 b4 = del12 * del3 c1 = del23 + del24 + del34 c2 = del13 + del41 + del34 c3 = del12 + del41 + del24 c4 = del12 + del13 + del23 alf1 = b1/ddd(1,kd) alf2 = b2/ddd(2,kd) alf3 = b3/ddd(3,kd) alf4 = b4/ddd(4,kd) bet1 = c1/ddd(1,kd) bet2 = c2/ddd(2,kd) bet3 = c3/ddd(3,kd) bet4 = c4/ddd(4,kd) c..set up variables (P=1,rho*E=2,mu/kt=3) c..Pressure ff(1) = g(1,1) ff(2) = g(2,1) ff(3) = g(3,1) ff(4) = g(4,1) f(1) = alf1*ff(1) + alf2*ff(2) + alf3*ff(3) + alf4*ff(4) pee = f(1) c..dP/d(rho) df(1) = bet1*ff(1) + bet2*ff(2) + bet3*ff(3) + bet4*ff(4) ped = df(1) c..dP/dT ff(1) = dg(1,1) ff(2) = dg(2,1) ff(3) = dg(3,1) ff(4) = dg(4,1) f(1) = alf1*ff(1) + alf2*ff(2) + alf3*ff(3) + alf4*ff(4) peet = f(1) c..rho*E ff(1) = g(1,2) ff(2) = g(2,2) ff(3) = g(3,2) ff(4) = g(4,2) f(2) = alf1*ff(1) + alf2*ff(2) + alf3*ff(3) + alf4*ff(4) erho = f(2) c..d(rho*E)/d(rho) df(2) = bet1*ff(1) + bet2*ff(2) + bet3*ff(3) + bet4*ff(4) erhod = df(2) c..d(rho*E)/dT ff(1) = dg(1,2) ff(2) = dg(2,2) ff(3) = dg(3,2) ff(4) = dg(4,2) f(2) = alf1*ff(1) + alf2*ff(2) + alf3*ff(3) + alf4*ff(4) erhot = f(2) c..mu/kT ff(1) = g(1,3) ff(2) = g(2,3) ff(3) = g(3,3) ff(4) = g(4,3) f(3) = alf1*ff(1) + alf2*ff(2) + alf3*ff(3) + alf4*ff(4) etaele = f(3) c..d(eta)/d(rho) df(3) = bet1*ff(1) + bet2*ff(2) + bet3*ff(3) + bet4*ff(4) detadd = df(3) c..d(eta)/dT ff(1) = dg(1,3) ff(2) = dg(2,3) ff(3) = dg(3,3) ff(4) = dg(4,3) f(3) = alf1*ff(1) + alf2*ff(2) + alf3*ff(3) + alf4*ff(4) detadt = f(3) c------------------------------------------------------------------ pe(ii) = pee pet(ii) = peet pev(ii) = - rho(ii) * d(ii) * ped c..ete and eve are for E, elec here is rho*E as is e ee(ii) = erho / rho(ii) eet(ii) = erhot / rho(ii) eev(ii) = erho - d(ii) * erhod c pe(ii) = pee + pdeg(ii) c ee(ii) = ee(ii) + edeg(ii) c pev(ii) = pev(ii) - pdeg(ii)*rho(ii)*gamdeg c eev(ii) = eev(ii) - pdeg(ii) pdv(ii) = - pdeg(ii)*rho(ii)*gamdeg edv(ii) = - pdeg(ii) eta(ii) = etaele etabd(ii) = detadd etabt(ii) = detadt 50 continue 1000 format(2i3,1p7e15.7) return end subroutine nrfd(dz,tz,pz,ptz,vpvz,ez,etz,vevz,s,si,k) implicit none save include 'arnett_eos.dek' c.. c..fermi-dirac integrals n = 1/2,3/2 for nonrelativistic c..electron eos c..e+e- pair gas added assuming zero chemical potential c..(degeneracy effects must be included through eos table) c.. c..declare integer k double precision dz,tz,pz,ptz,vpvz,ez,etz,vevz,s,si, 1 b,b3h,xj2,xj3,phi,b5h,a1, 2 fact,ppair,ptpair,pvpair,epair,etpair,evpair,spair, 3 cfact,bfact,afact,yp,ypt,gam,pradi,chem, 4 rgas,pairc,root2 parameter( rgas=8.3170000E+07, pairc = 4.93d17 ) parameter( root2 = 1.414213562d0 ) c---------------------------------------------------------------- c..si is ue/kT c..s is S/avagadro*k b = 5.929d9 / tz b3h = b*sqrt( b ) c..estimate fermi integral j2 for no positrons xj2 = b3h*dz / 2.9208d6 /root2 if( xj2 .lt. bj2( 1) )then call loi(xj2,xj3,si,phi) elseif( xj2 .gt. bj2(nerfd) )then call hii(xj2,xj3,si,phi) else call midi(xj2,xj3,si,phi) endif b5h = b*b3h a1 = 1.43391d24 fact = xj3*a1/b5h *2.0d0*root2 if( b .gt. 100.0d0 )then ppair = 0.0d0 ptpair = 0.0d0 pvpair = 0.0d0 epair = 0.0d0 etpair = 0.0d0 evpair = 0.0d0 spair = 0.0d0 else c..approximate relativistic pair gas component c..1.803 is J2(0) c..0.678 is J1.5(0) c..2.506 is sqrt(2*pi) cfact = 0.678094d0*sqrt(2.0d0*3.14159d0) bfact = dz*b**1.5d0/( cfact * 2.922d6 * exp(-b) ) c bfact = d*b**3/( 1.803d0 * 2.922d6 ) afact = 0.5d0*( bfact + sqrt( bfact**2 + 4.0d0) ) c..yp is Y*rho for pairs * 2 = rho*( Y+ + Y- ) c yp = 2.0d0 * 1.803d0 * 2.922d6/ afact /b**3 yp = 2.0d0 * cfact * 2.922d6/ afact /b**3 *exp(-b) c..NEED TO FIX DERIVATIVES HERE ypt = yp* ( 2.0d0*b + 3.0d0 )/tz c..J3(0)/[3 J2(0)] = 1.05051 yp = 1.05051d0*yp ppair = rgas * yp * tz ptpair = (ypt*tz + yp)*rgas pvpair = 0.0d0 gam = 1.5d0 epair = yp*( gam*rgas*tz + pairc ) etpair = ypt*(gam*rgas*tz + pairc ) + gam*rgas*yp evpair = 0.0d0 spair = (epair + ppair)/(rgas*dz*tz) endif pz = fact/3.0d0 + ppair ptz = 0.5d0*pz*(5.0d0 - 3.0d0*phi)/tz +ptpair vpvz = -pz * phi +pvpair c..energy per unit volume ez = 1.5d0 * pz + epair etz = 1.5d0 * ptz + etpair c..derivative per unit mass vevz = 1.5d0 * vpvz + ez + evpair c..entropy per unit Ye s = 5.0d0*xj3/(3.0d0*xj2) - si + spair c if( t .ge. 1.0d10 .and. d .le. 1.0d-4 )then c if( t .eq. 1.0d9 .and. d .eq. 1.0d-4 )then c write(*,'(8a12)')'fact/3','ppair','yp','beta','e(-2b)' c 1 ,'si','e(si)' c write(*,'(1p8e12.3)')fact/3.0d0, ppair,yp,b, c 1 exp(-2.0d0*b),si,exp(si) c write(*,'(8a12)')'fact','epair','e/p','T','d' c write(*,'(1p8e12.3)')fact, epair,epair/ppair,t,d c write(*,'(8a12)')'e(mu/kT)','bfact','Prad','pair/rad' c 1 ,'chem/kT' c pradi = 0.75641d-14/3.0*t**4 c chem = log( afact ) c write(*,'(1p8e12.3)')afact,bfact,pradi,ppair/pradi,chem c stop c endif return end subroutine hii(xj2,xj3,si,phi) implicit none save c.. c..degeneracy limit c..nonrelativistic fermi dirac functions j1/2, j3/2 c.. c..declare integer n double precision xj2,xj3,si,phi,pio8,u,a,bx,cx,si2,tsi2, 1 pi2,twothr parameter (pi2 = 9.8696044d0, twothr = 0.6666666666667d0 ) c..input: xj2 c..output: si = (mu -mc**2)/kT, xj3, phi=(3)j2**2/j1*j3 pio8 = pi2/8.0d0 u = 0.0d0 a = 1.5d0 * xj2 u = a**twothr do n = 1,5 u = ( a /(1.0d0 + pio8/(u*u) ) )**twothr enddo si = u si2 = si*si tsi2 = 3.0d0*si2 bx = 1.0d0 + pio8/ si2 cx = 1.0d0 + pio8/tsi2 xj3 = 0.6d0*xj2*si*bx/cx c..exact.coef. phi = ( 1.0d0 - pi2/tsi2 )/0.6d0 return end subroutine midi(xj2,xj3,si,phi) implicit none save include 'arnett_eos.dek' c..intermediate values of fermi dirac functions j1/2, j3/2 c..from table c..input xj2, asi, bj2, bj3, bphi c..output si, xj3, phi c.. c..declare integer j,k,l,i double precision xj2,xj3,si,phi,a,b,c,dj2p,dj2m,dj2,del,dj3p, 1 dj3m,dasip,dasim j = 1 do 100 k = 2, 29 j = k if( xj2 .lt. bj2(k) )go to 101 100 continue 101 continue c..linear interpol. for (3)j2**2/(j1*j3) l = j - 1 a = bphi(j) - bphi(l) b = bj2(j) - bj2(l) c = xj2 - bj2(l) phi = a*c/b + bphi(l) if( j .gt. 28 )j = 28 l = j - 1 i = j + 1 c..quadratic interpolation (17.2.81) for j3 and si dj2p = bj2(i) - bj2(j) dj2m = bj2(j) - bj2(l) dj2 = dj2p + dj2m del = xj2 - bj2(j) c..j3 interpolation dj3p = (bj3(i) - bj3(j))/dj2p dj3m = (bj3(j) - bj3(l))/dj2m a = bj3(j) b = (dj3m*dj2p + dj3p*dj2m)/dj2 c = (dj3p - dj3m)/dj2 xj3 = a + (b + c*del)*del dasip = (bsi(i) - bsi(j))/dj2p dasim = (bsi(j) - bsi(l))/dj2m a = bsi(j) b = (dasim*dj2p + dasip*dj2m)/dj2 c = (dasip - dasim)/dj2 si = a + (b + c*del)*del return end subroutine loi(xj2,xj3,si,phi) implicit none save c..nondegenerate, nonrelativistic fermi dirac j1/2, j3/2 c.. c..declare double precision xj2,xj3,si,phi,esi,a,b,root2,pi,rootpi parameter( root2 = 1.414213562d0, pi = 3.14159265359d0 ) parameter( rootpi = 1.772453851d0 ) c--------------------------------------------------------- c..input: xj2 (J1/2) c..output: si, xj3 (J3/2), phi (3[J1/2]**2/Jm1/2*J3/2) esi = 2.0d0 * xj2/( rootpi - xj2/root2 ) si = log( esi ) a = 1.0d0 - esi/root2/4.0d0 b = 1.0d0 - esi/root2/2.0d0 xj3 = 1.5d0 * xj2 * a / b phi = 1.003205d0 return end