program sahaCa
! ktk Jan.2007

       implicit none

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

       double precision c, h, k, me, eV, A1, A2, A3, pi, sigma
       double precision ne, Pe, T, Tincr, Tmin, Tmax
       double precision Zp3, Zpp, Zp, Zo, chi1, chi2, chi3
       double precision saha3(npts), 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 = 6.113158d+00        ! I.P. in eV for Ca ground state
       chi2 = 11.87172d+00       ! I.P. in eV for Ca+ ground state
       chi3 = 50.9131d+00         ! I.P. in eV for Ca++ ground state
       Zp3 = 1.d+00              ! Aprrox. partition function for Ca+3 ? 
       Zpp = 1.d+00              ! Aprrox. partition function for Ca++ ?
       Zp = 2.30d+00             ! Aprrox. partition function for Ca+
       Zo = 1.32d+00             ! Aprrox. partition function for Ca0

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

       print *
       write(6,*)'computes adjacent stages of ionization of Ca at T'
       write(6,*)'assumes Ca 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
       A3 = 2.d+00 * (Zp3/Zpp) * ( 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) Pe
79     format('# the assumed cgs electron density is ', 1p, e9.3)
       write(16,80)
80     format('#  T      CaIII/CaII   CaII/CaI   CaIII/Ca   CaII/Ca    
     cCaI/Ca')

         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))
           saha3(i) =  (A3/ne) * (k*T)**1.5 * dexp(-chi3*eV/(k*T))

           ratio1(i) = (1.d0/saha1(i)) / ( (1.d0/saha1(i)) + 1.d0
     -                  + saha2(i) + (saha3(i) * saha2(i)) )
           ratio2(i) = ratio1(i) * saha1(i)
           ratio3(i) = ratio2(i) * saha2(i)


         write(16,100) T, saha2(i), saha1(i), ratio3(i), ratio2(i), 
     -                 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      CaIII/CaII   CaII/CaI   CaIII/Ca   CaII/Ca    
     cCaI/Ca')

         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.d0/saha1(i)) / ( (1.d0/saha1(i)) + 1.d0
     -                  + saha2(i) + (saha3(i) * saha2(i)) )
           ratio2(i) = ratio1(i) * saha1(i)
           ratio3(i) = ratio2(i) * saha2(i)


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

         enddo

       endif
       print *
       close(16)

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

       stop
       end