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
      subroutine jeexc1(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     allowed transitions.
c     for details see doc=h-he-plasma , used for reactions 2.3.1,
c     2.3.5a, 2.3.5b, 2.3.6a, 2.3.6b, and 2.3.7.
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)  = nin(n), principal quantum number of the initital state
c
c     pcf(5)  = lin(l), orbital angular momentum of the initital state
c
c     pcf(6)  = lfin(l') orbital angular momentum of the final state
c
c     pcf(7)  = mult, spin multiplicity (2s+1) of the inital state
c
c     pcf(8) =   j, the total angular momentum of the initial state
c                only required for fine structure transitions
c                (see routine oscsthe for details of definition)
c
c     pcf(9) = nfin(n'), principal quantum number of the initital 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(10) = threshold energy for the transition (ev)
c
c     pcf(11) = oscillator strength for the transition
c
c     pcf(12) = the coefficient xi (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 nin, lin, nfin, lfin, mult, j
      dimension pcf(12)
      character*(*) kermsg
      data ry/13.58/
c
      beta = pcf(1)
      gamma = pcf(2)
      delta = pcf(3)
      nin = pcf(4)
      lin = pcf(5)
      lfin = pcf(6)
      mult = pcf(7)
      j = pcf(8)
      nfin = pcf(9)
c
      kermsg = ' '
      if (kncf .lt. 10) then
c
c        first call to jeexc1 determine energy independent
c        parameters and place in pcf for further use
c
c        determine the excitation energy of the final state
c
         eexcl = 0.0
         if (nin .ne.  1) then
           call heexcen (nin, lin, mult, 0, eexcl, kermsg)
           if (kermsg .ne.  ' ') return
         endif
c
         call heexcen (nfin, lfin, mult, 0, eexcu, kermsg)
         if (kermsg .ne.  ' ') return
         eexc = eexcu - eexcl
c
c        determine the oscillator strength for the transition
c

c Yuri R.: there're some problems here
c Denis H.: Problem solved: size of array pcf passed to the function was not fixed

         call oscsthe(nin, lin, nfin, lfin, mult, j, fnm, kermsg)
         if (kermsg .ne. ' ') return
c
         xi = beta * ((ry * fnm /eexc)** (-gamma))
c
c        place energy independent parameters in coefficient array and
c        update kncf
c
         pcf(10) = eexc
         pcf(11) = fnm
         pcf(12) = xi
         kncf = 12
c
      else if (kncf .eq. 12) then
          eexc = pcf(10)
          fnm  = pcf(11)
          xi   = pcf(12)
c
      else
          kermsg = ' incorrect number of coefficients passed to jeexc1'
          return
      endif
c
c        check for impact energies below threshold
c
      if (pe. lt. eexc ) then
        pxs = 0.0
        return
      endif
c
      u = pe / eexc
c
      cexp =  - xi * (u + 1.0)
      if (cexp .lt. -20.0) then
        vexp = 0.0
      else
        vexp = exp(cexp)
      endif
c
      if (delta .eq. 0.0) then
         y  = (1.0 - vexp) * log (u) / u
      else
         y  = (1.0 - vexp) * log (u+delta) / u
      endif
c
      pxs = 3.52e-16 * ( (ry/eexc)**2 ) * fnm * 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 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 jeexc1(pe, pcf, kncf):
    """
    This function calculates the cross section in cm^2 as a function of 
    electron energy in eV for dipole allowed transitions.

    pe: electron energy (eV)
    pcf: coefficient data array
    kncf: number of coefficients

    Returns:
    pxs: cross section in cm^2
    """
    ry = 13.58  # Rydberg energy in eV

    beta, gamma, delta, nin, lin, lfin, mult, j, nfin = pcf

    if kncf < 10:
        eexcl = 0.0
        if nin != 1:
            eexcl = heexcen(nin, lin, mult, 0)  
            if np.isnan(eexcl):
                raise ValueError('Invalid heexcen result')
        eexcu = heexcen(nfin, lfin, mult, 0) 
        if np.isnan(eexcu):
            raise ValueError('Invalid heexcen result')
        eexc = eexcu - eexcl
        fnm = oscsthe(nin, lin, nfin, lfin, mult, j)  
        if np.isnan(fnm):
            raise ValueError('Invalid oscsthe result')

        xi = beta * ((ry * fnm / eexc) ** (-gamma))

        pcf[9] = eexc
        pcf[10] = fnm
        pcf[11] = xi
        kncf = 12
    elif kncf == 12:
        eexc, fnm, xi = pcf[9:12]
    else:
        raise ValueError('Incorrect number of coefficients passed to jeexc1')

    if pe < eexc:
        pxs = 0.0
        raise ValueError('Impact energy below threshold')
    else:
        u = pe / eexc
        cexp = -xi * (u + 1.0)
        vexp = np.where(cexp < -20.0, 0.0, np.exp(cexp))
    
        if delta == 0.0:
            y = (1.0 - vexp) * np.log(u) / u
        else:
            y = (1.0 - vexp) * np.log(u + delta) / u
        
        pxs = 3.52e-16 * ((ry / eexc) ** 2) * fnm * y
    
    return pxs