      Program STATSTAR
c
c  This program is listed in Appendix H of both "An Introduction
c  to Modern Astrophysics," Bradley W. Carroll and Dale A. Ostlie,
c  Addison-Wesley Publishing Company, copyright 1996, and "An 
c  Introduction to Modern Stellar Astrophysics," Dale A. Ostlie
c  and Bradley W. Carroll, Addison-Wesley Publishing Company,
c  copyright 1996.
c
c
c  This program will calculate a static stellar model using the 
c  equations developed in the text.  The user is expected to supply the
c  star's mass, luminosity, effective temperature, and composition 
c  (X and Z).  If the choices for these quantities are not consistent 
c  with the central boundary conditions, an error message will be 
c  generated and the user will then need to supply a different set of
c  initial values.
c
      real*8 r(999), P(999), M_r(999), L_r(999), T(999), deltar, Te
      real*8 rho(999), kappa(999), epslon(999), dlPdlT(999)
      real*8 X, Y, Z, XCNO, mu
      real*8 Ms, Ls, Rs, T0, P0
      real*8 Pcore, Tcore, rhocor, epscor, rhomax
      real*8 Rsolar, Msolar, Lsolar, Qm, Rcrat, Mcrat, Lcrat
      real*8 deltam, dlPlim
      real*8 Rsun, Msun, Lsun
      real*8 f_im1(4), f_i(4), dfdr(4)
      real*8 dMdr, dPdr, dLdr, dTdr
      real*8 sigma, c, a, G, k_B, m_H, pi, gamma, gamrat, kPad,  tog_bf,
     1 g_ff
      character clim, rcf
      common /cnstnt/ sigma, c, a, G, k_B, m_H, pi, gamma, gamrat, kPad,
     1 g_ff
c
c  deltar = radius integration step
c  idrflg = set size flag
c         = 0 (initial surface step size of Rs/1000.)
c         = 1 (standard step size of Rs/100.)
c         = 2 (core step size of Rs/5000.)
c  
c  Nstart = number of steps for which starting equations are to be used
c           (the outermost zone is assumed to be radiative)
c  Nstop = maximum number of allowed zones in the star
c  Igoof = final model condition flag
c        = -1 (number of zones exceeded; also the initial value)
c        =  0 (good model)
c        =  1 (core density was extreme)
c        =  2 (core luminosity was extreme)
c        =  3 (extrapolated core temperature is too low)
c        =  4 (mass became negative before center was reached)
c        =  5 (luminosity became negative before center was reached)
c  X, Y, Z = mass fractions of hydrogen, helium, and metals
c  T0, P0 = surface temperature and pressure (T0 = P0 = 0 is assumed)
c  Ms, Ls, Rs = mass, luminosity, and radius of the star (cgs units)
c
      data Nstart, Nstop, Igoof, ierr / 10, 999, -1, 0 /                    
      data P0, T0, dlPlim / 0.0d0, 0.0d0, 99.9d0 /
      data Rsun, Msun, Lsun / 6.9599d+10, 1.989d+33, 3.826d+33 /
      data sigma, c, a, G, k_B, m_H, pi, gamma, tog_bf, g_ff
     1 / 5.67051d-5, 2.99792458d+10, 7.56591d-15, 6.67259d-8,
     2   1.380658d-16, 1.673534d-24, 3.141592654d0, 1.6666667d0, 0.01d0,
     3   1.0d0 /
c
c  Assign values to constants (cgs units)
c     Rsun = radius of the Sun
c     Msun = mass of the Sun
c     Lsun = luminosity of the Sun
c     sigma = Stefan-Boltzmann constant
c     c = speed of light in vacuum
c     a = 4*sigma/c (radiation pressure constant)
c     G = universal gravitational constant
c     K_B = Boltzmann constant
c     m_H = mass of hydrogen atom
c     pi = 3.141592654
c     gamma = 5/3 (adiabatic gamma for a monatomic gas)
c     gamrat = gamma/(gamma-1)
c     kPad = P/T**(gamma/(gamma-1)) (adiabatic gas law constant)
c     tog_bf = bound-free opacity constant 
c             (ratio of guillotine to gaunt factors)
c     g_ff = free-free opacity gaunt factor (assumed to be unity)
c
c      open(unit=20,file='starmodl.dat',form='formatted',
c     1 status='unknown')
c                             
c  Enter desired stellar parameters
c
      write(*,*) ' Enter the mass of the star (in solar units):'
      read(*,*) Msolar
      write(*,*) ' Enter the luminosity of the star (in solar units):'
      read(*,*) Lsolar
      write(*,*) ' Enter the effective temperature of the star (in K):'
      read(*,*) Te
   10 continue
          write(*,*) ' Enter the mass fraction of hydrogen (X):'
          read(*,*) X
          write(*,*) ' Enter the mass fraction of metals (Z):'
          read(*,*) Z
          Y = 1.d0 - X - Z
          if (Y.lt.0.0d0) then
              write(*,100) 
              go to 10
          end if
c
c  Select the mass fraction CNO to be 50% of Z.
c
      XCNO = Z/2.0d0
c     
c  Calculate the mass, luminosity, and radius of the star.
c  The radius is calculated from Eq. (3.17).
c
      Ms = Msolar*Msun
      Ls = Lsolar*Lsun
      Rs = sqrt(Ls/(4.d0*pi*sigma))/Te**2
      Rsolar = Rs/Rsun
c
c  Begin with a very small step size since surface conditions vary
c  rapidly.
c
      deltar = -Rs/1000.0d0
      idrflg = 0
c     
c  Calculate mean molecular weight mu assuming complete ionization
c  (see Eq. 10.21).
c
      mu = 1.0d0/(2.0d0*X + 0.75d0*Y + 0.5d0*Z)
c
c  Calculate the delimiter between adiabatic convection and radiation
c  (see Eq. 10.87).
c
      gamrat = gamma/(gamma - 1.0d0)
c
c  Initialize values of r, P, M_r, L_r, T, rho, kappa, and epslon at
c  the surface.  The outermost zone is assumed to be zone 1.  The zone 
c  number increases toward the center.
c
      r(1)   = Rs
      M_r(1) = Ms
      L_r(1) = Ls
      T(1)   = T0
      P(1)   = P0
      if (P0.le.0.0d0 .or. T0.le.0.0d0) then
          rho(1)    = 0.0d0
          kappa(1)  = 0.0d0
          epslon(1) = 0.0d0
      else
          call EOS(X, Z, XCNO, mu, P(1), T(1), rho(1), kappa(1),
     1        epslon(1), tog_bf, 1 ,ierr)
          if (ierr.ne.0) stop
      end if
c
c  Apply approximate surface solutions to begin the integration,
c  assuming radiation transport in the outermost zone (the do 20 loop).
c  irc = 0 for radiation, irc = 1 for convection.
c  Assume arbitrary initial values for kPad, and dlPdlT.
c  dlPdlT = dlnP/dlnT (see Eq. 10.87)
c
      kPad = 0.3d0
      irc = 0
      dlPdlT(1) = 4.25d0
      do 20 i = 1, Nstart
          ip1 = i + 1
          call STARTMDL(deltar, X, Z, mu, Rs, r(i), M_r(i), L_r(i), 
     1        r(ip1), P(ip1), M_r(ip1), L_r(ip1), T(ip1), tog_bf, irc)
          call EOS(X, Z, XCNO, mu, P(ip1), T(ip1), rho(ip1), kappa(ip1),
     1        epslon(ip1), tog_bf, ip1, ierr)
          if (ierr.ne.0) then 
              write(*,400) r(i)/Rs, rho(i), M_r(i)/Ms, kappa(i), T(i), 
     1            epslon(i), P(i), L_r(i)/Ls
              stop
          end if
c
c  Determine whether convection will be operating in the next zone by
c  calculating dlnP/dlnT numerically between zones i and i+1 (ip1).
c  Update the adiabatic gas constant if necessary.
c
          if (i.gt.1) then
              dlPdlT(ip1) = log(P(ip1)/P(i))/log(T(ip1)/T(i))
          else
              dlPdlT(ip1) = dlPdlT(i)
          end if
          if (dlPdlT(ip1).lt.gamrat) then
              irc = 1
          else
              irc = 0
              kPad = P(ip1)/T(ip1)**gamrat
          end if
c
c  Test to see whether the surface assumption of constant mass is still
c  valid.
c
          deltaM = deltar*dMdr(r(ip1), rho(ip1))
          M_r(ip1) = M_r(i) + deltaM
          if (abs(deltaM).gt.0.001d0*Ms) then
              write(*,200)
              if (ip1.gt.2) ip1 = ip1 - 1
              go to 30
          end if
   20 continue
c
c  This is the main integration loop.  The assumptions of constant
c  interior mass and luminosity are no longer applied.
c
   30 Nsrtp1 = ip1 + 1
      do 40 i = Nsrtp1, Nstop
          im1 = i - 1
c
c  Initialize the Runge-Kutta routine by specifying zone i-1 quantities
c  and their derivatives.  Note that the pressure, mass, luminosity,
c  and temperature are stored in the memory locations f_im1(1),
c  f_im1(2), f_im1(3), and f_im1(4), respectively.  The derivatives of
c  those quantities with respect to radius are stored in dfdr(1),
c  dfdr(2), dfdr(3), and dfdr(4).  Finally, the resulting values for
c  P, M_r, L_r, and T are returned from the Runge-Kutta routine in
c  f_i(1), f_i(2), f_i(3), and f_i(4), respectively.
c
c  The stellar structure equations dPdr (Eq. 10.7), dMdr (Eq. 10.8),
c  dLdr (Eq. 10.45), and dTdr (Eq. 10.61 or Eq. 10.81) are calculated 
c  in function calls, defined later in the code.
c
          f_im1(1) = P(im1)
          f_im1(2) = M_r(im1)
          f_im1(3) = L_r(im1)
          f_im1(4) = T(im1)
          dfdr(1)  = dPdr(r(im1), M_r(im1), rho(im1))
          dfdr(2)  = dMdr(r(im1), rho(im1))
          dfdr(3)  = dLdr(r(im1), rho(im1), epslon(im1))
          dfdr(4)  = dTdr(r(im1), M_r(im1), L_r(im1), T(im1), rho(im1),
     1        kappa(im1), mu, irc)
          call RUNGE(f_im1, dfdr, f_i, r(im1), deltar, irc, X, Z, XCNO,
     1        mu, i, ierr)
          if (ierr.ne.0) then           
              write(*,300)
              write(*,400) r(im1)/Rs, rho(im1), M_r(im1)/Ms, kappa(im1), 
     1            T(im1), epslon(im1), P(im1), L_r(im1)/Ls
              stop
          end if
c     
c  Update stellar parameters for the next zone, including adding
c  dr to the old radius (note that dr <  0 since the integration is
c  inward).
c
          r(i)   = r(im1) + deltar
          P(i)   = f_i(1)
          M_r(i) = f_i(2)
          L_r(i) = f_i(3)
          T(i)   = f_i(4)
c
c  Calculate the density, opacity, and energy generation rate for
c  this zone.
c
          call EOS(X, Z, XCNO, mu, P(i), T(i), rho(i), kappa(i),
     1        epslon(i), tog_bf, i, ierr)
          if (ierr.ne.0) then 
              write(*,400) r(im1)/Rs, rho(im1), M_r(im1)/Ms, kappa(im1),
     1            T(im1), epslon(im1), P(im1), L_r(im1)/Ls
              stop
          end if
c
c  Determine whether convection will be operating in the next zone by
c  calculating dlnP/dlnT and comparing it to gamma/(gamma-1)
c  (see Eq. 10.87).  Set the convection flag appropriately.
c
          dlPdlT(i) = log(P(i)/P(im1))/log(T(i)/T(im1))
          if (dlPdlT(i).lt.gamrat) then
              irc = 1
          else
              irc = 0
          end if
c
c  Check to see whether the center has been reached.  If so, set Igoof and 
c  estimate the central conditions rhocor, epscor, Pcore, and Tcore.
c  The central density is estimated to be the average density of the
c  remaining central ball, the central pressure is determined by
c  applying Eq. (H.4), and the central value for the energy
c  generation rate is calculated to be the remaining interior
c  luminosity divided by the mass of the central ball.  Finally, the
c  central temperature is computed by applying the ideal gas law
c  (where radiation pressure is neglected).
c
          if (r(i).le.abs(deltar) .and.
     1      (L_r(i).ge.0.1d0*Ls .or. M_r(i).ge.0.01d0*Ms)) then
              Igoof = 6
          else if (L_r(i).le.0.0d0) then
              Igoof = 5
              rhocor = M_r(i)/(4.0d0/3.0d0*pi*r(i)**3)
              if (M_r(i).ne.0.0d0) then
                  epscor = L_r(i)/M_r(i)
              else
                  epscor = 0.0d0
              end if
              Pcore = P(i) + 2.0d0/3.0d0*pi*G*rhocor**2*r(i)**2
              Tcore = Pcore*mu*m_H/(rhocor*k_B)
          else if (M_r(i).le.0.0d0) then
              Igoof  = 4
              Rhocor = 0.0d0
              epscor = 0.0d0
              Pcore  = 0.0d0
              Tcore  = 0.0d0
          else if (r(i).lt.0.02d0*Rs .and. M_r(i).lt.0.01d0*Ms .and.
     1      L_r(i).lt.0.1d0*Ls) then
              rhocor = m_r(i)/(4./3.*pi*r(i)**3)
              rhomax = 10.0d0*(rho(i)/rho(im1))*rho(i)
              epscor = L_r(i)/M_r(i)
              Pcore  = P(i) + 2.0d0/3.0d0*pi*G*rhocor**2*r(i)**2
              Tcore  = Pcore*mu*m_H/(rhocor*k_B)
              if (rhocor.lt.rho(i) .or. rhocor.gt.rhomax) then
                  Igoof = 1
              else if (epscor.lt.epslon(i)) then
                  Igoof = 2
              else if (Tcore.lt.T(i)) then
                  Igoof = 3
              else
                  Igoof = 0
              end if
          end if
          if (Igoof.ne.-1) then
              istop = i
              go to 50
          end if
c
c  Is it time to change the step size?
c
          if (idrflg.eq.0 .and. M_r(i).lt.0.99d0*Ms) then
              deltar = -Rs/100.0d0
              idrflg = 1
          end if
          if (idrflg.eq.1 .and.deltar.ge.0.5*r(i)) then
              deltar = -Rs/5000.0d0
              idrflg = 2
          end if
          istop = i
   40 continue
c
c  Generate warning messages for the central conditions.
c
      rhocor = M_r(istop)/(4.0d0/3.0d0*pi*r(istop)**3)
      epscor = L_r(istop)/M_r(istop)
      Pcore  = P(istop) + 2.0d0/3.0d0*pi*G*rhocor**2*r(istop)**2
      Tcore  = Pcore*mu*m_H/(rhocor*k_B)
   50 continue
      if (Igoof.ne.0) then
          if (Igoof.eq.-1) then
              write(*,5000) 
              write(*,5100)
          else if (Igoof.eq.1) then
              write(*,6000)
              write(*,5200) rho(istop)
              if (rhocor.gt.1.0d10) write(*,5300)
          else if (Igoof.eq.2) then
              write(*,6000)
              write(*,5400) epslon(istop)
          else if (Igoof.eq.3) then
              write(*,6000)
              write(*,5500) T(istop)
          else if (Igoof.eq.4) then
              write(*,5000) 
              write(*,5600)
          else if (Igoof.eq.5) then
              write(*,5000) 
              write(*,5700)
          else if (Igoof.eq.6) then
              write(*,5000)
              write(*,5800)
          end if
      else
          write(*,7000)
      end if
c
c  Print the central conditions.  If necessary, set limits for the
c  central radius, mass, and luminosity if necessary, to avoid format
c  field overflows.
c
      Rcrat = r(istop)/Rs
      if (Rcrat.lt.-9.999d0) Rcrat = -9.999d0
      Mcrat = M_r(istop)/Ms
      if (Mcrat.lt.-9.999d0) Mcrat = -9.999d0
      Lcrat = L_r(istop)/Ls
      if (Lcrat.lt.-9.999d0) Lcrat = -9.999d0
      write( *,2000) Msolar, Mcrat, Rsolar, Rcrat, Lsolar, Lcrat, Te,
     1 rhocor, X, Tcore, Y, Pcore, Z, epscor, dlPdlT(istop)
      write(20,1000)
      write(20,2000) Msolar, Mcrat, Rsolar, Rcrat, Lsolar, Lcrat, Te,
     1 rhocor, X, Tcore, Y, Pcore, Z, epscor, dlPdlT(istop)
      write(20,2500) Ms
c
c  Print data from the center of the star outward, labeling convective
c  or radiative zones by c or r, respectively.  If abs(dlnP/dlnT)
c  exceeds 99.9, set a print warning flag (*) and set the output limit 
c  to +99.9 or -99.9 as appropriate to avoid format field overflows.
c
      write(20,3000)
      do 60 ic = 1, istop
          i = istop - ic + 1
          Qm = 1.0d0 - M_r(i)/Ms
          if (dlPdlT(i).lt.gamrat) then
              rcf = 'c'
          else
              rcf = 'r'
          end if
          if (abs(dlPdlT(i)).gt.dlPlim) then
              dlPdlT(i) = sign(dlPlim,dlPdlT(i))
              clim = '*'
          else
              clim = ' '
          end if
          write(20,4000) r(i), Qm, L_r(i), T(i), P(i), rho(i), kappa(i),
     1        epslon(i), clim, rcf, dlPdlT(i)
   60 continue
      write(*,9000)
c
c  Format statements
c
  100 format(' ',/,' You must have X + Z <= 1',/,
     1 ' please reenter composition',/)
  200 format(' ',/,
     1 ' The variation in mass has become larger than 0.001*Ms',/,
     2 ' leaving the approximation loop before Nstart was reached',/)
  300 format(' ',/,' The problem occurred in the Runge-Kutta routine',/)
  400 format(' Values from the previous zone are:',/,
     1 10x,'r/Rs   = ',1pe12.5,'  ',
     2 12x,'rho     = ',1pe12.5,' g/cm**3',/,
     3 10x,'M_r/Ms = ',1pe12.5,'  ',
     4 12x,'kappa   = ',1pe12.5,' cm**2/g',/,
     5 10x,'T      = ',1pe12.5,' K',
     6 12x,'epsilon = ',1pe12.5,' ergs/g/s',/,
     7 10x,'P      = ',1pe12.5,' dynes/cm**2',/,
     8 10x,'L_r/Ls = ',1pe12.5)
 1000 format(' ',15x,'A Homogeneous Main-Sequence Model',/)
 2000 format(' ',/,
     1 ' The surface conditions are:',10x,'The central conditions are:',
     2 //,
     3 ' Mtot = ',0pf13.6,' Msun',12x,'Mc/Mtot     = ',1pe12.5,/,
     4 ' Rtot = ',0pf13.6,' Rsun',12x,'Rc/Rtot     = ',1pe12.5,/,
     5 ' Ltot = ',0pf13.6,' Lsun',12x,'Lc/Ltot     = ',1pe12.5,/,
     6 ' Teff = ',0pf13.6,' K   ',12x,'Density     = ',1pe12.5,
     7 ' g/cm**3',/,
     8 ' X    = ',0pf13.6,'     ',12x,'Temperature = ',1pe12.5,' K',/,
     9 ' Y    = ',0pf13.6,'     ',12x,'Pressure    = ',1pe12.5,
     1 ' dynes/cm**2',/,
     2 ' Z    = ',0pf13.6,'     ',12x,'epsilon     = ',1pe12.5,
     3 ' ergs/s/g',/,
     4 '        ',   13x ,'     ',12x,'dlnP/dlnT   = ',1pe12.5,//)
 2500 format(' ',/,' Notes: ',/,
     1 ' (1) Mass is listed as Qm = 1.0 - M_r/Mtot, where Mtot = ',
     2 1pe13.6,' g',/,
     3 ' (2) Convective zones are indicated by c, radiative zones by r',
     4 /,
     5 ' (3) dlnP/dlnT may be limited to +99.9 or -99.9; if so it is',
     6 ' labeled by *',//)
 3000 format(' ',5x,'r',7x,'Qm',7x,'L_r',7x,'T',8x,'P',7x,
     1 'rho',6x,'kap',6x,'eps',3x,'dlPdlT')
 4000 format(' ',1p8e9.2,2a1,0pf5.1)
 5000 format(' ',/,15x,'Sorry to be the bearer of bad news, but...',/,
     1 15x,'       Your model has some problems',/)
 5100 format(' ',   8x,'The number of allowed shells has been exceeded',
     1 /)
 5200 format(' ', 14x,'The core density seems a bit off,'/,
     1 5x,' density should increase smoothly toward the center.',/,
     2 5x,' The density of the last zone calculated was rho = ',
     3 1pe10.3,' gm/cm**3',/)
 5300 format(' ',   1x,'It looks like you will need a degenerate',
     1 ' neutron gas and general relativity',/,
     2 ' to solve this core.  Who do you think I am, Einstein?',/)
 5400 format(' ',  14x,'The core epsilon seems a bit off,',/,
     1 9x,' epsilon should vary smoothly near the center.',/,
     2 9x,' The value calculated for the last zone was eps =',
     3 1pe10.3,' ergs/g/s',/)
 5500 format(' ',8x,' Your extrapolated central temperature is too low',
     1 /,8x,' a little more fine tuning ought to do it.',/,
     2 8x,' The value calculated for the last zone was T = ',
     3 1pe10.3,' K',/)
 5600 format(' ',  10x,'You created a star with a hole in the center!',
     1 /)
 5700 format(' ',  10x,'This star has a negative central luminosity!',/)
 5800 format(' ',  5x,'You hit the center before the mass and/or ',
     1 'luminosity were depleted!',/)
 6000 format(' ',///,15x,'It looks like you are getting close,',/,
     1 12x,'however, there are still a few minor errors',/)
 7000 format(' ',///,15x,'CONGRATULATIONS, I THINK YOU FOUND IT!',/,
     1 9x,'However, be sure to look at your model carefully.',//)
 9000 format(' ',10x,'***** The integration has been completed *****',
     1 /,10x,'   The model has been stored in starmodl.dat',
     2 /)
      stop
      end
c
c  Subroutine STARTMDL computes values of M_r, L_r, P, and T, near the
c  surface of the star using expansions of the stellar structure
c  equations (M_r and L_r are assumed to be constant).
c
      Subroutine STARTMDL(deltar, X, Z, mu, Rs, r_i, M_ri, L_ri, r,
     1    P_ip1, M_rip1, L_rip1, T_ip1, tog_bf, irc)
      real*8 deltar, X, Z, mu, Rs, r_i, M_ri, L_ri, r, P_ip1, M_rip1,
     1    L_rip1, T_ip1, tog_bf
      real*8 A_bf, A_ff, Afac
      real*8 sigma, c, a, G , k_B, m_H, pi, gamma, gamrat, kPad, g_ff
      common /cnstnt/ sigma, c, a, G, k_B, m_H, pi, gamma, gamrat, kPad,
     1    g_ff
c
      r = r_i + deltar
      M_rip1 = M_ri
      L_rip1 = L_ri
c
c  This is the radiative approximation (neglect radiation pressure
c  and electron scattering opacity); see Eqs. (H.1), (H.2), (9.19),
c  and (9.20).
c
      if (irc.eq.0) then
          T_ip1 = G*M_rip1*mu*m_H/(4.25d0*k_B)*(1.0d0/r - 1.0d0/Rs)
          A_bf = 4.34d25*Z*(1.0d0 + X)/tog_bf
          A_ff = 3.68d22*g_ff*(1.0d0 - Z)*(1.0d0 + X)
          Afac = A_bf + A_ff
          P_ip1 = sqrt((1.0d0/4.25d0)*(16.0d0/3.0d0*pi*a*c)*
     1        (G*M_rip1/L_rip1)*(k_B/(afac*mu*m_H)))*T_ip1**4.25d0
c
c  This is the convective approximation; see Eqs. (H.3) and (10.75).
c
      else
          T_ip1 = G*M_rip1*mu*m_H/k_B*(1.0d0/r - 1.0d0/Rs)/gamrat
          P_ip1 = kPad*T_ip1**gamrat
      end if
      return
      end
c
c  Subroutine EOS calculates the values of density, opacity, the 
c  guillotine-to-gaunt factor ratio, and the energy generation rate at
c  the radius r.
c
      Subroutine EOS(X, Z, XCNO, mu, P, T, rho, kappa, epslon, tog_bf,
     1    izone, ierr)
      real*8 X, Z, mu, P, T, rho, kappa, epslon, tog_bf, Prad, Pgas
      real*8 k_bf, k_ff, k_e
      real*8 T6, fx, fpp, epspp, epsCNO, XCNO, oneo3, twoo3, psipp, Cpp
      real*8 sigma, c, a, G, k_B, m_H, pi, gamma, gamrat, kPad, g_ff
      common /cnstnt/ sigma, c, a, G, k_B, m_H, pi, gamma, gamrat, kPad,
     1    g_ff
      data oneo3, twoo3 / 0.333333333d0, 0.666666667d0 /
c
c  Solve for density from the ideal gas law (remove radiation pressure);
c  see Eq. (10.26).
c
      if (T.le.0.0d0 .or. P.le.0.0d0) then
        ierr = 1
        write(*,100) izone, T, P
        return
      end if
      Prad = a*T**4/3.0d0
      Pgas = P - Prad
      rho=(mu*m_H/k_B)*(Pgas/T)
      if (rho.lt.0.0d0) then
        ierr = 1
        write(*,200) izone, T, P, Prad, Pgas, rho
        return
      end if
c
c  Calculate opacity, including the guillotine-to-gaunt factor ratio;
c  see Novotny (1973), p. 469. k_bf, k_ff, and k_e are the bound-free,
c  free-free, and electron scattering opacities, given by Eqs. (9.19),
c  (9.20), and (9.21), respectively.
c
      tog_bf = 2.82d0*(rho*(1.0d0 + X))**0.2d0
      k_bf = 4.34d25/tog_bf*Z*(1.0d0 + X)*rho/T**3.5d0
      k_ff = 3.68d22*g_ff*(1.0d0 - Z)*(1.0d0 + X)*rho/T**3.5d0
      k_e = 0.2d0*(1.0d0 + X)
      kappa = k_bf + k_ff + k_e
c
c  Compute energy generation by the pp chain and the CNO cycle.  These
c  are calculated using Eqs. (10.49) and (10.53), which come from 
c  Fowler, Caughlan, and Zimmerman (1975). The screening factor for the
c  pp chain is calculated as fpp; see Clayton (1968), p. 359ff.
c
      T6 = T*1.0d-06
      fx = 0.133d0*X*sqrt((3.0d0 + X)*rho)/T6**1.5d0
      fpp = 1.0d0 + fx*X
      psipp = 1.0d0 + 1.412d8*(1.0d0/X - 1.0d0)*exp(-49.98*T6**(-oneo3))
      Cpp = 1.0d0 + 0.0123d0*T6**oneo3 + 0.0109d0*T6**twoo3
     1    + 0.000938d0*T6
      epspp = 2.38d6*rho*X*X*fpp*psipp*Cpp*T6**(-twoo3)
     1    *exp(-33.80d0*T6**(-oneo3))
      CCNO = 1.0d0 + 0.0027d0*T6**oneo3 - 0.00778d0*T6**twoo3
     1    - 0.000149d0*T6
      epsCNO = 8.67d27*rho*X*XCNO*CCNO*T6**(-twoo3)
     1    *exp(-152.28d0*T6**(-oneo3))
      epslon = epspp + epsCNO
c
c  Formats
c
  100 format(' ',/,' Something is a little wrong here.',
     1 /,' You are asking me to deal with either a negative temperature'
     2 /,' or a negative pressure.  I am sorry but that is not in my'
     3  ,' contract!',/,' You will have to try again with different',
     4   ' initial conditions.',/,
     5   ' In case it helps, I detected the problem in zone ',i3,
     6   ' with the following',/,' conditions:',/,
     7 10x,'T = ',1pe10.3,' K',/,
     8 10x,'P = ',1pe10.3,' dynes/cm**2')
  200 format(' ',/,' I am sorry, but a negative density was detected.',
     1 /,' my equation-of-state routine is a bit baffled by this new',
     2 /,' physical system you have created.  The radiation pressure',
     3 /,' is probably too great, implying that the star is unstable.'
     4 /,' Please try something a little less radical next time.',/,
     5   ' In case it helps, I detected the problem in zone ',i3,
     6   ' with the following',/,' conditions:',/,
     7 10x,'T       = ',1pe10.3,' K',/,
     8 10x,'P_total = ',1pe10.3,' dynes/cm**2',/,
     9 10x,'P_rad   = ',1pe10.3,' dynes/cm**2',/,
     1 10x,'P_gas   = ',1pe10.3,' dynes/cm**2',/,
     2 10x,'rho     = ',1pe10.3,' g/cm**3')
      return
      end        
c 
c  The following four function subprograms calculate the gradients of
c  pressure, mass, luminosity, and temperature at r.
c    
      real*8 Function dPdr(r, M_r, rho)
      real*8 r, M_r, rho
      real*8 sigma, c, a, G, k_B, m_H, pi, gamma, gamrat, kPad, g_ff
      common /cnstnt/ sigma, c, a, G, k_B, m_H, pi, gamma, gamrat, kPad,
     1    g_ff
c  Eq. (10.7)
      dPdr = -G*rho*M_r/r**2
      return
      end
c     
      real*8 Function dMdr(r, rho)
      real*8 r, rho
      real*8 sigma, c, a, G, k_B, m_H, pi, gamma, gamrat, kPad, g_ff
      common /cnstnt/ sigma, c, a, G, k_B, m_H, pi, gamma, gamrat, kPad,
     1    g_ff
c  Eq. (10.8)
      dMdr = 4.0d0*pi*rho*r**2
      return
      end
c 
      real*8 Function dLdr(r, rho, epslon)
      real*8 r, rho, epslon
      real*8 sigma, c, a, G, k_B, m_H, pi, gamma, gamrat, kPad, g_ff
      common /cnstnt/ sigma, c, a, G, k_B, m_H, pi, gamma, gamrat, kPad,
     1    g_ff
c  Eq. (10.45)
      dLdr = 4.0d0*pi*rho*epslon*r**2
      return
      end
c    
      real*8 Function dTdr(r, M_r, L_r, T, rho, kappa, mu, irc)
      real*8 r, M_r, L_r, T, rho, kappa, mu
      real*8 sigma, c, a, G, k_B, m_H, pi, gamma, gamrat, kPad, g_ff
      common /cnstnt/ sigma, c, a, G, k_B, m_H, pi, gamma, gamrat, kPad,
     1    g_ff
c  This is the radiative temperature gradient (Eq. 10.61).
      if (irc.eq.0) then
        dTdr = -(3.0d0/(16.0d0*pi*a*c))*kappa*rho/T**3*L_r/r**2
c  This is the adiabatic convective temperature gradient (Eq. 10.81).
      else
        dTdr = -1.0d0/gamrat*G*M_r/r**2*mu*m_H/k_B
      end if
      return
      end
c
c  This is a fourth-order Runge-Kutta integration routine.
c
      Subroutine RUNGE(f_im1, dfdr, f_i, r_im1, deltar, irc, X, Z, XCNO,
     1    mu, izone, ierr)
      real*8 f_im1(4), dfdr(4), f_i(4), df1(4), df2(4), df3(4),
     1    f_temp(4)
      real*8 r_im1, r_i, deltar, dr12, dr16, r12, X, Z, XCNO, mu
c      
      dr12 = deltar/2.0d0
      dr16 = deltar/6.0d0
      r12  = r_im1 + dr12
      r_i  = r_im1 + deltar
c
c  Calculate intermediate derivatives from the fundamental stellar
c  structure equations found in Subroutine FUNDEQ.
c
      do 10 i = 1, 4                     
          f_temp(i) = f_im1(i) + dr12*dfdr(i)
   10 continue
      call FUNDEQ(r12, f_temp, df1, irc, X, Z, XCNO, mu, izone, ierr)
      if (ierr.ne.0) return
c
      do 20 i = 1, 4
          f_temp(i) = f_im1(i) + dr12*df1(i)
   20 continue
      call FUNDEQ(r12, f_temp, df2, irc, X, Z, XCNO, mu, izone, ierr)
      if (ierr.ne.0) return
c      
      do 30 i = 1, 4
          f_temp(i) = f_im1(i) + deltar*df2(i)
   30 continue
      call FUNDEQ(r_i, f_temp, df3, irc, X, Z, XCNO, mu, izone, ierr)
      if (ierr.ne.0) return
c
c  Calculate the variables at the next shell (i + 1).
c
      do 40 i = 1, 4
          f_i(i) = f_im1(i) + dr16*(dfdr(i) + 2.0d0*df1(i)
     1        + 2.0d0*df2(i) + df3(i))
   40 continue
      return
      end
c
c  This subroutine returns the required derivatives for RUNGE, the
c  Runge-Kutta integration routine.
c
      Subroutine FUNDEQ(r, f, dfdr, irc, X, Z, XCNO, mu, izone, ierr)
      real*8 f(4), dfdr(4), X, Z, XCNO, r, M_r, P, T, L_r, rho, kappa,
     1    epslon, tog_bf, mu
      real*8 dPdr, dMdr, dLdr, dTdr
c
      P   = f(1)
      M_r = f(2)
      L_r = f(3)
      T   = f(4)
      call EOS(X, Z, XCNO, mu, P, T, rho, kappa, epslon, tog_bf, izone,
     1    ierr)
      dfdr(1) = dPdr(r, M_r, rho)
      dfdr(2) = dMdr(r, rho)
      dfdr(3) = dLdr(r, rho, epslon)
      dfdr(4) = dTdr(r, M_r, L_r, T, rho, kappa, mu, irc)
      return
      end
