program boltzmann123
! ktk Jan 2008

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

       double precision k, eV, Z, part
       double precision g1,g2,g3,E1,E2,E3,T, Tincr, Tmin, Tmax
       double precision n21(npts),n31(npts),n32(npts)
       double precision n1rat(npts), n2rat(npts), n3rat(npts)

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

       print *
       write(6,*)'computes H-like rel.level pops (1,2,3) for range in T'
       print *

       write(6,*)'enter number of protons in H-like specie; Z=1 for H'
       read(5,*) Z

* some constants CODATA (2006)

       k = 1.3806504d-16     ! Boltzmann's constant
       eV = 1.602176487d-12  ! ergs/eV
       E1 = -13.5984340d+00 * Z * Z  
! ignore contributions from reduced mass which differ from H
       E2 = E1 / (2.d+00*2.d+00)
       E3 = E1 / (3.d+00*3.d+00)
       g1 = 2.d+00
       g2 = 8.d+00
       g3 = 18.d+00

       print *
       write(6,*)'Also computes N_n/Ntot for n = 1,2,3 w/ less accuracy'

      write(6,*)'because simplified H-like partition function used here'

       write(6,*)'blows up at large n'
       print *

       write(6,*)'enter max n for H-like partition function sum.'
       write(6,*)'I suggest limiting nmax < 20; Try nmax = 10.'

       read(5,*) ncut

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

       write(16,96) Z
96     format('# relative level populations for H-like Z = ', f3.0)
       write(16,97)
97     format('# assumed form of H-like part.func uncertain at high T')
       write(16,98) ncut
98     format('# max. n level considered in partition function: ', i3)
       write(16,99)
99     format('# T(K)       N21        N31        N32        N1/HI    
     -N2/HI      N3/HI      PFn')

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

       do i = 1, npts

       T = T + Tincr

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

! compute absolute level population; need partition function...
! compute the partition function, cutoff at ncut = n < 100
       part = 0.
       do n = 1, ncut
        part = part + 2.d+00*n*n * 
     -  dexp( E1*eV*( 1.d+00 - 1.d+00/(n*n) )/(k*T) )
       enddo

       n1rat(i) = (g1/part) * 
     - dexp( E1*eV*( 1.d+00 - 1.d+00/(1.d+00*1.d+00) )/(k*T) )
       n2rat(i) = (g2/part) * 
     - dexp( E1*eV*( 1.d+00 - 1.d+00/(2.d+00*2.d+00) )/(k*T) )
       n3rat(i) = (g3/part) * 
     - dexp( E1*eV*( 1.d+00 - 1.d+00/(3.d+00*3.d+00) )/(k*T) )

       write(16,100)T,n21(i),n31(i),n32(i),n1rat(i),n2rat(i),n3rat(i),
     -              part
100    format(1p, 7(e9.3, 2x), e10.4 )

       enddo
       print *
       close(16)

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

       stop
       end