program derfer_main c This program calls dfrmintnew (calcdfi) to calculate the c thermodynamics variables and their derivatives c (dP/dT)_rho, (du/dT)_rho and (dP/drho)_T c for the fermi gas implicit double precision (a-h,o-z) dimension ff(3), dff(6), ffth(3), dffth(6) integer igo, ivar, eval real*8 lbstart, lbend, lbinc integer interp_order real*8 beta, betalog, bmass, cboltz, clight, \$ dedr_t, dedt_r, dedt_r1, dedt_r2, dedt_r3, \$ dedt_r4, dff, dffth, dpdr_t, \$ dpdt_r, dpdt_r1, dpdt_r2, dpdt_r3, dpdt_r4, \$ e, emass, en, eta, etaend, \$ etainc, etastart, ff, ffth, hplanck, \$ p, pi, x1, x2, x3, \$ x4, x5 common /evalt/ eval, interp_order interp_order = 6 emass = 9.11 ! d - 28 gr pi = 2.*asin(1.) cboltz = 1.381 ! d - 16 erg deg - 1 bmass = 1.673 !d - 24 gr hplanck = 6.626 !d - 27 erg sec clight = 2.9979 !d10 cm sec - 1 igo = 1 10 continue if (igo.ne.2) then write (6,'(a,\$)')' 0 stop, 1 constant beta , 2 constant eta ?' read (5,*) ivar if (ivar.eq.0) then stop elseif (ivar.eq.1) then write (6,'(a)') \$ ' LOG 10 beta value, eta-start, eta-end, d eta ? ' read (5,*) lbstart, etastart, etaend, etainc lbend = lbstart lbinc = 1 elseif (ivar.eq.2) then write (6,'(a)') ' Eta value, log beta start, log beta end, \$ log beta inc ? ' read (5,*) etastart, lbstart, lbend, lbinc etaend = etastart etainc = 1 endif endif write (6,'(2a)')' Eta logBeta n P E ', \$ ' dP/dT dE/dT dP/dn dE/dn X4 X5' do eta = etastart, etaend, etainc do betalog = lbstart, lbend, lbinc beta = 10.**betalog call calcdfi (eta, beta, \$ ff, \$ dff, \$ ffth, \$ dffth) p = 16.*pi*sqrt(2.)*emass**4*clight**5/3./hplanck**3 p = p*beta**(5./2.)*(ff(2)+(1./2.)*beta*ff(3)) p = p*1.d19 !cgs if (eta.gt.70.)then dpdt_r1 = (5.*ff(2)/2.+beta*dff(5)+7.*beta*ff(3)/4.+ \$ beta*beta*dff(6)/2.)*(dffth(1) + \$ beta*dffth(2)) dpdt_r2 = -(3.*ff(1)/2.+beta*dff(4)+5.*beta*ff(2)/2.+ \$ beta*beta*dff(5))*(dffth(2) + \$ beta*dffth(3)/2.) dpdt_r3 = (dff(1)+beta*dff(2)-dffth(1)-beta*dffth(2))* \$ (5.*ffth(2)/2.+beta*dffth(5)+7.*beta \$ *ffth(3)/4.+ \$ beta*beta*dffth(6)/2.) dpdt_r4 = -(dff(2)+beta*dff(3)/2.-dffth(2)-beta* \$ dffth(3)/2.)*(3.*ffth(1)/2. + \$ beta*dffth(4) + \$ 5.*beta*ffth(2)/2. + beta*beta*dffth(5)) dpdt_r = (dpdt_r1+dpdt_r2+dpdt_r3+dpdt_r4)/(dff(1)+ \$ beta*dff(2)) else dpdt_r = (5./2.)*ff(2) + beta*dff(5) + \$ (7./4.)*beta*ff(3) + \$ (1./2.)*beta*beta*dff(6) dpdt_r = dpdt_r - (dff(2)+(1./2.)*beta*dff(3))* \$ ((3./2.)*ff(1)+beta*dff(4)+(5./2.)*beta \$ *ff(2)+ \$ beta*beta*dff(5))/(dff(1) + beta*dff(2)) endif dpdt_r = 16.*pi*sqrt(2.)*emass**3*clight**3/3./hplanck**3* \$ cboltz*dpdt_r*beta**(3./2.) dpdt_r = dpdt_r*1.d11 ! cgs e = 8.*pi*sqrt(2.)*emass**4*clight**5/hplanck**3 e = e*beta**(5./2.)*(ff(2)+beta*ff(3)) e = e*1.d19 !cgs if (eta.gt.70.)then dedt_r1 = (5.*ff(2)/2.+beta*dff(5)+7.*beta*ff(3)/2.+ \$ beta*beta*dff(6))*(dffth(1) + beta*dffth(2)) dedt_r2 = -(3.*ff(1)/2.+beta*dff(4)+5.*beta*ff(2)/2.+ \$ beta*beta*dff(5))*(dffth(2) + beta*dffth(3)) dedt_r3 = (dff(1)+beta*dff(2)-dffth(1)-beta*dffth(2))* \$ (5.*ffth(2)/2.+beta*dffth(5)+7.*beta \$ *ffth(3)/2.+beta*beta*dffth(6)) dedt_r4 = -(dff(2)+beta*dff(3)-dffth(2)-beta \$ *dffth(3))* \$ (3.*ffth(1)/2.+beta*dffth(4)+5.*beta \$ *ffth(2)/2.+beta*beta*dffth(5)) dedt_r = (dedt_r1+dedt_r2+dedt_r3+dedt_r4)/(dff(1)+ \$ beta*dff(2)) else dedt_r = (5./2.)*ff(2) + beta*dff(5) + \$ (7./2.)*beta*ff(3) + \$ beta*beta*dff(6) dedt_r = dedt_r - (dff(2)+beta*dff(3))* \$ ((3./2.)*ff(1)+beta*dff(4)+(5./2.)*beta \$ *ff(2)+ \$ beta*beta*dff(5))/(dff(1) + beta*dff(2)) endif dedt_r = 8.*pi*sqrt(2.)*emass**3*clight**3/hplanck**3* \$ cboltz*dedt_r*beta**(3./2.) dedt_r = dedt_r*1.d11 !cgs en = 8.*pi*sqrt(2.)/hplanck**3*emass**3*clight**3 en = en*beta**(3./2.)*(ff(1)+beta*ff(2)) en = en*1.d27 !cgs dpdr_t = (2./3.)*emass*clight*clight* \$ beta*(dff(2)+(1./2.)*beta*dff(3))/ \$ (dff(1)+beta*dff(2)) dpdr_t = dpdr_t*1.d-8 !cgs dedr_t = emass*clight*clight* \$ beta*(dff(2)+beta*dff(3))/ \$ (dff(1)+beta*dff(2)) dedr_t = dedr_t*1.d-8 !cgs x1 = beta*emass*clight**2/cboltz*dpdt_r/p c x2 = beta*emass*clight**2/cboltz*dedt_r/e x3 = en*dpdr_t/p x4 = dedt_r/dpdt_r x5 = (x1+x4*x3)/(x1+x4*x3-1.) write (6,'(2f8.3,9(1pe10.3))') eta,betalog,en,p,e, \$ dpdr_t,dedr_t,dpdt_r,dedt_r, \$ x4,x5 enddo enddo write (6,'(a)') ' ' go to 10 end