program boltsaha
! ktk Jan.2007

       implicit none
       
       integer npts, i
       parameter (npts = 1001)

       double precision k, eV, Z
       double precision g1,g2,g3,g4,E1,E2,E3,E4
       double precision T, Tincr, Tmin, Tmax
       double precision n21(npts),n31(npts),n41(npts)
       double precision n32(npts),n42(npts),n43(npts)
       double precision n2saha(npts), n3saha(npts)

       double precision c, h, k, me, eV, A1, A2, pi, sigma
       double precision Pe, T, Tincr, Tmin, Tmax
       double precision Zp, Zm, Zo, chi1, chim       ! for H
       double precision Z4pp, Z4p, Z4o, chi41, chi42  ! for He+
       double precision saha1(npts), saha2(npts)


       open(16,file='bolt2saha.out',status='unknown')

       print *
       write(6,*)'compute Nx/Ntot for H or He+, assumes 4-level atom '
       print *

       write(6,*)'enter Z, equal to 1 for H, or 2 for He+'
       read(5,*) Z
       print *

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

       if (Z .eq. 1) then ! H
         write(16,94)
94       format('# approximate N2/Ntot for H ')
         write(16,95) Pe
95       format('# the assumed cgs electron pressure is ', 1p, e9.3)
         write(16,96)
96       format('# T(K)      N2/Ntot')

       else

         write(16,97)
97       format('# approximate N3/Ntot for He+ ')
         write(16,98) Pe
98       format('# the assumed cgs electron pressure is ', 1p, e9.3)
         write(16,99)
99       format('# T(K)      N3/Ntot')

       endif

* some constants CODATA (2006)
* for Boltzmann
       k = 1.3806504d-16     ! Boltzmann's constant
       eV = 1.602176487d-12  ! ergs/eV
! HI, HeII both
       E1 = -13.5984340d+00 * Z * Z
       E2 = E1 / (2.d+00*2.d+00)
       E3 = E1 / (3.d+00*3.d+00)
       E4 = E1 / (4.d+00*4.d+00)
       g1 = 2.d+00
       g2 = 8.d+00
       g3 = 18.d+00
       g4 = 32.d+00

* for Saha
       c = 2.99792458d+10
       h = 6.62606896d-27
       sigma = 5.670400d-05
       me = 9.10938215d-28
       pi = 3.1415926535897932d+00
!       eV = 1.602176487d-12 ! already defined
* the ions
! HI
       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
! HeII
       chi41 = 24.5873876d+00  ! I.P. in eV for He ground state
       chi42 = 54.417763d+00  ! I.P. in eV for He+ ground state
       Z4pp = 1.d+00         ! partition function for He++
       Z4p = 2.d+00          ! partition function for He+
       Z4o = 1.d+00          ! partition function for He in ground state

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

! the temperature array increment
       Tincr = (Tmax - Tmin)/(npts - 1)
       T = Tmin - Tincr

! the temperature loop
       do i = 1, npts

       T = T + Tincr

! the relative Boltzmann populations in the 4-level atom(s)

       n21(i) = (g2/g1) * dexp( -(E2 - E1)*eV/(k*T) )
       n31(i) = (g3/g1) * dexp( -(E3 - E1)*eV/(k*T) )
       n41(i) = (g4/g1) * dexp( -(E4 - E1)*eV/(k*T) )
       n32(i) = (g3/g2) * dexp( -(E3 - E2)*eV/(k*T) )
       n42(i) = (g4/g2) * dexp( -(E4 - E2)*eV/(k*T) )
       n43(i) = (g4/g3) * dexp( -(E4 - E3)*eV/(k*T) )

! now the Saha ion ratios

       if (Z .eq. 1.) then  ! hydrogen

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

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

! the requested ratio
       n2saha(i) = ( n21(i) / ( 1.d0 + n21(i) + n31(i) + n41(i) ) ) / 
     -             ( ( 1.d0/saha1(i) ) + 1.d0 + saha2(i) )

       write(16,100) T, n2saha(i)
100    format(1p, 2(e10.4, 2x))

       else  ! He+

       A1 = 2.d+00* (Z4p/Z4o) * ( 2.d+00*pi*me/(h*h) )**1.5
       A2 = 2.d+00* (Z4pp/Z4p) * ( 2.d+00*pi*me/(h*h) )**1.5

! ionized He relative to neutral He
           saha1(i) =  (A1/Pe) * (k*T)**2.5 * dexp(-chi41*eV/(k*T))
! He++ relative to He+
           saha2(i) =  (A2/Pe) * (k*T)**2.5 * dexp(-chi42*eV/(k*T))

! the requested ratio
       n3saha(i) = ( n32(i) / ( (1.d0/n21(i)) + 1.d0 + n32(i) + n42(i)  
     -           ) ) / ( ( 1.d0/saha1(i) ) + 1.d0 + saha2(i) )

       write(16,101) T, n3saha(i)
101    format(1p, 2(e10.4, 2x))

       endif

       enddo
       print *
       close(16)

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

       stop
       end