program sahaHe
! ktk Jan.2007

       implicit none

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

       double precision c, h, k, me, eV, A1, A2, pi, sigma
       double precision ne, Pe, T, Tincr, Tmin, Tmax
       double precision Zpp, Zp, Zo, chi1, chi2
       double precision saha2(npts), saha1(npts)
       double precision ratio1(npts), ratio2(npts), ratio3(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 = 24.5873876d+00     ! I.P. in eV for He ground state
       chi2 = 54.417763d+00     ! I.P. in eV for He+ ground state
       Zpp = 1.d+00            ! partition function for He++
       Zp = 2.d+00             ! partition function for He+
       Zo = 1.d+00             ! partition function for He0


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


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

       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
       A2 = 2.d+00* (Zpp/Zp) * ( 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      HeIII/HeII   HeII/HeI   HeIII/He   HeII/He
     cHeI/He')

         do i = 1, npts

           T = T + Tincr

           saha1(i) =  (A1/ne) * (k*T)**1.5 * dexp(-chi1*eV/(k*T))
           saha2(i) =  (A2/ne) * (k*T)**1.5 * dexp(-chi2*eV/(k*T))

           ratio1(i) = (1./saha1(i))/(1./saha1(i) + 1.+ saha2(i))
           ratio2(i) = 1./(1./saha1(i) + 1.+ saha2(i))
           ratio3(i) = 1.d+00 - ( ratio1(i) + ratio2(i) )


           write(16,100) T, saha2(i), saha1(i), ratio3(i), ratio2(i), 
     c                   ratio1(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        HeIII/HeII   HeII/HeI   HeIII/He    HeII/He
     c   HeI/He')

         do i = 1, npts

           T = T + Tincr

           saha1(i) =  (A1/Pe) * (k*T)**2.5 * dexp(-chi1*eV/(k*T))
           saha2(i) =  (A2/Pe) * (k*T)**2.5 * dexp(-chi2*eV/(k*T))

           ratio1(i) = (1.d+00/saha1(i))/(1./saha1(i) + 1.+ saha2(i))
           ratio2(i) = 1.d+00/(1.d+00/saha1(i) + 1.d+00+ saha2(i))
           ratio3(i) = 1.d+00 - ( ratio1(i) + ratio2(i) )

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

         enddo

       endif
       print *
       close(16)

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

       stop
       end