c###################################################################
c
      subroutine heionen(n, l, mult, sumen, eion, kermsg)
c
c     this subroutine passes the ionization energy for excited states
c     of  helium taken from tables a.4 and a.5 given by janev et al.
c     (see doc=h-he-plasma.)
c
c     the input subroutine parameters are
c
c    n = principal quantum number of the excited electronic state
c
c    l = orbital angular momentum quantum number of the exited state
c
c    mult = the spin multiplicity (2s+1) of the state
c
c    sumen =  indicates choice of ionization energy to be returned.
c             if sumen=0  the ionization energy for the specific
c              state (quantified by n, l amd mult) is returned
c             if sumen=1  the ionization energy taken as an average
c             over angular momentum and toatal spin is returned
c
c     the output subroutine parameters are
c
c    eion = ionization energy
c
c    kermsg = blank if no errors
c
c     written by j. j. smith , iaea atomic and molecular data unit
c
c------------------------------------------------------------------------
c
      character*(*) kermsg
      integer n, l, mult, sumen, lp, multp
      dimension enl(4,4,2), en(7)
      data (((enl(i,j,k),i=1,4),j=1,4),k=1,2) /32*0.0/
c
c     data enl(1,1,1)
c    *    / 24.588 /,
c
c----   n=2 excited states
c
c    1  enl(2,1,1),enl(2,1,2),enl(2,2,1),enl(2,2,2)
c    *    / 3.973, 4.769, 3.371, 3.625 /,
c
c----   n=3 excited states
c
c    2   enl(3,1,1),enl(3,1,2),enl(3,2,1),enl(3,2,2),enl(3,3,1),
c    *   enl(3,3,2)
c    *    / 1.699, 1.871, 1.502, 1.582, 1.515, 1.515 /,
c
c----   n=4 excited states
c
c    3   enl(4,1,1),enl(4,1,2),enl(4,2,1),enl(4,2,2),enl(4,3,1),
c    *   enl(4,3,2),enl(4,4,1),enl(4,4,2)
c    *    / 0.9155, 0.9951, 0.8470, 0.8812, 0.8529, 0.8530, 0.8520,
c    *      0.8520 /
c
c-----   ionization energies summed over angular momentum and
c-----   total spin
c
      data (en(k),k=1,2) /2*0.0/
      data en(3),en(4),en(5),en(6),en(7)
     *    / 1.609 , 0.8811, 0.52, 0.37, 0.29/
c
      enl(1,1,1)=24.588
c----   n=2 excited states
      enl(2,1,1)=3.973
      enl(2,1,2)=4.769
      enl(2,2,1)=3.371
      enl(2,2,2)=3.625
c----   n=3 excited states
      enl(3,1,1)=1.699
      enl(3,1,2)=1.871
      enl(3,2,1)=1.502
      enl(3,2,2)=1.582
      enl(3,3,1)=1.515
      enl(3,3,2)=1.515
c----   n=4 excited states
      enl(4,1,1)=0.9155
      enl(4,1,2)=0.9951
      enl(4,2,1)=0.8470
      enl(4,2,2)=0.8812
      enl(4,3,1)=0.8529
      enl(4,3,2)=0.8530
      enl(4,4,1)=0.8520
      enl(4,4,2)=0.8520
      kermsg =' '
      if (sumen .eq. 1) then
        if (n .ge. 8) then
          eion =13.58/(n*n)
          return
        else
           eion =en(n)
           if (eion .eq. 0.0) kermsg =
     *        'ionization energy for n value not in table in heionen'
        endif
      else
           lp = l+1
           if (mult .eq. 1) then
             multp = 1
           else if (mult .eq. 3) then
               multp = 2
             else
                 kermsg = 'invalid spin multipclity as input in heionen'
             endif
           eion = enl(n,lp,multp)
           if (eion .eq. 0.0) kermsg =
     *        'ionization energy not in table in heionen'
      endif
      return
      end
c
c###################################################################
c
      subroutine heexcen(n, l, mult, sumen, eexc, kermsg)
c
c     this subroutine passes the excitation energy for excited states
c     of  helium taken from tables a.4 and a.5 given by janev et al.
c     (see doc=h-he-plasma.)
c
c     the input subroutine parameters are
c
c    n = principal quantum number of the excited electronic state
c
c    l = orbital angular momentum quantum number of the exited state
c
c    mult = the spin multiplicity (2s+1) of the state
c
c    sumen =  indicates choice of excitation enegy to be returned.
c             if sumen=0  the excitation energy for the specific
c              state (quantified by n, l amd mult) is returned
c             if sumen=1  the excitation energy taken as an average
c             over angular momentum and toatal spin is returned
c
c     the output subroutine parameters are
c
c    eexc = excitation energy
c
c    kermsg = blank if no errors
c
c     written by j. j. smith , iaea atomic and molecular data unit
c
c------------------------------------------------------------------------
c
      character*(*) kermsg
      integer n, l, mult, sumen, lp, multp
      dimension enl(4,4,2), en(7)
      data (((enl(i,j,k),i=1,4),j=1,4),k=1,2) /32*0.0/
c
c----   n=1 state
c
c     data enl(1,1,1)
c    *    / 24.588 /,
c
c----   n=2 excited states
c
c    1   enl(2,1,1),enl(2,1,2),enl(2,2,1),enl(2,2,2)
c    *    / 20.614, 19.818, 21.217, 20.963 /,
c
c----   n=3 excited states
c
c    2   enl(3,1,1),enl(3,1,2),enl(3,2,1),enl(3,2,2),enl(3,3,1),
c    *   enl(3,3,2)
c    *    / 22.919, 22.717, 23.086, 23.006, 23.073, 23.072 /,
c
c----   n=4 excited states
c
c    3   enl(4,1,1),enl(4,1,2),enl(4,2,1),enl(4,2,2),enl(4,3,1),
c    *   enl(4,3,2),enl(4,4,1),enl(4,4,2)
c    *    / 23.672, 23.529, 23.741, 23.706, 23.735, 23.735, 23.736,
c    *     23.736 /
c
c-----   excitation energies summed over angular momentum and
c-----   total spin
c
      enl(1,1,1)=24.588
c----   n=2 excited states
      enl(2,1,1)=20.614
      enl(2,1,2)=19.818
      enl(2,2,1)=21.217
      enl(2,2,2)=20.963
c----   n=3 excited states
      enl(3,1,1)=22.919
      enl(3,1,2)=22.717
      enl(3,2,1)=23.086
      enl(3,2,2)=23.006
      enl(3,3,1)=23.073
      enl(3,3,2)=23.072
c----   n=4 excited states
      enl(4,1,1)=23.672
      enl(4,1,2)=23.529
      enl(4,2,1)=23.741
      enl(4,2,2)=23.706
      enl(4,3,1)=23.735
      enl(4,3,2)=23.735
      enl(4,4,1)=23.736
      enl(4,4,2)=23.736
      data (en(k),k=1,2) /2*0.0/
      data en(3),en(4),en(5),en(6),en(7)
     *    / 22.9799 , 23.699, 24.07, 24.30, 24.71/
c
      kermsg =' '
      if (sumen .eq. 1 .or. n .gt. 4) then
        if (n .ge. 8) then
          eexc =enl(1,1,1) - 13.58/(n*n)
          return
        else
           eexc =en(n)
           if (eexc .eq. 0.0) kermsg =
     *        'excitation energy for n value not in table in heexcen'
        endif
      else
           lp = l+1
           if (mult .eq. 1) then
             multp = 1
           else if (mult .eq. 3) then
               multp = 2
             else
                 kermsg = 'invalid spin multipclity as input in heexcen'
             endif
           eexc = enl(n,lp,multp)
           if (eexc .eq. 0.0) kermsg =
     *        'excitation energy not in table in heexcen'
      endif
      return
      end
c##################################################################
c###################################################################
c
      subroutine oscsthe(nin, lin, nfin, lfin, mult, j, fnm, kermsg)
c
c     this subroutine passes the oscillator srtrength for helium for
c     allowed transitions taken from the table of transitions
c     given by janev et al, see doc=h-he-plasma.
c
c     the subroutine parameters are
c
c    nin = principal quantum number of the initial state
c
c    lin = orbital angular momentum quantum number of the initial state
c
c    nfin = principal quantum number of the final state
c
c    lfin = orbital angular momentum quantum number of the final state
c
c    mult = the spin multiplicity (2s+1) of the initial state
c
c    j  = the total quantum number of the initial state
c         (required for fine structure transitions only)
c
c    fnm = oscillator strength
c
c    kermsg = blank if no errors
c
c     written by j. j. smith , iaea atomic and molecular data unit
c
c------------------------------------------------------------------------
c
      character*(*) kermsg
      dimension tfnm(2,2,5,3,2)
      dimension t2s2p(4)
c
      data (((((tfnm(i,j,k,l,m),i=1,2),j=1,2),k=1,5),l=1,3),m=1,2)
     &      /120*0.0/
c
c-----   singlet transitions 1s1s - 1snp
c
c     data tfnm(1,1,2,2,1),tfnm(1,1,3,2,1),tfnm(1,1,4,2,1)
c    *     , tfnm(1,1,5,2,1)
c    *    / 0.276, 7.34e-02, 3.02e-02, 1.53e-02 /,
c
c-----   singlet transitions 1s2s - 1snp
c
c    1   tfnm(2,1,2,2,1),tfnm(2,1,3,2,1),tfnm(2,1,4,2,1),tfnm(2,1,5,2,1)
c    *    / 0.376, 0.151, 0.0507, 0.0221 /,
c
c-----  triplet transitions 1s2s - 1s2p  fine structure - in array t2s2p
c
c-----   triplet transitions 1s2s - 1snp
c
c    2   tfnm(2,1,3,2,2),tfnm(2,1,4,2,2),tfnm(2,1,5,2,2)
c    *    / 0.0645, 0.0231, 0.0114 /,
c
c-----   singlet transitions 1s2p - 1sns
c
c    3   tfnm(2,2,3,1,1),tfnm(2,2,4,1,1),tfnm(2,2,5,1,1)
c    *    / 0.0480, 0.834e-02, 0.308e-02 /,
c
c-----   singlet transitions 1s2p - 1snd
c
c    4   tfnm(2,2,3,3,1),tfnm(2,2,4,3,1),tfnm(2,2,5,3,1)
c    *    / 0.711, 0.122, 0.0436 /,
c
c-----   triplet transitions 1s2p - 1sns
c
c    5   tfnm(2,2,3,1,2),tfnm(2,2,4,1,2),tfnm(2,2,5,1,2)
c    *    / 0.0692, 0.0118, 3.65e-02 /,
c
c-----   transitions 1s2pd - 1snd
c
c    6   tfnm(2,2,3,3,2),tfnm(2,2,4,3,2),tfnm(2,2,5,3,2)
c    *    / 0.609, 0.125, 0.0474 /
c
c-----   triplet fine structure transitions  1s2s-1s2p
c-----   the first value is an average over the fine structure levels
c-----   the following three are for j = 0, 1, 2
c
      data t2s2p(1),t2s2p(2),t2s2p(3),t2s2p(4)/0.539,0.300,0.180,0.060/
c
c---  check for fine structure transitions 1s2s - 1s2p
c---  if j value is >= 0 the fine structure oscillator strength is taken
c---  if j value is < 0 the avearge value is taken.
c
c-----   singlet transitions 1s1s - 1snp
      tfnm(1,1,2,2,1)=0.276
      tfnm(1,1,3,2,1)=7.34e-02
      tfnm(1,1,4,2,1)=3.02e-02
      tfnm(1,1,5,2,1)=1.53e-02
c-----   singlet transitions 1s2s - 1snp
      tfnm(2,1,2,2,1)=0.376
      tfnm(2,1,3,2,1)=0.151
      tfnm(2,1,4,2,1)=0.0507
      tfnm(2,1,5,2,1)=0.0221
c-----   triplet transitions 1s2s - 1snp
      tfnm(2,1,3,2,2)=0.0645
      tfnm(2,1,4,2,2)=0.0231
      tfnm(2,1,5,2,2)=0.0114
c-----   singlet transitions 1s2p - 1sns
      tfnm(2,2,3,1,1)=0.0480
      tfnm(2,2,4,1,1)=0.834e-02
      tfnm(2,2,5,1,1)=0.308e-02
c-----   singlet transitions 1s2p - 1snd
      tfnm(2,2,3,3,1)=0.711
      tfnm(2,2,4,3,1)=0.122
      tfnm(2,2,5,3,1)=0.0436
c-----   triplet transitions 1s2p - 1sns
      tfnm(2,2,3,1,2)=0.0692
      tfnm(2,2,4,1,2)=0.0118
      tfnm(2,2,5,1,2)=3.65e-02
c-----   transitions 1s2pd - 1snd
      tfnm(2,2,3,3,2)=0.609
      tfnm(2,2,4,3,2)=0.125
      tfnm(2,2,5,3,2)=0.0474
      kermsg = ' '
      if ( (nin .eq. 2) .and. (lin .eq. 0)  .and. (nfin .eq. 2)
     *    .and. (lfin .le. 1) .and. (mult .eq. 3)) then
        if ( (j .ge. 0) .and. (j .le. 2)) then
          fnm=t2s2p(j+2)
        else
          fnm=t2s2p(1)
        endif
 
      else if ( (nin .le. 2) .and. (lin .le. 1) .and. (nfin .le. 5)
     *                       .and. (lfin .le. 2) ) then
c
             if (mult .eq. 3) then
               imult = 2
             else
               imult = 1
             endif
c
             fnm = tfnm (nin, lin+1, nfin, lfin+1, imult)
             if (fnm .eq. 0.0) then
               kermsg ='transition not included in table in oscsthe'
             endif
      else
           kermsg ='transition not included in table in oscsthe'
           fnm=0.0
      endif
      return
      end

c
c###################################################################
c
      subroutine oscsth2(n, m, fnm)
c
c     this function calculates the oscillator strength fnm for
c     hydrogen from the formulae of johnson , l. c., astrophys. j. 174,
c     227, (1972).
c
c     n = the principle quantum number of the initial state
c     m = the principle quantum number of the final state
c
      data pi /3.14159/
c
      an = n
      am = m
      n2= an * an
      if (n. eq. 1) then
        g0 = 1.1330
        g1 = -0.4059
        g2 = 0.07014
      else if(n .eq. 2) then
        g0 = 1.0785
        g1 = -0.2319
        g2 = 0.02947
      else
        g0 = 0.9935 + 0.2328/an  -  0.1296/n2
        g1 = -(0.6282 - 0.5598/an + 0.5299/n2)/an
        g2 = (0.3887 - 1.181/an + 1.47/n2)/n2
      endif
c
      x = 1.0 -((an/am)**2)
      gnx = g0 + g1/x + g2/(x*x)
      fnm = 32.0*an*gnx /(3.0*sqrt(3.0)*pi*(am**3.0)*(x**3.0))
c
      return
      end
c
c###################################################################
c
      function dbeta (beta)
c
c     this function calculates the d(beta) function of janev, r. k.,
c     and presnyakov, l. p., j. phys. b, 13, 4233, (1980)
c
c------------------------------------------------------------------------
c
      dimension binvt (36), dbetat (36)
c
c-----  table of 1/beta
c
      data binvt
     & /  0.1,    0.2,    0.3,    0.4,    0.5,    0.6,    0.8,    1.0,
     &    1.25,   1.5,    2.0,    2.5,    3.0,    3.5,    4.0,    4.5,
     &    5.0,    6.0,    8.0,   10.0,   12.5,   15.0,   20.0,   25.0,
     &   30.0,   40.0,   60.0,   80.0,  100.0,  150.0,  200.0,  300.0,
     &  400.0,  600.0,  800.0, 1000.0 /
c
c-----  table of function d(beta)
c
      data dbetat
     & / 0.057, 0.104,  0.135,  0.157,  0.175,  0.194,  0.230, 0.264,
     &   0.296, 0.328,  0.367,  0.388,  0.399,  0.405,  0.410, 0.405,
     &   0.399, 0.380,  0.345,  0.318,  0.285,  0.263,  0.227, 0.205,
     &   0.190, 0.168,  0.141,  0.124,  0.110,  0.092,  0.080, 0.064,
     &   0.054, 0.042,  0.035,  0.0289 /
c
      if (beta .lt. 1.0e-03) then
        dbeta = 4.0 * beta * log (1.4/beta)
c
        else if (beta .gt. 1.0e+01) then
          dbeta = beta * exp (-sqrt (2.0*beta) ) / 2.0
          else
c
c--   determine d(beta) using a linear interpolation from values in
c--   the data table  dbetat
c
            binv = 1.0/beta
            do 10 i = 1, 35
              if (binv .ge. binvt(i) .and. binv .le. binvt(i+1)) then
                dbeta =  ( (binvt(i+1) - binv) * dbetat(i) +
     &                     ( binv - binvt(i)) * dbetat (i+1) )
     &                    / ( binvt(i+1) - binvt(i) )
              endif
 10         continue
         endif
c
         return
         end
c###################################################################
c
      subroutine pexch4(pe, pcf, kncf, pxs, kermsg)
c
c     this is a subroutine to calculate cross sections (cm[2])
c     versus energy (ev) for heavy particle impact excitation
c     in the dipole-approximation close-coupling method (dacc).
c
c     pe = collision energy in ev
c
c     the coefficient data array passed should contain
c
c     pcf(1)  =   zt, the atomic number of target
c
c     pcf(2)  =   zp, the atomic number of projectile
c
c     pcf(3)  =  itrans,  integer which defines the type of transition.
c                 for transitions defined only in terms of the initial
c                 and final principal quantum numbers (n,m), itrans=1.
c                 for tranitions between (nl,ml') states itrans=2
c
c     pcf(4) =   n, the principal quantum number of the initial state
c
c     pcf(5) =   m, the principal quantum number of the final state
c
c     if itrans=2, the remainder of the coefficient array should
c                 contain
c
c     pcf(6) =   sumen, this indicates choice of ionization energy to be
c                returned.  if sumen=0  the ionization energy for the
c                specific state (quantified by n, l amd mult) is
c                returned.  if sumen=1  the ionization energy taken as
c                an average over angular momentum and
c                total spin is returned
c
c
c     pcf(7) =   lin, the orbital angular monentum of the inital state
c
c     pcf(8) =   lfin, the orbital angular monentum of the final state
c
c     pcf(9) =   mult, spin multiplicity (2s+1) of the inital state
c
c     pcf(10) =   j, the total angular momentum of the initial state
c                only required for fine structure transitions
c
c    - warning- .
c
c        the coefficient array pcf is updated by this routine to
c     include energy independent constants. these coefficients can be
c     used in subsequent calls for the same entry. the coefficeients
c     added are:
c
c     pcf(11) = eth,threshold excitation energy for the transition (ev)
c
c     pcf(12) = wnm, energy difference between initial and final states
c
c     pcf(13) = parameter lambda
c
c     kermsg = blank if no errors
c
c     pxs = cross section in cm[2]
c
c      this evaluation function can only derive cross sections for those
c    transitions for which the oscillator strengths have been included
c    in subroutines oscsth1 and oscthe. oscsth1 contains the oscillator
c    strengths for allowed transitions for hydrogen from table a.1. of
c    janev et al. oscthe contains the oscillator strengths for allowed
c    transitions for helium  from table a.4 of janev et al.
c      for transitions beteen excited states of helium the excitation
c    energies are returned by subroutine heexcen and the ionization
c    energies are returned by subroutine heionen.
c
c------------------------------------------------------------------------
c
      double precision pe, pcf, pxs
c
      integer zt, zp, n, m, sumen, itrans
      real lambda
      dimension pcf(13)
      character*(*) kermsg
c
      data ry /1.358e+01/
c
      zt = pcf(1)
      zp = pcf(2)
      rzt = zt
      rzp = zp
      itrans  = pcf(3)
      n = pcf(4)
      m = pcf(5)
      an = n
      am = m
c
      if (m .le. n) then
        kermsg = 'error  principal quantum number m <= n '
        return
      endif
c
      if (kncf .lt. 11) then
c
c        first call to pexch4 determine energy independent
c        parameters and place in pcf for further use
c
c        (n,m) transitions - determine energies, wnm and oscillator
c                            strength.
c
        if (itrans .eq. 1) then
          if (zt .eq. 1) then
            c= 1.0/(an*an) - 1.0/(am*am)
            eth = ry * c
            wnm = 0.5 * c
          else if (zt .eq. 2) then
            call heionen(n, 0, 0, 1, eionn, kermsg)
            call heionen(m, 0, 0, 1, eionm, kermsg)
            eth = eionn - eionm
            wnm = eth /(2.0*ry)
          else
            kermsg = 'target not covered by pexch4 '
            return
          endif
c
c     calculate the oscillator strength, fnm, for the transition from
c     the formula of johnson
c
          if (zt .eq. 1) then
            call oscsth2(n, m, fnm)
          else
            kermsg = 'average oscillator strength missing in pexch4'
            return
          endif
c
c        (nl,ml') transitions - determine energies, wnm and oscillator
c                               strength.
c
        else if (itrans .eq. 2 .and. zt .eq. 2) then
          if (kncf .ne. 10) then
            kermsg =
     1           ' invalid number of coefficients passed to pexch4'
            return
          endif
c
          sumen =  pcf(6)
          lin =  pcf(7)
          lfin = pcf(8)
          mult = pcf(9)
          j = pcf(10)
c
          if (n .eq .1) then
             call heexcen(m, lfin, mult, sumen, eth, kermsg)
             wnm = eth /(2.0*ry)
          else
              call heionen(n, lin, mult, sumen, ethn, kermsg)
              call heionen(m, lfin, mult, sumen, ethm, kermsg)
              eth = ethn - ethm
              wnm = eth /(2.0*ry)
          endif
        endif
c
        if (itrans .eq. 2) then
          call oscsthe(n, lin, m, lfin, mult, j, fnm, kermsg)
        endif
c
c
c     determine the potential strength, lambda
c
      lambda= sqrt(fnm/(2.0*wnm))
c
c        place energy independent parameters in coefficient array and
c        update kncf
c
        pcf(11) = eth
        pcf(12) = wnm
        pcf(13) = lambda
        kncf = 13
      else if (kncf .eq. 13) then
          eth  = pcf(11)
          wnm  = pcf(12)
          lambda = pcf(13)
c
      else
          kermsg = ' incorrect number of coefficients passed to pexch4'
          return
      endif
c
      if(pe .lt. eth) then
        pxs=0.0
        return
      endif
c
c     determine the velocity, v
c     corrected 1/zp factor -- dh -- 13 March 2008
        v = 6.3246e-03* sqrt(pe) / zp
c
c     determine the value of beta
c
      beta = zp* lambda *  wnm / (v * v)
c
c     determine the value of the cross section  pxs
c
      pxs = 1.76e-16 *  zp * lambda * dbeta(beta) / wnm
c
      return
      end

def heionen(n, l, mult, sumen):
    """
    This function passes the ionization energy for excited states of helium 
    taken from tables a.4 and a.5 given by janev et al.
 
    n: principal quantum number of the excited electronic state
    l: orbital angular momentum quantum number of the exited state
    mult: the spin multiplicity (2s+1) of the state
    sumen: indicates choice of ionization energy to be returned.
    """
    enl = np.zeros((4, 4, 2))
    en = np.zeros(7)
    
    # Define the ionization energy values
    enl[0, 0, 0] = 24.588
    enl[1, 0, 0] = 3.973
    enl[1, 0, 1] = 4.769
    enl[1, 1, 0] = 3.371
    enl[1, 1, 1] = 3.625
    enl[2, 0, 0] = 1.699
    enl[2, 0, 1] = 1.871
    enl[2, 1, 0] = 1.502
    enl[2, 1, 1] = 1.582
    enl[2, 2, 0] = 1.515
    enl[2, 2, 1] = 1.515
    enl[3, 0, 0] = 0.9155
    enl[3, 0, 1] = 0.9951
    enl[3, 1, 0] = 0.8470
    enl[3, 1, 1] = 0.8812
    enl[3, 2, 0] = 0.8529
    enl[3, 2, 1] = 0.8530
    enl[3, 3, 0] = 0.8520
    enl[3, 3, 1] = 0.8520
    en[2:5] = [1.609, 0.8811, 0.52, 0.37, 0.29]
    
    eion = 0.0
    
    if sumen == 1:
        if n >= 8:
            eion = 13.58 / (n * n)
        else:
            try:
                eion = en[n - 1]
            except IndexError:
                raise ValueError('Ionization energy for n value not in table in heionen')
    else:
        lp = l + 1
        if mult == 1:
            multp = 1
        elif mult == 3:
            multp = 2
        else:
            raise ValueError('Invalid spin multiplicity as input in heionen')
        
        try:
            eion = enl[n - 1, lp - 1, multp - 1]
        except IndexError:
            raise ValueError('Ionization energy not in table in heionen')
    
    return eion
##############################################################

def heexcen(n, l, mult, sumen):
   """
    This function passes the excitation energy for excited states of  helium 
    taken from tables a.4 and a.5 given by janev et al. (Elementary processes in H-He plasmsas.
 
    n: principal quantum number of the excited electronic state
    l: orbital angular momentum quantum number of the exited state
    mul: the spin multiplicity (2s+1) of the state
    sumen: indicates choice of excitation enegy to be returned.
    """
    enl = np.zeros((4, 4, 2))
    en = np.zeros(7)
    
    # Define the excitation energy values
    enl[0, 0, 0] = 24.588
    enl[1, 0, 0] = 20.614
    enl[1, 0, 1] = 19.818
    enl[1, 1, 0] = 21.217
    enl[1, 1, 1] = 20.963
    enl[2, 0, 0] = 22.919
    enl[2, 0, 1] = 22.717
    enl[2, 1, 0] = 23.086
    enl[2, 1, 1] = 23.006
    enl[2, 2, 0] = 23.073
    enl[2, 2, 1] = 23.072
    enl[3, 0, 0] = 23.672
    enl[3, 0, 1] = 23.529
    enl[3, 1, 0] = 23.741
    enl[3, 1, 1] = 23.706
    enl[3, 2, 0] = 23.735
    enl[3, 2, 1] = 23.735
    enl[3, 3, 0] = 23.736
    enl[3, 3, 1] = 23.736
    en[2:5] = [22.9799, 23.699, 24.07, 24.30, 24.71]
    
    if sumen == 1 or n > 4:
        if n >= 8:
            eexc = enl[0, 0, 0] - 13.58 / (n * n)
        else:
            try:
                eexc = en[n - 1]
            except IndexError:
                raise ValueError('Excitation energy for n value not in table in heexcen')
    else:
        lp = l + 1
        if mult == 1:
            multp = 1
        elif mult == 3:
            multp = 2
        else:
            raise ValueError('Invalid spin multiplicity as input in heexcen')
        
        try:
            eexc = enl[n - 1, lp - 1, multp - 1]
        except IndexError:
            raise ValueError('Excitation energy not in table in heexcen')
    
    return eexc

##############################################################

def oscsthe(nin, lin, nfin, lfin, mult, j):
    """
    This function calculates the oscillator strength for helium allowed transitions.

    nin: principal quantum number of the initial state
    lin: orbital angular momentum quantum number of the initial state
    nfin: principal quantum number of the final state
    lfin: orbital angular momentum quantum number of the final state
    mult: the spin multiplicity (2s+1) of the initial state
    j: the total quantum number of the initial state (required for fine structure transitions only)

    Returns:
    fnm: oscillator strength
    """
    tfnm = np.zeros((2, 2, 5, 3, 2))
    t2s2p = np.array([0.539, 0.300, 0.180, 0.060])

    # Initialize tfnm array with zeros
    tfnm[:] = 0.0

    # Set the relevant values in tfnm array
    tfnm[0, 0, 1, 1, 0] = 0.276
    tfnm[0, 0, 2, 1, 0] = 0.0734
    tfnm[0, 0, 3, 1, 0] = 0.0302
    tfnm[0, 0, 4, 1, 0] = 0.0153

    tfnm[1, 0, 1, 1, 0] = 0.376
    tfnm[1, 0, 2, 1, 0] = 0.151
    tfnm[1, 0, 3, 1, 0] = 0.0507
    tfnm[1, 0, 4, 1, 0] = 0.0221

    tfnm[1, 0, 2, 1, 1] = 0.0645
    tfnm[1, 0, 3, 1, 1] = 0.0231
    tfnm[1, 0, 4, 1, 1] = 0.0114

    tfnm[1, 1, 2, 0, 0] = 0.0480
    tfnm[1, 1, 3, 0, 0] = 0.00834
    tfnm[1, 1, 4, 0, 0] = 0.00308

    tfnm[1, 1, 2, 2, 0] = 0.711
    tfnm[1, 1, 3, 2, 0] = 0.122
    tfnm[1, 1, 4, 2, 0] = 0.0436

    tfnm[1, 1, 2, 0, 1] = 0.0692
    tfnm[1, 1, 3, 0, 1] = 0.0118
    tfnm[1, 1, 4, 0, 1] = 0.0365

    tfnm[1, 1, 2, 2, 1] = 0.609
    tfnm[1, 1, 3, 2, 1] = 0.125
    tfnm[1, 1, 4, 2, 1] = 0.0474

    fnm = 0.0  # Default value
    imult = 0  # Default value

    if (nin == 2 and lin == 0 and nfin == 2 and lfin <= 1 and mult == 3):
        if (j >= 0 and j <= 2):
            fnm = t2s2p[j + 2]
        else:
            fnm = t2s2p[1]
    elif (nin <= 2 and lin <= 1 and nfin <= 5 and lfin <= 2):
        if (mult == 3):
            imult = 2
        else:
            imult = 1

        fnm = tfnm[nin-1, lin, nfin-1, lfin, imult]

        if (fnm == 0.0):
            raise ValueError("Transition not included in table in oscsthe")
    else:
        raise ValueError("Transition not included in table in oscsthe")
        fnm = 0.0

    return fnm

##############################################################
def oscsth2(n, m):
    """ 
    This function calculates the oscillator strength fnm for hydrogen 
    using the formulae of Johnson, L. C., Astrophys. J. 174, 227, (1972).
    """

    pi = 3.14159
    
    an = n
    am = m
    n2 = an * an

    if n == 1:
        g0 = 1.1330
        g1 = -0.4059
        g2 = 0.07014
    elif n == 2:
        g0 = 1.0785
        g1 = -0.2319
        g2 = 0.02947
    else:
        g0 = 0.9935 + 0.2328 / an - 0.1296 / n2
        g1 = -(0.6282 - 0.5598 / an + 0.5299 / n2) / an
        g2 = (0.3887 - 1.181 / an + 1.47 / n2) / n2
    
    x = 1.0 - ((an / am) ** 2)
    gnx = g0 + g1 / x + g2 / (x * x)
    fnm = 32.0 * an * gnx / (3.0 * np.sqrt(3.0) * pi * (am ** 3.0) * (x ** 3.0))
    
    return fnm

##############################################################
def dbeta(beta):
    """
    this function calculates the d(beta) function of janev, r. k.,
     and presnyakov, l. p., j. phys. b, 13, 4233, (1980)
    """
    binvt = np.array([
        0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.8, 1.0,
        1.25, 1.5, 2.0, 2.5, 3.0, 3.5, 4.0, 4.5,
        5.0, 6.0, 8.0, 10.0, 12.5, 15.0, 20.0, 25.0,
        30.0, 40.0, 60.0, 80.0, 100.0, 150.0, 200.0, 300.0,
        400.0, 600.0, 800.0, 1000.0
    ])
    
    dbetat = np.array([
        0.057, 0.104, 0.135, 0.157, 0.175, 0.194, 0.230, 0.264,
        0.296, 0.328, 0.367, 0.388, 0.399, 0.405, 0.410, 0.405,
        0.399, 0.380, 0.345, 0.318, 0.285, 0.263, 0.227, 0.205,
        0.190, 0.168, 0.141, 0.124, 0.110, 0.092, 0.080, 0.064,
        0.054, 0.042, 0.035, 0.0289
    ])
    
    if beta < 1.0e-03:
        dbeta = 4.0 * beta * np.log(1.4 / beta)
    elif beta > 1.0e+01:
        dbeta = beta * np.exp(-np.sqrt(2.0 * beta)) / 2.0
    else:
        binv = 1.0 / beta
        for i in range(35):
            if binv >= binvt[i] and binv <= binvt[i + 1]:
                dbeta = ((binvt[i + 1] - binv) * dbetat[i] +
                         (binv - binvt[i]) * dbetat[i + 1]) / (binvt[i + 1] - binvt[i])
                break
    
    return dbeta

##############################################################

def pexch4(pe, pcf, kncf):
    """
    This function calculates  calculate cross sections (cm[2]) versus energy (ev) for heavy particle 
     impact excitation in the dipole-approximation close-coupling method (dacc).
    
    pe: collision energy in eV
    pcf: parameters for fit to the cross section
    """
    zt = pcf[0]
    zp = pcf[1]
    rzt = zt
    rzp = zp
    itrans = pcf[2]
    n = pcf[3]
    m = pcf[4]
    an = n
    am = m
    
    ry = 1.358e+01
    
    if m <= n:
        raise ValueError("Principal quantum number m <= n")
    
    if kncf < 11:
        if itrans == 1:
            if zt == 1:
                c = 1.0 / (an * an) - 1.0 / (am * am)
                eth = ry * c
                wnm = 0.5 * c
            elif zt == 2:
                eionn = heionen(n, 0, 0, 1)  # heionen function call
                eionm = heionen(m, 0, 0, 1)  # heionen function call
                eth = eionn - eionm
                wnm = eth / (2.0 * ry)
            else:
                raise ValueError("Target not covered by pexch4")
            
            if zt == 1:
                fnm = oscsth2(n, m)  # oscsth2 function call
            else:
                raise ValueError("Average oscillator strength missing in pexch4")
            
        elif itrans == 2 and zt == 2:
            if kncf != 10:
                raise ValueError("Invalid number of coefficients passed to pexch4")
            
            sumen = pcf[5]
            lin = pcf[6]
            lfin = pcf[7]
            mult = pcf[8]
            j = pcf[9]
            
            if n == 1:
                eth = heexcen(m, lfin, mult, sumen)  #heexcen function call
                wnm = eth / (2.0 * ry)
            else:
                ethn = heionen(n, lin, mult, sumen)  # heionen function call
                ethm = heionen(m, lfin, mult, sumen)  # heionen function call
                eth = ethn - ethm
                wnm = eth / (2.0 * ry)
            
            fnm = oscsthe(n, lin, m, lfin, mult, j)  #  oscsthe function call
        
        lambda_val = np.sqrt(fnm / (2.0 * wnm))
        
        pcf[10] = eth
        pcf[11] = wnm
        pcf[12] = lambda_val
        kncf = 13
        
    elif kncf == 13:
        eth = pcf[10]
        wnm = pcf[11]
        lambda_val = pcf[12]
    
    else:
        raise ValueError("Incorrect number of coefficients passed to pexch4")
    
    if pe < eth:
        pxs = 0.0
        return pxs
    
    v = 6.3246e-03 * np.sqrt(pe) / zp
    beta = zp * lambda_val * wnm / (v * v)
    pxs = 1.76e-16 * zp * lambda_val * dbeta(beta) / wnm
    
    return pxs