program sahaH
! ktk Jan.2007

       implicit none

       integer npts, icode, i
       parameter (npts = 1001)

       double precision c, h, k, me, eV, A1, Am, pi, sigma
       double precision ne, Pe, T, Tincr, Tmin, Tmax
       double precision Zp, Zm, Zo, chi1, chim
       double precision saha1(npts), saham(npts)
       double precision ratio1(npts), ratio2(npts), ratiom(npts)

! some constants CODATA (2006) 
       c = 2.99792458d+10
       h = 6.62606896d-27
       k = 1.3806504d-16
       sigma = 5.67040d-05
       me = 9.10938215d-28
       eV = 1.602176487d-12
       pi = 3.1415926535897932d+00

       chi1 = 13.5984340d+00  ! I.P. in eV for H ground state
       chim = 0.754d+00     ! I.P. in eV for H-minus
       Zp = 1.d+00          ! partition function for H+
       Zo = 2.d+00          ! partition function for Ho in ground state
       Zm = 1.d+00          ! partition function for H-minus


       open(16,file='sahaH.out',status='unknown')
       write(16,68)
68     format('# the Saha ionization fractions for H')

       print *
       write(6,*)'computes adjacent stages of ionization of H at T'
       write(6,*)'assumes neutral H in the ground state (primarily)'

       print *
       write(6,*)'enter min, max T(K)'
       read(5,*) Tmin, Tmax

       Tincr = (Tmax - Tmin)/(npts - 1)
       T = Tmin - Tincr

       A1 = 2.d+00* (Zp/Zo) * ( 2.d+00*pi*me/(h*h) )**1.5
       Am = 2.d+00* (Zo/Zm) * ( 2.d+00*pi*me/(h*h) )**1.5

       print *
       write(6,*)'enter electron density(1)? or pressure(2)? '
       read(5,*) icode

       if (icode .eq. 1) then
         write(6,*)'enter the electron number density cgs '
         read(5,*) ne
         Pe = ne*k*T

       write(16,79) ne
79     format('# the assumed cgs electron density is ', 1p, e9.3)
       write(16,80)
80     format('#  T      HII/HI   HI/Hminus   HII/H   HI/H   Hminus/H')

         do i = 1, npts

           T = T + Tincr

! ionized relative to neutral H
           saha1(i) =  (A1/ne) * (k*T)**1.5 * dexp(-chi1*eV/(k*T))
! HI relative to Hminus
           saham(i) =  (Am/ne) * (k*T)**1.5 * dexp(-chim*eV/(k*T))

           ratio2(i) = saha1(i)/(1.d+00 + saha1(i) + 1.d+00/saham(i))
           ratio1(i) = 1.d+00/(1.d+00 + saha1(i) + 1.d+00/saham(i))
           ratiom(i) = ratio1(i) / saham(i)

         write(16,100) T,saha1(i),saham(i),ratio2(i),ratio1(i),ratiom(i)
100        format(1p, 6(e10.4, 2x))

         enddo

       else
         write(6,*)'enter the electron pressure cgs '
         read(5,*) Pe

       write(16,89) Pe
89     format('# the assumed cgs electron pressure is ', 1p, e9.3)
       write(16,90)
90     format('#   T       HII/HI      HI/Hminus     HII/H       HI/H
     -   Hminus/H')

         do i = 1, npts

           T = T + Tincr

! ionized relative to neutral H
           saha1(i) =  (A1/Pe) * (k*T)**2.5 * dexp(-chi1*eV/(k*T))
! HI relative to Hminus
           saham(i) =  (Am/Pe) * (k*T)**2.5 * dexp(-chim*eV/(k*T))

           ratio2(i) = saha1(i)/(1.d+00 + saha1(i) + 1.d+00/saham(i))
           ratio1(i) = 1.d+00/(1.d+00 + saha1(i) + 1.d+00/saham(i))
           ratiom(i) = ratio1(i) / saham(i)

         write(16,101) T,saha1(i),saham(i),ratio2(i),ratio1(i),ratiom(i)
101        format(1p, 6(e10.4, 2x))

         enddo

       endif
       print *
       close(16)

       write(6,*)' output sent to sahaH.out'
       print *

       stop
       end