c    program for saha solution of H, He mix at different rho, T.
c    prints out to the screen, not a file- when running, redirect
c    output to desired filename.
c
c    To compile and run in unix: 
c    rename file so it has a .f extension
c    f77 -o saha -r8 f77saha.f
c    saha > tablefile
c
      dimension rn(5),gn(5),f(2)
c
c    rn: computed number fractions for HI,HII,HeI,HeII,HeIII
c    gn: statistical weights for the above
c    f:  number fraction abundance of H, He
c
      data  gn/2.,1.,1.,2.,1./
      data  f/.8,.2/
c
      rk=8.6174e-5
c       =boltzmann const in ev.
      pk=2.0704e-16
c       =combination of constants in saha equation.
c        note: usually a good idea to set variables rather than write
c              in numbers each time.  program runs slightly faster.
c
c    do calculations at three chosen electron densities.
      do 30 n=1,3
       enl=14+n*2
       en=10**enl
c        =electron number density.
c     next line prints column titles for the output table.
c
       print 5, enl
    5  format(' enl=',f5.1/'tl, rn(5), ef, anl')
c
c     set initial HII fraction equal to total H fraction.
       rn(2)=f(1)
       an=en
c
c     for each density, do calculations for 50 temperatures.
       do 20 m=1,50
        rm=m-1
        tl=4.+rm/30
        t=10**(tl)
        rkt=rk*t
        rn(2)=rn(2)*an
        rn(1)=en*rn(2)*pk*gn(1)*t**(-1.5)*exp(13.595/(rkt))/gn(2)
             =HI density given HII density.
c
c     (diagnostic prints):
c      print 6, t,en,rn(2),rkt,gn(1),gn(2),rn(1)
c   6  format('t,e,r2,kt,g1,g2,r1'/7e11.2)
c
        hn=rn(1)+rn(2)
c         =total H given original HII.
        he=hn*f(2)/f(1)
c         =total He given above total H.
        rn(4)=he
c            =assume all He is HeII.
        rn(3)=en*rn(4)*pk*gn(3)*t**(-1.5)*exp(24.587/(rkt))/gn(4)
c            =HeI given HeII.
        rn(5)=rn(4)*gn(5)*t**(1.5)*exp(-54.416/(rkt))/(en*pk*gn(4))
c            =HeIII given HeII.
c
c     (diagnostic prints):
c      print 7, hn,he,rn(3),rn(4),rn(5)
c   7  format('hn,he,r3,r4,r5'/5e11.2)
c
        sum=rn(3)+rn(4)+rn(5)
c          =new revised total He
        fnorm=he/sum
c            =scale factor to make He correct relative to H.
c
c     (diagnostic prints):
c      print 8, sum,he,fnorm
c   8  format('sum,he,fnorm'/5e11.2)
c
        rn(3)=fnorm*rn(3)
        rn(4)=fnorm*rn(4)
        rn(5)=fnorm*rn(5)
c             all He's now correct relative to H.
        csum=rn(2)+rn(4)+2.*rn(5)
c           =charge sum. should be = to electron density.
        fnorm=en/csum
c            =scale factor to correct all abundances to give correct electron density.
c
        an=0.
c     do correction scaling, and sum nuclear density.
        do 10 i=1,5
         rn(i)=fnorm*rn(i)
         an=an+rn(i)
   10    continue
c
c     divide by an to produce fractional abundances.
        do 11 i=1,5
         rn(i)=rn(i)/an
   11    continue
        ef=en/an
        anl=alog10(an)
c
c     print table entries.
        print 15, tl,(rn(i),i=1,5),ef,anl
   15   format(f7.4,5f8.5, 2f8.4)
   20   continue
   30  continue
      stop
      end