program teos implicit none save include 'vector_eos.dek' c.. c..tests the weaver 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 eskep c..write out the results call pretty_eos_out('weaver:') 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 eskep implicit none save include 'vector_eos.dek' c.. c..given the temperature tt (deg k), the density dd (gr/cc), the average ion c..charge zbar, the average ion weight abar (amu) and a coulomb correction c..multiplier coulmult this routine computes the pressure (erg/cc), c..energy (erg/gr), entropy (dimensionless) and their derivatives with c..respect to temperature and density. c.. c..this routine assumes planck photons, non-relativistic but perhaps c..degenerate gas ions, possible relativistic and/or degenerate gas c..electrons and electron pairs. c.. c..the non-relativistic ion degeneracy is taken into account in terms of an c..assumed effective ion mass of xmimult times the neutron mass and a c..statistical weight of 2. c.. c..the effects of derivatives of the electron abundance (gotten from sdot) c..are also taken into account in calculating the electron energy and pressure c..derivatives. c.. c..other output quantities... c.. pt = dp/dt (ergs/cc/deg) c.. pd = dp/dd (ergs/g) c.. et = de/dt (ergs/g/deg) c.. ed = de/dd (ergs-cc/g/g) c.. c..declare external funcne integer nsev parameter (nsev=7) integer j double precision clight,pi,h,me,kerg,hion,ev2erg,avo,ssol,asol, 1 con2 parameter (clight = 2.99792458d10, 1 pi = 3.1415926535897932384d0, 2 h = 6.6260755d-27, 3 me = 9.1093897d-28, 4 kerg = 1.380658d-16, 5 hion = 13.605698140d0, 6 ev2erg = 1.602d-12, 7 avo = 6.0221367d23, 8 ssol = 5.67051d-5, 9 asol = 4.0d0 * ssol / clight, & con2 = clight * clight) c.. c..note b=1.0/sqrt(8.0*pi) double precision cpf,cbeta,fa0,cxne,c2mec2,b parameter (fa0=64.0/(9.0*pi), + cpf=3.*h**3/(8.*pi*me**3*clight**3), + cbeta=kerg/(me*clight**2), + cxne=2.488668349711d+30, + c2mec2=2.*me*clight**2, + b=0.199471140) c.. integer npflag,icalcne,jj double precision tt,dd,zbar,abar,ccmult double precision eee,pee,xnebe,ept, + pe,ped,pet,pd,pt,xnee,xnebd,xnebt, + xnpp,xnpbep,xnpbt,xnpbe, + epep,ppt,ppep,pp,ep, + ionmult,radmult,elemult,dzero,xmimult, + ab,zb,ye,yebt,yebd,t,d,abar1,abar1bd,abar1bt, + xnecon,xneconv, + sig1,p,pe1,e1,sigion,fa,ga,gae,sigp, + xnetop,xnemc,xnemcbt,xnemcbe,etatop c..fot the derivatives double precision chit,chid,cv,cp,gam1,gam2,gam3,nabad,sound,zz,z c..eff charge double precision xip,x,f,s,r,ze,zlt,zld c..rad component double precision pr,prt,prd,er,ert,erd,sigr c..ion component double precision xni,xnibt,xnibd,xmi,tmkt,f12,sigi,f32,g1,pi1, + pi1t,pi1d,ef,rx,xx,ei,eit,eid,eionmult c..coulomb correction double precision bcoul,abar123,pcoul,ecoul,pcouldd,pcouldt, + ecouldd,ecouldt c..electron component double precision beta,betai,xnem,xnembd,xnembt,xeta,xetal,cdegl, + funcgx,pf2,fxx,ee,eed,eet,et,ed c..eta root find integer npfla466,jroot double precision xncom,tcom,rootnr,sige,etap,etapbeta,eta, + etabd,etabt,etapbt,etaconv,xnpair c..thermodynamic consistency check double precision dse,dpe,dsp c.. c..communication to kappa and funcne common /kpar/ pe,ped,pet,xnee,xnebd,xnebt,eta,etabd,etabt c..communication to funcne the eta root find routine common /compari/ npfla466 common /comparr/ xncom,tcom c.. c..communication to a top level routine c common /answer/ e1,ei,ee,er,p,pi1,pe1,pr,sig1,sigi,sige,sigr, c 1 ze,xnetop,etatop,ga,fa,gae,radmult, c 2 ionmult,elemult,pd,pt,ed,et,pp,ep,xnpp,etap, c 3 sigp,pcoul,ecoul,sound,ccmult c.. c..these are the fermi integrals, denoted f12a and 2/3 * f32a are are given c..by clayton pg 94-100; the integer i stores the last table entry used; c..powi contains the interpolating power integer iat c common /esdat0/ iat double precision f12a(200),f32a(200),powi(199) c common /esdat1/ f12a,f32a,powi data iat /0/ data (f12a(jj), jj=1,50) / 1 0.016128, 0.017812, 0.019670, 0.021721, 0.023984, 2 0.026480, 0.029233, 0.032269, 0.035615, 0.039303, 3 0.043366, 0.047842, 0.052770, 0.058194, 0.064161, 4 0.070724, 0.077938, 0.085864, 0.094566, 0.104116, 5 0.114588, 0.126063, 0.138627, 0.152373, 0.167397, 6 0.183802, 0.201696, 0.221193, 0.242410, 0.265471, 7 0.290501, 0.317630, 0.346989, 0.378714, 0.412937, 8 0.449793, 0.489414, 0.531931, 0.577470, 0.626152, 9 0.678094, 0.733403, 0.792181, 0.854521, 0.920505, + 0.990209, 1.063694, 1.141015, 1.222215, 1.307327/ data (f12a(jj), jj=51,100) / 1 1.396375, 1.489372, 1.586323, 1.687226, 1.792068, 2 1.900833, 2.013496, 2.130027, 2.250391, 2.374548, 3 2.502458, 2.634072, 2.769344, 2.908224, 3.050659, 4 3.196598, 3.345988, 3.498775, 3.654905, 3.814326, 5 3.976985, 4.142831, 4.311811, 4.483876, 4.658977, 6 4.837066, 5.018095, 5.202020, 5.388795, 5.578378, 7 5.77072, 5.96580, 6.16356, 6.36396, 6.56698, 8 6.77257, 6.98070, 7.19134, 7.40445, 7.62001, 9 7.83797, 8.05832, 8.28103, 8.50606, 8.73339, + 8.96299, 9.19485, 9.42893, 9.66521, 9.90367/ data (f12a(jj), jj=101,150) / 1 10.14428, 10.38703, 10.63190, 10.87886, 11.12789, 2 11.37898, 11.63211, 11.88726, 12.14440, 12.40354, 3 12.66464, 12.92769, 13.19267, 13.45958, 13.72839, 4 13.99910, 14.27168, 14.54612, 14.82241, 15.10053, 5 15.38048, 15.66224, 15.94580, 16.23114, 16.51826, 6 16.80714, 17.09776, 17.39013, 17.68423, 17.98004, 7 18.27756, 18.57677, 18.87768, 19.18026, 19.48451, 8 19.79041, 20.09796, 20.40715, 20.71797, 21.03042, 9 21.34447, 21.66013, 21.97738, 22.29622, 22.61664, + 22.93862, 23.26217, 23.58728, 23.91393, 24.24212/ data (f12a(jj), jj=151,200) / 1 24.57184, 24.90309, 25.23586, 25.57013, 25.90591, 2 26.24319, 26.58195, 26.92220, 27.26393, 27.60712, 3 27.95178, 28.29789, 28.64545, 28.99446, 29.34491, 4 29.69679, 30.05009, 30.40482, 30.76096, 31.11851, 5 31.47746, 31.83781, 32.19956, 32.56268, 32.92720, 6 33.29308, 33.66034, 34.02896, 34.39894, 34.77028, 7 35.14297, 35.51700, 35.89238, 36.26908, 36.64712, 8 37.02649, 37.40718, 37.78918, 38.17250, 38.55712, 9 38.94304, 39.33027, 39.71879, 40.10859, 40.49969, + 40.89206, 41.28571, 41.68064, 42.07683, 42.47429/ data (f32a(jj), jj=1,50) / 1 0.016179, 0.017875, 0.019748, 0.021816, 0.024099, 2 0.026620, 0.029404, 0.032476, 0.035868, 0.039611, 3 0.043741, 0.048298, 0.053324, 0.058868, 0.064981, 4 0.071720, 0.079148, 0.087332, 0.096347, 0.106273, 5 0.117200, 0.129224, 0.142449, 0.156989, 0.172967, 6 0.190515, 0.209777, 0.230907, 0.254073, 0.279451, 7 0.307232, 0.337621, 0.370833, 0.407098, 0.446659, 8 0.489773, 0.536710, 0.587752, 0.643197, 0.703351, 9 0.768536, 0.839082, 0.915332, 0.997637, 1.086358, + 1.181862, 1.284526, 1.394729, 1.512858, 1.639302/ data (f32a(jj), jj=51,100) / 1 1.774455, 1.918709, 2.072461, 2.236106, 2.410037, 2 2.594650, 2.790334, 2.997478, 3.216467, 3.447683, 3 3.691502, 3.948298, 4.218438, 4.502287, 4.800202, 4 5.112536, 5.439637, 5.781847, 6.139503, 6.512937, 5 6.902476, 7.308441, 7.731147, 8.170906, 8.628023, 6 9.102801, 9.595535, 10.106516, 10.636034, 11.184369, 7 11.75180, 12.33860, 12.94505, 13.57140, 14.21793, 8 14.88489, 15.57253, 16.28111, 17.01088, 17.76208, 9 18.53496, 19.32976, 20.14671, 20.98604, 21.84799, + 22.73279, 23.64067, 24.57184, 25.52653, 26.50495/ data (f32a(jj), jj=101,150) / 1 27.50733, 28.53388, 29.58481, 30.66033, 31.76065, 2 32.88598, 34.03652, 35.21247, 36.41404, 37.64142, 3 38.89481, 40.17441, 41.48041, 42.81301, 44.17239, 4 45.55875, 46.97227, 48.41315, 49.88156, 51.37769, 5 52.90173, 54.45385, 56.03424, 57.64307, 59.28052, 6 60.94678, 62.64201, 64.36639, 66.12009, 67.90329, 7 69.71616, 71.55886, 73.43157, 75.33445, 77.26768, 8 79.23141, 81.22582, 83.25106, 85.30730, 87.39471, 9 89.51344, 91.66365, 93.84552, 96.05918, 98.30481, + 100.58256, 102.89259, 105.23505, 107.61010, 110.01789/ data (f32a(jj), jj=151,200) / 1 112.45857, 114.93231, 117.43924, 119.97953, 122.55332, 2 125.16076, 127.80201, 130.47720, 133.18650, 135.93004, 3 138.70797, 141.52044, 144.36760, 147.24958, 150.16654, 4 153.11861, 156.10594, 159.12868, 162.18696, 165.28092, 5 168.41071, 171.57646, 174.77831, 178.01642, 181.29090, 6 184.60190, 187.94956, 191.33401, 194.75540, 198.21385, 7 201.70950, 205.24249, 208.81295, 212.42101, 216.06681, 8 219.75048, 223.47215, 227.23196, 231.03003, 234.86650, 9 238.74150, 242.65515, 246.60759, 250.59895, 254.62936, + 258.69893, 262.80781, 266.95612, 271.14398, 275.37153/ c.. c..the powi the first time through if (iat .eq. 0) then do 5 iat=1,199 powi(iat)=log(f32a(iat+1)/f32a(iat))/log(f12a(iat+1)/f12a(iat)) 5 continue iat=1 end if radmult = 1.0d0 ionmult = 1.0d0 elemult = 1.0d0 ccmult = 0.0d0 eionmult = 0.0d0 c..pipeline loop eosfail = .false. do j=jlo_eos,jhi_eos tt = temp_row(j) dd = den_row(j) abar = abar_row(j) zbar = zbar_row(j) c.. c..initialize eet=0.0 eee=eet pet=eee pee=pet xnebt=pee xnebe=xnebt ee=xnebe pe=ee pe1 = pe xnee=pe ept=0.0 epep=ept ppt=epep ppep=ppt xnpbep=ppep xnpbt=xnpbep ep=xnpbt pp=ep xnpp=pp xnecon=xneconv c..turns on/off the density (pressure) ionization term in the saha equation c.. dzero=0.1 dzero=1.0d-20 c.. npflag=1 icalcne=0 xmimult=1.0 ab=abar zb=zbar ye=zb/ab yebt=0. yebd=0. t=tt d=dd abar1=ab abar1bd=0. abar1bt=0. pcoul=0. ecoul=0. pcouldd=0. pcouldt=0. ecouldd=0. ecouldt=0. etap = 0.0d0 xnpp = 0.0d0 c.. c..calculate the effective charge ze and its temp and density derivatives c..the parameter constant c1 = 0.5 * ( (h*h)/(2.0*pi*me) )**(1.5) c..saha = z0/z1 * ne * c1 * beta**(1.5) * exp(eion*beta) c..set dzero = 1e-20 for complete ionization; default is 0.1 xip=hion * ev2erg*zb x=(h*h)/(2.*pi*me*kerg*t) f=4.*zb*avo/ab*d*.5*x*sqrt(x)*exp(xip/(kerg*t+.005*xip)-d/dzero) s=sqrt(1.+f) r=1./(1.+s) ze=2.*zb*r zlt=r*.5/s*f*(1.5+xip*kerg*t/(kerg*t+.005*xip)**2)/t zld=r*.5/s*f*(1./dzero-1./d) c.. c..radiation component pr=1./3.*asol*t**4 prt=4.*pr/t prd=0. er=asol/d*t**4 ert=4.*er/t erd=-er/d sigr=1.213d-22*t*t*t/d c.. c..ion component, assumes a statistical weight average of 2 and that the c..particle mass involved in calculating the degeneracy is xmimult times c..the neutron mass xni=avo/abar1*d xnibt=-avo*d*abar1bt/(abar1*abar1) xnibd=avo/abar1*(1.-d*abar1bd/abar1) xmi=1.67482d-24*xmimult tmkt=2.*xmi*kerg*t f12=(h*h*h)/(4.*pi)*xni/(tmkt*sqrt(tmkt)) c.. c.. c..note 2.379 = 5/2 + ln(sqrt(pi1/4)) sigi=(log(abar*sqrt(abar)/f12)+2.379)/abar c..turn off ion degeneracy, comment this line if you want partial c..degeneracy of the ions taken into account f12 = 0.0d0 c.. c.... check if table hit is still good if(f12-f12a(iat)) 41,42,42 42 if(f12-f12a(iat+1)) 43,43,44 c.... use table data 43 f32=f32a(iat)*(f12/f12a(iat))**powi(iat) g1=f32/f12 pi1=xni*kerg*t*g1 pi1t=xni*kerg*g1*(2.5-1.5*powi(iat)+t*powi(iat)*xnibt/xni) pi1d=kerg*t*g1*powi(iat)*xnibd go to 60 41 if(iat-1) 45,46,45 45 iat=iat-1 if(f12-f12a(iat)) 41,43,43 44 if(iat-199) 47,48,47 47 iat=iat+1 if(f12-f12a(iat+1)) 43,43,44 c.. c..near perfect gas 46 x=b*f12 c..turn off ion degeneracy, comment this line if you want partial c..degeneracy of the ions taken into account x = 0.0d0 pi1=xni*kerg*t*(1.+x) pi1t=xni*kerg*(1.-.5*x+(1.+2.*x)*t*xnibt/xni) pi1d=kerg*t*(1.+2.*x)*xnibd go to 60 c.. c.... near complete degeneracy 48 ef=(h*h)/(8.*xmi)*(3./pi*xni)**(2./3.) x=ef/(kerg*t) rx=(kerg*t)/ef pi1=xni*kerg*t*(.4*x+(pi*pi)/6.*rx) xx=2./3.*x+(pi*pi)/18.*rx pi1t=xni*kerg*((pi*pi)/3.*rx+xx*t*xnibt/xni) pi1d=kerg*t*xx*xnibd c.. c..this is valid only in the hard degeneracy region c..coulomb corrections to pressure and energy. c..use wigner-seitz approximation p. 152 clayton c..approximate z53bar by zbar53 c..note bcoul=-0.3*(4 pi1/3)**(1/3)*na**(4/3)*e**2 60 continue bcoul=-5.6739e+12 abar123=abar1**(2./3.) pcoul=bcoul*ye*ye*abar123*d**(4./3.) ecoul=3.*pcoul/d pcouldd=pcoul*(4./(3.*d)+(2./ye)*yebd + +(2./(3.*abar1))*abar1bd) pcouldt=pcoul*((2./ye)*yebt+(2./(3.*abar1))*abar1bt) ecouldd=3.*pcouldd/d-3.*pcoul/(d*d) ecouldt=3.*pcouldt/d c.. c.... figure ion energy ei=1.5*pi1/d eit=1.5*pi1t/d eid=1.5/d*(pi1d-pi1/d) c.. c.. c..electron component; calculate number of free matter electrons assuming c..that ionization is complete when beta gt 0.1 i.e. t gt about 10**8 beta=cbeta*t betai=1./beta xnem=d*avo*ze*ye/zb xnetop = xnem xnembd=xnem*(1./d+zld+yebd/ye) xnembt=xnem*(zlt+yebt/ye) c.. c..estimate the fermi degeneracy parameter, eta, from its c..asymptotic limits (c.f. cox and giuli, chapt. 24) if(beta.ge.1.)xnpair=2.345d+30*beta**2 if(beta.lt.1.)xnpair=2.546d+30*beta*(1.+.75*beta)*exp(-betai) if(xnpair.lt.xnem) go to 90 xeta=-0.5 go to 100 90 cdegl=log10(xnem**0.6/t) 95 if(cdegl.gt.9.5) go to 110 funcgx=(1.+fa0*beta)*sqrt(1.+fa0*beta*.5)-1. xeta=-log(cxne*beta*sqrt(beta*pi)*0.5*(1.+1.5*funcgx/fa0)/xnem) if(cdegl.lt.8.5) go to 100 xetal=xeta 110 pf2=(xnem*cpf)**(2./3.) if(pf2.lt.1.d-6)xeta=pf2*(1.-pf2*0.25)*.5*betai if(pf2.ge.1.d-6)xeta=betai*(sqrt(pf2+1.)-1.) if(cdegl.gt.9.5) go to 100 120 fxx=9.5-cdegl xeta=xetal*fxx+xeta*(1.-fxx) c.. c..solve for the exact value of eta using the root-finding subroutine rootnr 100 xncom=xnem tcom=t npfla466=npflag jroot=49 etaconv = 1.0d-6 eta=rootnr(funcne,xeta,etaconv,jroot) etatop = eta c write(6,111) eta c 111 format(1x,1p6e14.6) c..calculate the negaton contribution call emeos(t,eta,nsev,xnee,pe,ee,xnebe,xnebt,pee,pet,eee,eet) c write(6,111) eta,xnee,xnem,xnee-xnem,xnee/xnem pe1 = pe sige=((pe+ee)/(kerg*t)-eta*xnee)/(d*avo) c..calculate the positron contribution and add to negaton values c..if temperature is not too low if((beta.le.0.02).or.(npflag.le.0)) go to 130 etap=-eta-2./beta etapbeta=-1. etapbt=2./(beta*t) if(etap.le.-200.)etapbeta=0. if(etap.le.-200.)etapbt=0. if(etap.le.-200.)etap=-200. call emeos(t,etap,nsev,xnpp,pp,ep,xnpbep,xnpbt,ppep,ppt,epep,ept) xnee=xnee+xnpp xnemc=xnee-xnpp pe=pe+pp pe1 = pe ee=ee+ep+xnpp*c2mec2 xnpbe=xnpbep*etapbeta xnpbt=xnpbt+xnpbep*etapbt xnemcbe=xnebe-xnpbe xnebe=xnebe+xnpbe xnemcbt=xnebt-xnpbt xnebt=xnebt+xnpbt pee=pee+ppep*etapbeta pet=pet+ppt+ppep*etapbt eee=eee+epep*etapbeta+c2mec2*xnpbe eet=eet+ept+epep*etapbt+c2mec2*xnpbt sigp=((pp+ep)/(kerg*t)-etap*xnpp)/(d*avo) go to 140 c.. c..no pair case 130 xnemc=xnee xnemcbe=xnebe xnemcbt=xnebt sigp=0. c.. c..calculate the final density and temperature derivatives 140 etabd=xnembd/xnemcbe etabt=(-xnemcbt+xnembt)/xnemcbe pet=pet+pee*etabt eet=eet+eee*etabt ped=pee*etabd eed=eee*etabd xnebt=xnebt+xnebe*etabt xnebd=xnebe*etabd c..convert energy density to specific energy c..and correct for ionization potential energy x=xip*avo/ab ee=ee/d+x*ze*eionmult eed=(eed-ee)/d+x*ze*(zld+1./d)*eionmult eet=eet/d+x*ze*zlt*eionmult sigion=xip*ze*eionmult/(kerg*t*abar) c..total all components p=pi1*ionmult+pr*radmult+pe*elemult + pcoul*ccmult e1=ei*ionmult+er*radmult+ee*elemult + ecoul*ccmult pt=pi1t*ionmult+prt*radmult+pet*elemult + pcouldt*ccmult pd=pi1d*ionmult+prd*radmult+ped*elemult + pcouldd*ccmult et=eit*ionmult+ert*radmult+eet*elemult + ecouldt*ccmult ed=eid*ionmult+erd*radmult+eed*elemult + ecouldd*ccmult sig1=sigi+sige+sigp+sigion+sigr c.. c..calculate adiabatic indices ga=d/p*(pt/et*(p/(d*d)-ed)+pd) fa=(p/(d*d)-ed)*d/(et*t*ga) gae=d/pe*(pet/eet*(pe/(d*d)-eed)+ped) c.. 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/dd chit = tt/p * pt chid = pd/zz cv = et x = zz * chit/(tt * 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 + (e1 + con2)/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.0e20 c dse = tt*st/et - 1.0d0 dpe = (ed*d**2 + tt*pt)/p - 1.0d0 dsp = 1.0e20 c dsp = -sd*(d**2/pt) - 1.0d0 c..end of vector loop ptot_row(j) = p dpt_row(j) = pt dpd_row(j) = pd dpa_row(j) = 0.0d0 dpz_row(j) = 0.0d0 etot_row(j) = e1 det_row(j) = et ded_row(j) = ed dea_row(j) = 0.0d0 dez_row(j) = 0.0d0 stot_row(j) = sig1 * 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) = pr erad_row(j) = er srad_row(j) = sigr * kerg * avo pion_row(j) = pi1 eion_row(j) = ei sion_row(j) = sigi * kerg * avo xni_row(j) = xni pele_row(j) = pe-pp ppos_row(j) = pp dpept_row(j) = pet dpepd_row(j) = ped dpepa_row(j) = 0.0d0 dpepz_row(j) = 0.0d0 eele_row(j) = (ee*dd - ep - xnpp*c2mec2)/dd epos_row(j) = (ep+xnpp*c2mec2)/dd deept_row(j) = eet deepd_row(j) = eed deepa_row(j) = 0.0d0 deepz_row(j) = 0.0d0 sele_row(j) = sige * kerg * avo spos_row(j) = sigp * kerg * avo dsept_row(j) = 0.0d0 dsepd_row(j) = 0.0d0 dsepa_row(j) = 0.0d0 dsepz_row(j) = 0.0d0 xnem_row(j) = xnem xne_row(j) = xnee - xnpp dxnet_row(j) = xnebt dxned_row(j) = xnebd dxnea_row(j) = 0.0d0 dxnez_row(j) = 0.0d0 xnp_row(j) = xnpp etaele_row(j) = eta detat_row(j) = etabt detad_row(j) = etabd detaa_row(j) = 0.0d0 detaz_row(j) = 0.0d0 etapos_row(j) = etap pcou_row(j) = pcoul ecou_row(j) = ecoul 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 emeos(t,eta,ixne,xne,pe,ee,xnebe,xnebt,pee,pet,eee,eet) implicit none save c..this routine computes the electron (or positron) number density, c..electron (or positron) pressure and the electron (or positron) c..kinetic energy density given a value for the degeneracy parameter (=mu/kt) c..and the temperature. the method is based on divine, ap. j. 142,1652(1965) c..and is valid for arbitrary degenerate and relativistic situations. c.. c..required input c..t = temperature (deg k) c..eta = parameter (mu/kt) for electron type (follows chiu's convention) c..ixne = integer flag; ixne=4 compute xne only; ixne=7, compute all c.. c..output quantities c..xne = electron number density (1/cc) c..pe = electron pressure (erg/cc) c..ee = electron energy density (erg/cc) c..xnebe = dxne/deta (1/cc) c..xnebt = dxne/dt (1/cc/degk) c..pee = dpe/deta (erg/cc) c..pet = dpe/dt (erg/cc/degk) c..eee = dee/deta (erg/cc) c..eet = dee/dt (erg/cc/degk) c.. c.. c..notes: eta is additive inverse of divine's and clayton's alpha. c.. fe(1)=f(1/2),fe(2)=f(1),fe(3)=f(3/2), etc. are fermi integrals c.. except fe(7)=f(-1/2). c.. c..declare integer ixne double precision t,eta,xne,pe,ee,xnebe,xnebt,pee,pet,eee,eet double precision clight,pi,h,me,kerg parameter (clight = 2.99792458d10, 1 pi = 3.1415926535897932384d0, 2 h = 6.6260755d-27, 3 me = 9.1093897d-28, 4 kerg = 1.380658d-16) integer nfpt,nfpt0,itab,ii,ij,i1,i2,i3,l,lf parameter (nfpt=67*7, nfpt0=41*7) double precision eqfm(nfpt),eqfmbe(nfpt),eqfx(nfpt0), + eqfxbe(nfpt0),gam(7),xfe(7),febe(7),fe(7), + afac(7),pie2,c2,c4,c6,c10,c12,c14,c16,c22, + c32,c34,c36,c42,c50,c52,c54,c56,c62,c64,c72, + c74,c76,cbeta,cxne,cee,cpe, + zero,eeta,xk,fbar,etar,eta2,eta2i,fe1,fe3,beta, + betar,beta32,u1,u1i,funu1r,ffunu1,betai, + beta52,u3,u3i,u3i2,funu3r,funu3l,ffunu3, + gfunu3,betabt,u1beta,f1bbeta,x,funu1l,f1beta, + u3beta,f3beta,g3beta,f3bbeta,g3bbeta parameter (pie2=pi*pi) parameter (c2=pie2/48.0, c4=c2*pie2*7.0/80.0, + c6=c4*pie2*1395.0/3528.0, + c10=2.0/3.0, c12=6.0*c2, + c14=6.0*c4, c16=14.0*c6, + c22=8.0*c2, c32=30.0*c2, + c34=-10.0*c4, c36=-10.0*c6, + c42=16.0*c2, c50=5.0/7.0, + c52=70.0*c2, c54=70.0*c4, + c56=14.0*c6, c62=24.0*c2, + c64=64.0*c4, c72=-2.0*c2, + c74=-10.0*c4, c76=-42.0*c6) parameter (cbeta=kerg/(me*clight**2), + cxne=2.488668349711d+30, + cee=cxne*me*clight**2, cpe=2.*cee/3.) c.. c..note: cxne=sqrt(2.)*pi*(me*c/(hbar*pi))**3 c.. c..common link to ecapnuc c..this common block initialization must be present for the routine c..to operate correctly. rigorous compilers will, of course, complain c..since its non-standard fortran integer it double precision etai(67),fm(67,7),fmbe(67,7),etax(41),fx(41,7), + fxbe(41,7) common /esdat/ fm,fmbe,fx,fxbe,etai,etax,it equivalence (fm,eqfm), (fmbe,eqfmbe), (fx,eqfx), + (fxbe,eqfxbe) c..here are the tables of reduced fermi intergrals. c..fm(i,j) is the fermi intergral of order j/2, evaluated c..at eta=etai(i), and divided by fbar=((eta+.16)**(j/2+1.))/(j/2+1.)+xfe(j) c..to remove most of the eta dependence. c.. data xfe /0.5966,0.7178,1.0000,1.5570,2.6527,4.8771,0./ data etai / -.16, -.14, -.12, -.10, -.08, -.06, -.04, -.02, 0., 1 .02, .04, .06, .08, .10, .15, .20, .25, .30, 2 .35, .40, .45, .50, .60, .70, .80, .90, 1.0, 3 1.2, 1.4, 1.6, 1.8, 2., 2.2, 2.4, 2.6, 2.8, 4 3., 3.2, 3.4, 3.6, 3.8, 4., 4.2, 4.4, 4.6, 5 4.8, 5., 5.5, 6., 6.5, 7., 7.5, 8., 8.5, 6 9., 9.5, 10., 11., 12., 13., 14., 15., 16., 7 17., 18., 19., 20./ data ((fm(ii,ij),ij=1,6),ii=1,8)/ 1 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 1 1.013117, 1.017019, 1.018019, 1.018598, 1.019019, 1.019332, 1 1.023702, 1.033711, 1.036246, 1.037507, 1.038379, 1.039023, 1 1.032646, 1.050053, 1.054641, 1.056720, 1.058084, 1.059080, 1 1.040272, 1.066021, 1.073171, 1.076230, 1.078138, 1.079507, 1 1.046775, 1.081595, 1.091809, 1.096028, 1.098541, 1.100310, 1 1.052296, 1.096753, 1.110525, 1.116105, 1.119296, 1.121496, 1 1.056942, 1.111478, 1.129293, 1.136452, 1.140403, 1.143068/ data ((fm(ii,ij),ij=1,6),ii=9,27)/ 1 1.060804, 1.125753, 1.148086, 1.157058, 1.161862, 1.165032, 1 1.063959, 1.139562, 1.166876, 1.177910, 1.183673, 1.187393, 1 1.066474, 1.152892, 1.185638, 1.198998, 1.205836, 1.210155, 1 1.068409, 1.165732, 1.204346, 1.220306, 1.228347, 1.233322, 1 1.069819, 1.178072, 1.222971, 1.241821, 1.251204, 1.256898, 1 1.070751, 1.189904, 1.241489, 1.263528, 1.274404, 1.280886, 1 1.071282, 1.217227, 1.287143, 1.318522, 1.333873, 1.342678, 1 1.069707, 1.241296, 1.331563, 1.374329, 1.395365, 1.407097, 1 1.066516, 1.262129, 1.374367, 1.430636, 1.458755, 1.474156, 1 1.062109, 1.279799, 1.415197, 1.487097, 1.523879, 1.543844, 1 1.056811, 1.294431, 1.453728, 1.543342, 1.590531, 1.616121, 1 1.050882, 1.306187, 1.489680, 1.598976, 1.658465, 1.690915, 1 1.044534, 1.315257, 1.522814, 1.653592, 1.727391, 1.768118, 1 1.037936, 1.321854, 1.552948, 1.706778, 1.796978, 1.847581, 1 1.024496, 1.328517, 1.603747, 1.807238, 1.936605, 2.012468, 1 1.011317, 1.327970, 1.641677, 1.897303, 2.073955, 2.183452, 1 0.998845, 1.321904, 1.667116, 1.974436, 2.205241, 2.357597, 1 0.987320, 1.311809, 1.681060, 2.036854, 2.326565, 2.531177, 1 0.976851, 1.298932, 1.684912, 2.083663, 2.434260, 2.699802/ data ((fm(ii,ij),ij=1,6),ii=28,46)/ 1 0.959138, 1.268614, 1.668821, 2.131264, 2.597315, 3.002927, 1 0.945423, 1.236474, 1.631461, 2.125634, 2.681284, 3.230589, 1 0.935134, 1.205512, 1.583147, 2.081449, 2.689582, 3.358925, 1 0.927654, 1.177145, 1.531015, 2.013821, 2.637851, 3.384739, 1 0.922422, 1.151882, 1.479436, 1.934962, 2.546243, 3.323774, 1 0.918961, 1.129751, 1.430792, 1.853310, 2.433192, 3.201724, 1 0.916878, 1.110550, 1.386189, 1.774017, 2.312554, 3.044725, 1 0.915856, 1.093979, 1.345970, 1.699853, 2.193350, 2.873722, 1 0.915647, 1.079714, 1.310064, 1.632041, 2.080732, 2.702955, 1 0.916049, 1.067447, 1.278184, 1.570878, 1.977186, 2.540828, 1 0.916909, 1.056893, 1.249958, 1.516152, 1.883535, 2.391537, 1 0.918105, 1.047805, 1.224990, 1.467396, 1.799653, 2.256594, 1 0.919542, 1.039966, 1.202896, 1.424043, 1.724925, 2.135948, 1 0.921147, 1.033194, 1.183323, 1.385507, 1.658518, 2.028736, 1 0.922863, 1.027331, 1.165956, 1.351227, 1.599544, 1.933733, 1 0.924646, 1.022245, 1.150515, 1.320691, 1.547138, 1.849610, 1 0.926464, 1.017824, 1.136755, 1.293438, 1.500503, 1.775072, 1 0.928291, 1.013972, 1.124466, 1.269064, 1.458921, 1.708925, 1 0.930108, 1.010610, 1.113464, 1.247213, 1.421762, 1.650096/ data ((fm(ii,ij),ij=1,6),ii=47,65)/ 1 0.931901, 1.007670, 1.103591, 1.227577, 1.388472, 1.597648, 1 0.936219, 1.001808, 1.082993, 1.186513, 1.319187, 1.489339, 1 0.940240, 0.997555, 1.066964, 1.154422, 1.265376, 1.406086, 1 0.943934, 0.994446, 1.054322, 1.128990, 1.222936, 1.340998, 1 0.947305, 0.992161, 1.044231, 1.108573, 1.188987, 1.289306, 1 0.950372, 0.990477, 1.036087, 1.091988, 1.161479, 1.247666, 1 0.953160, 0.989239, 1.029449, 1.078374, 1.138932, 1.213693, 1 0.955697, 0.988333, 1.023989, 1.067089, 1.120258, 1.185658, 1 0.958009, 0.987678, 1.019463, 1.057654, 1.104645, 1.162284, 1 0.960120, 0.987214, 1.015683, 1.049704, 1.091482, 1.142618, 1 0.962052, 0.986897, 1.012507, 1.042958, 1.080298, 1.125934, 1 0.965455, 0.986575, 1.007538, 1.032243, 1.062491, 1.099400, 1 0.968346, 0.986528, 1.003918, 1.024252, 1.049147, 1.079518, 1 0.970828, 0.986649, 1.001236, 1.018177, 1.038939, 1.064293, 1 0.972976, 0.986871, 0.999224, 1.013483, 1.030995, 1.052418, 1 0.974850, 0.987152, 0.997699, 1.009808, 1.024719, 1.043007, 1 0.976499, 0.987466, 0.996535, 1.006895, 1.019698, 1.035448, 1 0.977958, 0.987796, 0.995641, 1.004564, 1.015635, 1.029305, 1 0.979258, 0.988131, 0.994954, 1.002682, 1.012316, 1.024260/ data ((fm(ii,ij),ij=1,6),ii=66,67)/ 1 0.980422, 0.988464, 0.994426, 1.001152, 1.009582, 1.020078, 1 0.981471, 0.988792, 0.994021, 0.999900, 1.007311, 1.016585/ data ((fm(ii,ij),ij=7,7),ii=1,67)/ + 1.0000000, 1.0199580, 1.0403098, 1.0610631, 1.0822256, 1.1038051, + 1.1258098, 1.1482479, 1.1711275, 1.1944574, 1.2182460, 1.2425021, + 1.2672348, 1.2924531, 1.3576847, 1.4261588, 1.4980308, 1.5734622, + 1.6526202, 1.7356781, 1.8228141, 1.9142114, 2.1105391, 2.3261794, + 2.5626401, 2.8213428, 3.1035136, 3.7412440, 4.4748891, 5.2861403, + 6.1300174, 6.9303090, 7.5886824, 8.0121513, 8.1485232, 8.0065453, + 7.6463371, 7.1499510, 6.5934455, 6.0326793, 5.5017276, 5.0175590, + 4.5859108, 4.2060742, 3.8741229, 3.5848527, 3.3328433, 2.8336372, + 2.4722872, 2.2047422, 2.0020915, 1.8453269, 1.7217470, 1.6226850, + 1.5420971, 1.4756789, 1.4203042, 1.3340125, 1.2706905, 1.2228846, + 1.1859375, 1.1568176, 1.1334818, 1.1145124, 1.0989003, 1.0859112, + 1.0750006/ c.. c..extended fermi integral table for eta=-4 to 0. data etax/-4.,-3.9,-3.8,-3.7,-3.6,-3.5,-3.4,-3.3,-3.2,-3.1,-3.0, 1 -2.9,-2.8,-2.7,-2.6,-2.5,-2.4,-2.3,-2.2,-2.1,-2.0, 2 -1.9,-1.8,-1.7,-1.6,-1.5,-1.4,-1.3,-1.2,-1.1,-1.0, 3 -.9,-.8,-.7,-.6,-.5,-.4,-.3,-.2,-.1,0./ data ((fx(ii,ij),ij=1,6),ii=1,15)/ 1 0.9935882, 0.9954580, 0.9967836, 0.9977229, 0.9983882, 0.9988594, 1 0.9929212, 0.9949845, 0.9964477, 0.9974848, 0.9982195, 0.9987399, 1 0.9921857, 0.9944622, 0.9960771, 0.9972220, 0.9980333, 0.9986080, 1 0.9913747, 0.9938861, 0.9956682, 0.9969320, 0.9978277, 0.9984623, 1 0.9904808, 0.9932508, 0.9952171, 0.9966119, 0.9976007, 0.9983014, 1 0.9894957, 0.9925503, 0.9947195, 0.9962587, 0.9973502, 0.9981238, 1 0.9884106, 0.9917782, 0.9941707, 0.9958690, 0.9970737, 0.9979278, 1 0.9872155, 0.9909273, 0.9935656, 0.9954392, 0.9967687, 0.9977114, 1 0.9858999, 0.9899899, 0.9928987, 0.9949652, 0.9964321, 0.9974726, 1 0.9844522, 0.9889576, 0.9921638, 0.9944426, 0.9960609, 0.9972091, 1 0.9828598, 0.9878212, 0.9913541, 0.9938665, 0.9956515, 0.9969184, 1 0.9811090, 0.9865705, 0.9904623, 0.9932317, 0.9952001, 0.9965978, 1 0.9791851, 0.9851948, 0.9894806, 0.9925323, 0.9947025, 0.9962442, 1 0.9770723, 0.9836823, 0.9884002, 0.9917620, 0.9941541, 0.9958543, 1 0.9747533, 0.9820201, 0.9872118, 0.9909140, 0.9935499, 0.9954245/ data ((fx(ii,ij),ij=1,6),ii=16,29)/ 1 0.9722100, 0.9801946, 0.9859050, 0.9899806, 0.9928846, 0.9949508, 1 0.9694225, 0.9781908, 0.9844689, 0.9889539, 0.9921520, 0.9944289, 1 0.9663700, 0.9759930, 0.9828916, 0.9878250, 0.9913457, 0.9938541, 1 0.9630303, 0.9735839, 0.9811602, 0.9865842, 0.9904587, 0.9932213, 1 0.9593799, 0.9709456, 0.9792609, 0.9852213, 0.9894833, 0.9925248, 1 0.9553940, 0.9680586, 0.9771789, 0.9837252, 0.9884113, 0.9917585, 1 0.9510467, 0.9649025, 0.9748985, 0.9820839, 0.9872338, 0.9909158, 1 0.9463112, 0.9614558, 0.9724028, 0.9802845, 0.9859409, 0.9899896, 1 0.9411597, 0.9576958, 0.9696741, 0.9783134, 0.9845225, 0.9889721, 1 0.9355637, 0.9535990, 0.9666935, 0.9761560, 0.9829674, 0.9878550, 1 0.9294942, 0.9491410, 0.9634414, 0.9737967, 0.9812637, 0.9866292, 1 0.9229221, 0.9442967, 0.9598970, 0.9712192, 0.9793986, 0.9852851, 1 0.9158184, 0.9390404, 0.9560391, 0.9684063, 0.9773587, 0.9838125, 1 0.9081548, 0.9333462, 0.9518454, 0.9653399, 0.9751298, 0.9822001/ data ((fx(ii,ij),ij=1,6),ii=30,41)/ 1 0.8999038, 0.9271882, 0.9472933, 0.9620012, 0.9726967, 0.9804364, 1 0.8910396, 0.9205407, 0.9423598, 0.9583707, 0.9700437, 0.9785088, 1 0.8815383, 0.9133786, 0.9370216, 0.9544284, 0.9671543, 0.9764043, 1 0.8713786, 0.9056778, 0.9312555, 0.9501539, 0.9640114, 0.9741089, 1 0.8605424, 0.8974157, 0.9250387, 0.9455263, 0.9605972, 0.9716083, 1 0.8490153, 0.8885715, 0.9183489, 0.9405248, 0.9568935, 0.9688874, 1 0.8367872, 0.8791267, 0.9111650, 0.9351287, 0.9528820, 0.9659305, 1 0.8238531, 0.8690656, 0.9034670, 0.9293177, 0.9485439, 0.9627216, 1 0.8102129, 0.8583758, 0.8952365, 0.9230718, 0.9438604, 0.9592442, 1 0.7958728, 0.8470484, 0.8864572, 0.9163723, 0.9388129, 0.9554815, 1 0.7808969, 0.8351 , 0.8771039, 0.9092181, 0.9333918, 0.9514126, 1 0.7651953, 0.8225 , 0.8671884, 0.9015595, 0.9275623, 0.9470281/ data ((fx(ii,ij),ij=7,7),ii=1,41)/ 1 0.9997143, 0.9996843, 0.9996511, 0.9996145, 0.9995741, 0.9995294, 1 0.9994801, 0.9994256, 0.9993654, 0.9992989, 0.9992254, 0.9991444, 1 0.9990549, 0.9989560, 0.9988469, 0.9987265, 0.9985936, 0.9984470, 1 0.9982852, 0.9981068, 0.9979099, 0.9976929, 0.9974536, 0.9971899, 1 0.9968993, 0.9965793, 0.9962269, 0.9958390, 0.9954121, 0.9949426, 1 0.9944264, 0.9938592, 0.9932363, 0.9925526, 0.9918025, 0.9909804, 1 0.9900797, 0.9890938, 0.9880156, 0.9868374, 0.9855511/ c c..calculate spline intepolation coefficients, gamma functions, c..and series constants, if this is the first pass. if (it) 10,10,144 10 it=1 itab=1 gam(7)=sqrt(pi) gam(1)=0.5*gam(7) gam(2)=1. gam(3)=1.5*gam(1) gam(4)=2. gam(5)=2.5*gam(3) gam(6)=6. afac(7)=1./sqrt(2.d0) afac(1)=.5*afac(7) do 20 i1=2,6 20 afac(i1)=afac(7)*afac(i1-1) lf=0 do 11 i1=1,7 do 12 i2=1,67 eqfmbe(lf+i2)=eqfm(lf+i2) 12 continue zero=0. call solve(etai,eqfmbe(lf+1),67,0,0,zero,zero) lf=lf+67 11 continue lf=0 do 14 i1=1,7 do 15 i2=1,41 eqfxbe(lf+i2)=eqfx(lf+i2) 15 continue zero=0. call solve(etax,eqfxbe(lf+1),41,0,0,zero,zero) lf=lf+41 14 continue 144 continue c..check range of eta if(eta+3.8) 31,30,30 30 if(eta+.14) 41,40,40 40 if (eta-18.) 50,50,60 c..use asymptotic series to calculate fermi integrals c..if eta lt -3.8 c..see cox and giuli equ. 24.277. 31 eeta=exp(eta) do 17 l=1,7 fe(l)=gam(l)*eeta*(1.-afac(l)*eeta) 17 continue go to 70 c..calculate fermi integrals and derivatives using spline c..interpolated table for eta between -3.8 and -0.14 41 lf=0 eeta=exp(eta) call search(eta,itab,etax,41) do 16 i3=1,ixne if(i3.eq.2) go to 80 if(i3.eq.7)go to 16 call seval2(eta,fe(i3),febe(i3),etax,eqfx(lf+1), 1 eqfxbe(lf+1),41,itab) xk=dble(i3)*.5 fbar=gam(i3)*eeta febe(i3)=(febe(i3)+fe(i3))*fbar fe(i3)=fe(i3)*fbar 80 lf=lf+41 16 continue go to 70 c..calculate fermi integrals and derivatives using spline interpolated c..table for eta between -0.14 and 18. 50 lf=0 call search(eta,itab,etai,67) do 13 i3=1,ixne if(i3.eq.2) go to 90 if(i3.eq.7) go to 13 call seval2(eta,fe(i3),febe(i3),etai,eqfm(lf+1),eqfmbe(lf+1), 1 67,itab) xk=dble(i3)*.5 fbar=((eta+.16)**(xk+1.))/(xk+1.)+xfe(i3) febe(i3)=febe(i3)*fbar+fe(i3)*(eta+0.16)**xk fe(i3)=fe(i3)*fbar 90 lf=lf+67 13 continue go to 70 c..eta gt 18. use asymptotic series(eq.24.38 and 24.39 of cox giuli ) 60 etar=sqrt(eta) eta2=eta*eta eta2i=1./eta2 fe1=c10*eta*etar fe3=0.6*eta*fe1 fe(1)=fe1*(1.+eta2i*(c12+eta2i*(c14+eta2i*c16))) fe(2)=0.5*eta2+c22 fe(3)=fe3*(1.+eta2i*(c32+eta2i*(c34+eta2i*c36))) fe(4)=eta*(eta2/3.+c42) fe(7)=2.*etar*(1.+eta2i*(c72+eta2i*(c74+eta2i*c76))) if(ixne.eq.4) go to 70 fe(5)=c50*eta*fe3*(1.+eta2i*(c52+eta2i*(c54*eta2i*c56))) fe(6)=0.25*eta2*eta2+c62*eta2+c64 c..use nonrelativistic fermi functions to calculate state variables c.. for arbitrary degeneracy and for relativistic case. c..use approximation of divine,n., ap. j. 142, 1652 (1965). see also c.. chiu page 131. 70 beta=cbeta*t betar=sqrt(beta) beta32=beta*betar if(eta.le.18..and.eta.ge.-3.8) go to 100 febe(1)=0.5*fe(7) febe(3)=1.5*fe(1) febe(4)=2.*fe(2) febe(5)=2.5*fe(3) febe(6)=3.*fe(4) c.. 100 u1=(fe(4)/fe(3))**2 u1i=1./u1 funu1r=sqrt(1.+0.5*u1*beta) funu1l=1.+u1*beta ffunu1=funu1l*funu1r xne = cxne * ( fe(1) + fe(3)*(ffunu1 - 1.)*u1i)*beta32 if(ixne.eq.4) go to 110 betai=1./beta beta52=beta32*beta u3=(fe(6)/fe(5))**2 u3i=1./u3 u3i2=u3i*u3i funu3r=sqrt(1.+0.5*u3*beta) funu3l=1.+u3*beta ffunu3=funu3l*funu3r gfunu3=funu3r**3 ee = cee * ( fe(3) + fe(5)*(ffunu3 - 1.)*u3i)*beta52 pe = cpe * ( fe(3) + fe(5)*(gfunu3 - 1.)*u3i)*beta52 c..calculate derivatives: xnebe,xnebt,eee,eet,pee,pet betabt=cbeta 110 u1beta=2.*u1*(febe(4)/fe(4)-febe(3)/fe(3)) f1bbeta=u1*(funu1r+funu1l/(4.*funu1r)) f1beta=u1beta*beta*f1bbeta*u1i x=ffunu1-1. xnebe=cxne*beta32*(febe(1)+febe(3)*x*u1i+fe(3)*(f1beta*u1i 1 -x*u1i*u1i*u1beta)) if(ixne.eq.4) go to 200 xnebt=betabt*(1.5*betai*xne+cxne*beta32*fe(3)*f1bbeta*u1i) u3beta=2.*u3*(febe(6)/fe(6)-febe(5)/fe(5)) f3bbeta=u3*(funu3r+funu3l/(4.*funu3r)) f3beta=u3beta*beta*f3bbeta*u3i g3beta=0.75*funu3r*beta*u3beta g3bbeta=0.75*funu3r*u3 x=ffunu3-1. eee=cee*beta52*(febe(3)+febe(5)*x*u3i+fe(5)*(f3beta*u3i 1 -x*u3i2*u3beta)) eet=betabt*(2.5*betai*ee+cee*beta52*fe(5)*f3bbeta*u3i) x=gfunu3-1. pee=cpe*beta52*(febe(3)+febe(5)*x*u3i+fe(5)*(g3beta*u3i 1 -x*u3i2*u3beta)) pet=betabt*(2.5*betai*pe+cpe*beta52*fe(5)*g3bbeta*u3i) 200 return end subroutine solve(x,y,n,if1,if2,d1,d2) implicit none save c..this routine calculates the derivatives needed for 1d spline interpolation. c.. x = array of independent variable values. c.. y = array of dependent variable values. c.. n = number of values in the x and y arrays. c.. the derivatives of y with respect to x are returned in array y. c.. c..declare integer n,if1,if2,iflag,i,k,nm2 double precision y(1),x(1),d1,d2,xt(4),yt(4),dd1,dd2,funda, + cl,cm c.. c..communicate double precision al(67),ad(67),au(67),b(67),pq(134) common /scom1/ al,ad,au,b,pq c.. c..initialize iflag=0 dd1=d1 if((if1.eq.1).or.(n.lt.4)) go to 1 dd1=funda(x,y) 1 continue dd2=d2 if((if2.eq.1).or.(n.lt.4)) go to 3 do 2 i=1,4 k=n+1-i xt(i)=x(k) yt(i)=y(k) 2 continue dd2=funda(xt,yt) 3 continue c.. c..set up the matrix coefficients and solve nm2=n-2 do 4 i=1,nm2 cl=(x(i+2)-x(i+1))/(x(i+2)-x(i)) cm=1.-cl al(i)=cl ad(i)=2. au(i)=cm b(i)=3.*cl*(y(i+1)-y(i))/(x(i+1)-x(i))+ + 3.*cm*(y(i+2)-y(i+1))/(x(i+2)-x(i+1)) 4 continue b(1)=b(1)-al(1)*dd1 b(nm2)=b(nm2)-au(nm2)*dd2 call trislv(al,ad,au,nm2,b,pq,y(2),iflag) if(iflag.ne.0) stop 'iflag = 0 in routine solve' c.. y(1)=dd1 y(n)=dd2 return end double precision function funda(x,y) implicit none save c.. c..looks like a matrix determinant c.. double precision x(1),y(1),h1,h2,h3,h12,h123,h23,t1,t2,t3 h1=x(2)-x(1) h2=x(3)-x(2) h3=x(4)-x(3) h12=h1+h2 h123=h12+h3 h23=h2+h3 t1=((y(2)-y(1))*h12*h123)/(h1*h2*h23) t2=((y(3)-y(1))*h123*h1)/(h12*h2*h3) t3=((y(4)-y(1))*h1*h12)/(h123*h3*h23) funda=t1-t2+t3 return end subroutine trislv (al,ad,au,n,b,pq,x,iflag) implicit none save c..solution of the tridiagonal linear system a*x = b by gaussian elimination c..without pivoting. input is al of dimension n, the lower diagonal of a (note c..al(1) is arbitrary); ad also of dimension n, the main diagonal of a; au, c..the upper diagonal of a also of dimension n (note au(n) is arbitrary); c..b, the the right-hand side vector and pq, a working vector of length at c..least 2*n. on output x is filled with the solution vector and iflag is c..set to 0 for normal return and 1 if the matrix is singular. c.. c..written by: f. n. fritsch, lawrence livermore laboratory. c..date last changed: 13 march 1974. c..modified by t.a. weaver on 5/1/76 to give correct results for n=2 c.. c.. c..notes on memory management: c.. pq(i+1) = p(i), i = o(1)n-1 c.. pq(n+i) = q(i), i = 1(1)n c.. c..declare integer n,iflag,i,npi,k double precision al(1),ad(1),au(1),b(1),x(1),pq(1),denom c.. c..stop silliness if (n.lt.2) return c.. c..for n greater than 2 if (n .gt. 2) then pq(1) = 0. pq(n) = 0. do 10 i=1,n-1 denom = al(i)*pq(i) + ad(i) if (denom .eq. 0.0) then iflag = i return end if pq(i+1) = -au(i)/denom npi = n+i pq(npi) = ( b(i) - al(i)*pq(npi-1) )/denom 10 continue denom = al(n)*pq(n) + ad(n) if (denom .eq. 0.0) then i = n iflag = i return end if pq(npi+1) = ( b(n) - al(n)*pq(npi) )/denom c.. c..back substitution x(n) = pq(2*n) do 20 k = 2, n i = n-k+1 npi = n+i x(i) = pq(npi) + pq(i+1)*x(i+1) 20 continue iflag = 0 return c.. c..for n=2 else if (n .eq. 2) then denom=ad(1)*ad(2)-al(2)*au(1) i=222 if (denom .eq. 0.0) then iflag = i return end if x(1)=(ad(2)*b(1)-au(1)*b(2))/denom x(2)=(ad(1)*b(2)-al(2)*b(1))/denom iflag=0 return end if end subroutine seval2(yval,x,dx,rv,hv,fyh,npts,ily) implicit none c..this subroutine calculates the spline interpolated value x, and derivative c..dx at the point yval given the given the array of independent variables rv, c..the array of dependent variables hv and the functions derivatives fyh. npts c..is the number of points in the arrays. c.. c..rv must increase monotonically and yval must be in the range of rv. c..the derivative array fyh must be from the routine solve. c..the first table entry ily, less than yval must also be supplied. c.. c..declare integer npts,ily,ily1 double precision yval,x,dx,rv(1),hv(1),fyh(1),dy,hy,hyi,r,r2, + fh2,fh3,fh4,fh1d,fh3d,fh4d c.. if (ily.eq.0) call search(yval,ily,rv,npts) ily1=ily+1 dy=yval-rv(ily) hy=rv(ily1)-rv(ily) hyi=1./hy r=dy*hyi r2=r*r fh2=-2.*r2*r+3.*r2 fh3=dy*(r2-2.*r+1.) fh4=dy*(r2-r) fh1d=6.*(r2-r)*hyi fh3d=3.*r2-4.*r+1. fh4d=3.*r2-2.*r x=hv(ily)*(1.-fh2)+hv(ily1)*fh2+fyh(ily)*fh3+fyh(ily1)*fh4 dx=(hv(ily)-hv(ily1))*fh1d+fyh(ily)*fh3d+fyh(ily1)*fh4d return end subroutine search(xval,i,x,n) implicit none save c..this routine calculates the table entry i which falls just below xval in c..the monotonically increasing array x(i) which is n points long. c.. integer n,i double precision xval,x(1) c.. if((xval.lt.x(1)).or.(xval.gt.x(n))) stop 'bad xval in search' if((i.le.0).or.(i.gt.n)) i=n/2+1 if(xval-x(i)) 10,20,30 10 i=i-1 if(xval-x(i)) 10,20,40 30 i=i+1 if(xval-x(i))20,20,30 20 i=i-1 if(i.eq.0)i=1 40 return end double precision function funcne(eta,funcnbeta) implicit none save c.. c..funcne returns the normalized value of the difference between the required c..value of xnem (xncom) and its value as a function of eta (xnemc). c..called by rootnr in es c.. c..declare double precision eta,funcnbeta,xnp,xnpbe,xne,xnebe,beta,zer, + xnemc,denom,etap,c,cbeta parameter (c=2.997925e+10, + cbeta=1.38054e-16/(9.10908e-28*c*c)) c.. c..communication with funcne the eta root find routine integer npfla466 common /compari/ npfla466 double precision xncom,tcom common /comparr/ xncom,tcom c.. data xnp /0.0/ data xnpbe /0.0/ data xnebe /0.0/ c.. beta=cbeta*tcom xne=0.0 zer=0.0 call emeos(tcom,eta,4,xne,zer,zer,xnebe,zer,zer,zer,zer,zer) xnemc=xne if((beta.le.0.02).or.(npfla466.le.0)) go to 10 etap=max(-eta-2.d0/beta,-2.d+02) call emeos(tcom,etap,4,xnp,zer,zer,xnpbe,zer,zer,zer,zer,zer) xnemc=xne-xnp if(etap.gt.-200.) xnebe=xnebe+xnpbe 10 denom=1./(xncom+20.3*tcom*tcom*tcom) funcne=(xnemc-xncom)*denom funcnbeta=xnebe*denom c write(6,112) xne,xnp,xncom,xne-xncom,eta c 112 format(1x,1p6e12.4) return end double precision function rootnr(func,gess,xconv,j) implicit none save c..this routine uses newton raphson to find the root of the function, func, c..to a relative accuracy of xconv, given an initial guess, gess,for the root. c..if convergence is not reached in j cycles c..j is changed to minus j (initial) c..otherwise the returned j is the number of iterations performed. c.. c..func must be a function of two arguments, the first being the independent c..variable, and the second being the derivative of func with respect to the c..indep variable. c.. c..declare external func integer j,jroot double precision func,gess,xconv,x(50),y(50),dy(50),dydx,dx data dydx /0.0/ c.. c.. jroot=j x(1)=gess do 10 j=1,jroot y(j)=func(x(j),dydx) dy(j)=dydx dx=-y(j)/(dydx+1.d-199) x(j+1)=x(j)+dx if(abs(dx)-abs(x(1)*xconv)) 20,20,10 10 continue rootnr=x(jroot+1) j=-jroot return 20 rootnr=x(j+1) return end