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