c
c###################################################################
c
      subroutine albeli(pe, pcf, kncf, pxs, kermsg)
c
c     electron impact ionization cross section fits
c
c     reference: k. l. bell et al, j. phys. chem. ref. data 12, 891
c                (1983)
c
c     this is an iaea subroutine to calculate cross sections for
c     projectile energy (ev).
c
c     pe =  electron energy (ev)
c
c     the number of fitting parameters varies depending on the
c     number of terms taken in the numerical fitting and on the
c     allowance for excitation autoionization in the cross section
c     to fit cross sections with excitation autoionization two seperate
c     fits are defined. one from the ionization threshold and a second
c     fit for energies above the autoionization threshold.
c     the number of parameters in any entry is given by kncf
c
c     pcf(1) = ionization potential (ev)
c     pcf(2-7) = fitting parameters ( can be less than 6 parameters)
c
c     if cross section has excitation autoionization structure then
c     for the second fit
c
c     pcf(8) = autoionization threshold (ev)
c     pcf(9) = ionization potential (ev)
c     pcf(10-15) = fitting parameters for this fit (can be less than
c                  6 parameters)
c
c     kncf = number of parameters supplied in pcf (must be 8)
c     pxs = ionization cross section (cm[2])
c     kermsg = error message, ' ' is ok
c
c     written by j. j. smith , iaea atomic and molecular data unit
c
c======================================================================
c
      double precision pe, pcf, pxs
      double precision ion, power, power1, xs, a, x, x2
      dimension pcf(15)
      character*(*) kermsg
c
      kermsg = ' '
      if(pe .lt. pcf(1)) then
        pxs = 0.0d0
        return
      else
        kermsg = ' '
      endif
c
c---  determine parameters to be used
c
      if (kncf .gt. 7 .and. pe .gt. pcf(8) ) then
c
c---      autoionization included and energy > autoionization threshold
c
        ion = pcf(9)
        a=pcf(10)
        istart=11
        iend=kncf
      else
        ion = pcf(1)
        a=pcf(2)
        istart=3
        if (kncf .gt. 7) then
          iend=7
        else
          iend = kncf
        endif
      endif
c
c---  generate cross section
c
      x=ion/pe
      x2= 1.0d0/x
c
c---  contribution from bethe term
c
      xs = a*dlog(x2)
      if ( kncf .ge. istart) then
c
c---  contribution from least squares fit terms
c
        power1 = 1.0d0 - x
        power = power1
        do 10 i=istart,iend
          xs = xs + pcf(i)*power
          power = power*power1
  10    continue
      endif
c
c---  scale results to cm[2]
c
      pxs= 1.0d-13*xs /(pe*ion)
c
      return
      end

def albeli(pe, pcf):
    """
    This function calculates the cross section for electron impact ionization.

    pe: electron energy in eV
    pcf: parameter data array
        pcf[0]: ionization potential (eV)
        pcf[1:7]: fitting parameters for the first fit reprsenting direct ionization
        if cross section has excitation autoionization structure then
        pcf[7]: autoionization threshold (eV)
        pcf[8]: ionization potential for the second fit (eV)
        pcf[9:15]: fitting parameters for the second fit
    """

    if pe < pcf[0]:
        warnings.warn('Electron energy is below the ionization potential, cross section set to 0.')
        pxs = 0.0
        return pxs

    if len(pcf) > 7 and pe > pcf[7]:
        # Autoionization included and energy > autoionization threshold
        ion = pcf[8]
        a = pcf[9]
        istart = 10
        iend = len(pcf)
    else:
        a = pcf[1]
        istart = 2
        if len(pcf) > 7:
            iend = 7
        else:
            iend = len(pcf)

    x = ion / pe
    x2 = 1.0 / x

    xs = a * np.log(x2)

    if len(pcf) >= istart:
        power1 = 1.0 - x
        power = power1
        for i in range(istart, iend+1):
            xs += pcf[i] * power
            power *= power1

    pxs = 1.0e-13 * xs / (pe * ion)

    return pxs