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###################################################################
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
      subroutine jeexc2(pe, pcf, kncf, pxs, kermsg)
c
c     this is a subroutine to calculate a cross section
c     in cm[2] as function of electron energy in ev for dipole
c     forbidden transitions without a spin change.
c     for details see doc=h-he-plasma , used for reactions 2.3.2
c     dipole-forbidden transitions only
c
c     pcf is the coefficient data array, where
c
c     pcf(1)  =   coefficient beta
c
c     pcf(2)  =   coefficient gamma
c
c     pcf(3)  =   coefficient delta
c
c     pcf(4)  =   n, the principal quantum number of the final state
c
c     pcf(5)  =   l, the angular momentum of the the final state
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(6) = threshold excitation energy for the transition (ev)
c
c     pcf(7) = ionization energy of the state
c
c     pcf(8) = the coefficient t (see doc=h-he-plasma for definition)
c
c     pe = electron energy (ev)
c
c     kermsg = blank if no errors
c
c     pxs = cross section in cm[2]
c
c     written by j. j. smith , iaea atomic and molecular data unit
c
c------------------------------------------------------------------------
c
      double precision pe, pcf, pxs
      integer n, l
      dimension pcf(8)
      dimension b(2,3)
      character*(*) kermsg
      data ry/13.58/
      data eionhe/24.5876/
c
c         the following table contains the values for the constant
c         b(1s,nl) taken from the papers of kingston et al
c         proc. phys. soc. 87,399 (1966), ibid 88, 597 (1966)
c         (for helium the value of b must be divided by two to account
c         for the two electrons, this is carried in using the b table )
c
c----   n=2 excited states
c
      data b(1,1),b(1,2)/2.21972e-01,1.31230/,
c
c----   n=3 excited states
c
     1   b(2,1),b(2,2),b(2,3)/4.40709e-02,2.43642e-01,3.39007e-02/

      beta = pcf(1)
      gamma = pcf(2)
      delta = pcf(3)
      n = pcf(4)
      l = pcf(5)
c
      kermsg = ' '
      if (kncf .lt. 6) then
c
c        first call to jeexc2 determine energy independent
c        parameters and place in pcf for further use
c
c        determine the excitation energy of the final state
c
        call heexcen (n, l,  1, 0,  eexc, kermsg)
        if (kermsg .ne.  ' ') return
c
c        determine the ionization energy of the final state
c
         call heionen (n, l,  1, 0,  eion, kermsg)
         if ( kermsg .ne.  ' ') return
c
c        determine the parameter t
c
         t = 1.6 * beta * ((b(n,l+1)/2.0)**(-gamma))
     &        * ( ( ry * eexc / (4.0*eionhe*eion) ) ** (1.0 - gamma))
c
c        place energy independent parameters in coefficient array and
c        update kncf
c
         pcf(6) = eexc
         pcf(7) = eion
         pcf(8) = t
         kncf = 8
c
      else if (kncf .eq. 8) then
          eexc = pcf(6)
          eion = pcf(7)
          t   = pcf(8)
c
      else
          kermsg = ' incorrect number of coefficients passed to jeexc2'
          return
      endif
c
      if (pe .lt. eexc) then
        pxs=0.0
        return
      endif
c
      u = pe / eexc
      cexp = -t * u
      if (cexp .lt. -20.0) then
        vexp = 0.0
      else
        vexp = exp(cexp)
      endif
c
      y  = (1.0 - vexp) * (u - 1.0 + delta)/ (u*u)
      pxs = 3.52e-16 * (b(n,l+1)/2.0) * (ry/eexc) * y
c
      return
c
      end

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 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 jeexc2(pe, pcf, kncf):
    """
    This function is used to calculate cross section in cm^2 as a function of 
    electron energy in eV for dipole forbidden transitions without a spin change.    
         
    pe = electron energy (eV)
    pcf is the coefficient data array, where
    pcf[0] = coefficient beta
    pcf[1] = coefficient gamma
    pcf[2] = coefficient delta
    pcf[3] = n, the principal quantum number of the final state
    pcf[4] = l, the angular momentum of the final state
    """

    beta = pcf[0]
    gamma = pcf[1]
    delta = pcf[2]
    n = pcf[3]
    l = pcf[4]

    ry = 13.58
    eionhe = 24.5876

    b = np.array([[2.21972e-01, 1.31230],  # n=2 excited states
                  [4.40709e-02, 2.43642e-01, 3.39007e-02]])  # n=3 excited states

    if kncf < 6:
        # First call to jeexc2 determines energy-independent parameters and places 
        #them in pcf for further use

        # Determine the excitation energy of the final state
        eexc = heexcen(n, l, 1, 0)
        if eexc is None:
            raise ValueError("Error: Unable to calculate excitation energy.")

        # Determine the ionization energy of the final state
        eion = heionen(n, l, 1, 0)
        if eion is None:
            raise ValueError("Error: Unable to calculate ionization energy.")

        # Determine the parameter t
        t = 1.6 * beta * ((b[n-2, l] / 2.0) ** (-gamma)) * ((ry * eexc / (4.0 * eionhe * eion)) ** (1.0 - gamma))

        # Place energy-independent parameters in the coefficient array and update kncf
        pcf[5] = eexc
        pcf[6] = eion
        pcf[7] = t
        kncf = 8

    elif kncf == 8:
        eexc = pcf[5]
        eion = pcf[6]
        t = pcf[7]
    else:
        raise ValueError("Error: Incorrect number of coefficients passed to jeexc2.")

    # Check for impact energies below threshold
    if pe < eexc:
        pxs = 0.0
        return pxs

    u = pe / eexc

    cexp = -t * u
    vexp = np.where(cexp < -20.0, 0.0, np.exp(cexp))

    y = (1.0 - vexp) * (u - 1.0 + delta) / (u ** 2)
    pxs = 3.52e-16 * (b[n-2, l] / 2.0) * (ry / eexc) * y

    return pxs