1. E /eV
2. rate /cm³ s^-1

• B53: Y. Itikawa, S. Hara, T. Kato, S. Nakazaki, M. Pindzola, D. Crandall, "Electron-impact cross sections and rate coefficients for excitations of carbon and oxygen ions", Atomic Data and Nuclear Data Tables 33, 149-193 (1985). [10.1016/0092-640x(85)90025-7]

def eiexp(x):
    x2 = x * x
    x3 = x2 * x
    x4 = x3 * x
    x5 = x4 * x

    if x <= 1.0:
        a0 = -0.57721566
        a1 = 0.99999193
        a2 = -0.24991055
        a3 = 0.05519968
        a4 = -0.00976004
        a5 = 0.00107857
        eiexp = (a0 + a1 * x + a2 * x2 + a3 * x3 + a4 * x4 + a5 * x5 - np.log(x)) * np.exp(x)
    else:
        a1 = 8.5733287401
        a2 = 18.0590169730
        a3 = 8.6347608925
        a4 = 0.2677737343
        b1 = 9.5733223454
        b2 = 25.6329561486
        b3 = 21.0996530827
        b4 = 3.9584969228
        eiexp = (x4 + a1 * x3 + a2 * x2 + a3 * x + a4) / (x4 + b1 * x3 + b2 * x2 + b3 * x + b4) / x

    return eiexp

def fnagex(pt, pcf):
    """
    This function calculates the electron-impact excitation rate coefficients 
     as a function of the  electron temperature.

    pt: electron temperature in eV
    pcf: parameter data array
        pcf[0]: type of the fit, to represent resonace and non-resonace region
        pcf[1]: excitation energy
        pcf[2]: lower limit of fitting of reduced energy x
        pcf[3]: upper limit of fitting of reduced energy x
        pcf[4]: statistical weight of initial state (2s+1)*(2l+1)
        pcf[5:14]: fitting parameters, depending upon the type of fit
    """
    const = 8.010e-8 / (np.sqrt(pt) * pcf[5])
    y = pcf[2] / pt

    if len(pcf) < 12:
        if pcf[1] == 1:
            ra = pcf[6] / y + pcf[8] + 0.5 * pcf[9] * (1.0 - y)
            rb = pcf[7] - pcf[8] * y + 0.5 * pcf[9] * y * y + pcf[10] / y
            rnor = const * np.exp(-y) * y * (ra + eiexp(y) * rb)
            rres = 0.0
        elif pcf[1] == 2:
            ra = pcf[6] * (1.0 - eiexp(y) * y)
            f = pcf[11]
            rb = pcf[7] * np.exp(-f) / (f + y) + pcf[8] * np.exp(-2.0 * f) / (2.0 * f + y) \
                + pcf[9] * np.exp(-3.0 * f) / (3.0 * f + y) + pcf[10] * np.exp(-4.0 * f) / (4.0 * f + y)
            rnor = const * np.exp(-y) * y * (ra + rb)
            rres = 0.0
        else:
            raise ValueError('Invalid integer for fit type in fnagex - must be 1 or 2')
    else:
        if pcf[1] == 1:
            x1 = pcf[13]
            y1 = y * x1
            ra = pcf[6] / y + pcf[8] / x1 + s * pcf[9] * (1.0 / (x1 * x1) - y / x1) + pcf[10] * np.log(x1) / y
            rb = pcf[7] - pcf[8] * y + pcf[9] * y * y + pcf[10] / y
            rnor = const * y * np.exp(-y1) * (ra + eiexp(y1) * rb)
            rres = const * np.exp(-y) * (pcf[11] * (1.0 + 1.0 / y) *
                                         (1.0 - np.exp((1.0 - x1) * y) * (x1 + 1.0 / y) / (1.0 + 1.0 / y))
                                         + pcf[12] * (1.0 - np.exp((1.0 - x1) * y)))
        elif pcf[1] == 2:
            x1 = pcf[14]
            y1 = y * x1
            ra = pcf[6] * (1.0 / x1 - eiexp(y1) * y)
            f = pcf[11]
            rb = pcf[6] * np.exp(-f * x1) / (f + y) \
                + pcf[7] * np.exp(-2.0 * f * x1) / (2.0 * f + y) \
                + pcf[8] * np.exp(-3.0 * f * x1) / (3.0 * f + y) \
                + pcf[9] * np.exp(-4.0 * f * x1) / (4.0 * f + y)
            rnor = const * y * np.exp(-y1) * (ra + rb)
            rres = const * np.exp(-y) * (pcf[12] * (1.0 + 1.0 / y) \
                * (1.0 - np.exp((1.0 - x1) * y) * (x1 + 1.0 / y) / (1.0 + 1.0 / y)) \
                + pcf[13] * (1.0 - np.exp((1.0 - x1) * y)))
        else:
            raise ValueError('Invalid integer for fit type in fnagex - must be 1 or 2')

    pfit = rnor + rres

    if pfit <= 0.0:
        raise ValueError('Error: Reaction rate is negative. Check data and temperature range.')

    return pfit

c
c###################################################################
c
      function eiexp(x)
      double precision x, x2, x3, x4, x5
      double precision a0, a1, a2, a3, a4, a5, b1, b2, b3, b4
c
c    eiexp=e1(x)*exp(x)
c    handbook of mathematical functions, page 231
c    by m. abramowitz and i. a. stegun
c    (routine copied from ippj-am-27, y. itikawa et al 1983)
      x2=x*x
      x3=x2*x
      x4=x3*x
      x5=x4*x
c
      if (x.le.1.0) then
        a0=-0.57721566
        a1=0.99999193
        a2=-0.24991055
        a3=0.05519968
        a4=-0.00976004
        a5=0.00107857
        eiexp=(a0 +a1*x +a2*x2 +a3*x3 +a4*x4 +a5*x5 - log(x))*exp(x)
c
      else
        a1=8.5733287401
        a2=18.0590169730
        a3=8.6347608925
        a4=.2677737343
        b1=9.5733223454
        b2=25.6329561486
        b3=21.0996530827
        b4=3.9584969228
        eiexp=(x4 +a1*x3+ a2*x2+ a3*x+ a4) /
     1       (x4 +b1*x3 +b2*x2 +b3*x+ b4)/x
      endif
c
      return
      end function eiexp
c###################################################################
c
      subroutine fnagex(pt, pcf, kncf, pfit, kermsg)
c
c     this is an iaea subroutine to calculate the electron impact
c     excitation rate coefficients as a function of the electron
c     temperature.
c
c     pt = electron temperature in ev
c
c     pcf is the coefficient data array, where
c
c     pcf(1)  = itype, index for type of fit, either 1 or 2
c
c     pcf(2)  = excitation energy , v (also referred to as delta e)
c
c     pcf(3)  = lower limit of fitting of reduced energy x
c
c     pcf(4)  = upper limit of fitting of reduced energy x
c
c     pcf(5)  = statistical weight of initial state (2s+1)*(2l+1)
c
c     pcf(6)  = parameter a
c
c     pcf(7)  = parameter b
c
c     pcf(8)  = parameter c
c
c     pcf(9)  = parameter d
c
c     pcf(10)  = parameter e
c
c if ftype = 1 this can be followed by three more parameters for
c the region which contains resonances represented by a linear term
c
c     pcf(11)  = parameter p
c
c     pcf(12)  = parameter q
c
c     pcf(13)  = parameter x1, the upper limit of the range over which
c                the collision strength is represented by a linear
c                approximation.
c
c if ftype = 2, the pcf array element 10 can be followed by
c
c     pcf(11)  = parameter f
c
c and possibly three more parameters follow for the region which
c contains resonances represented by a linear term
c
c     pcf(12)  = parameter p
c
c     pcf(13)  = parameter q
c
c     pcf(14)  = parameter x1, the upper limit of the range over which
c                the collision strength is represented by a linear
c                approximation.
c
c     kermsg = blank if no errors
c
c     pfit = rate coefficient in cm[3]/s
c
c     written by j. j. smith , iaea atomic and molecular data unit
c     (taken from report ippj-am-27, y. itikawa et al, nagoya,
c      institute of plasma physics, nsgoya univ., (1983))
c
c------------------------------------------------------------------------
c
      double precision pt, pcf, pfit, eiexp
      dimension pcf(14)
      character*(*) kermsg
      data s/0.5/
c
c	DH, Dec, 16, 2003 ---- kermsg = ' ' to avoid loop on error message ----
	kermsg = ' '
      itype = pcf(1)
      const = 8.010e-8 / (dsqrt(pt)*pcf(5))
      y = pcf(2) / pt
c
      if (kncf .lt. 12) then
c
c---    rate coefficient without resonances
c
        if (itype .eq. 1) then
c
c---    rate coefficient represented by power-log type of fit
c
          ra = pcf(6)/y + pcf(8) + s * pcf(9) * (1.0 - y)
          rb = pcf(7) - pcf(8)*y + s * pcf(9) * y * y + pcf(10)/y
          rnor = const * exp(-y) * y *  (ra + eiexp(y) * rb)
          rres = 0.0
c
        else if (itype .eq. 2) then
c
c---    rate coefficient represented by exponential type of fit
c
          ra = pcf(6) * (1.0 - eiexp(y)*y)
          f = pcf(11)
          rb = pcf(7) * exp(-f) / (f+y)
     1         + pcf(8) * exp (-2.*f) / (2.*f+y)
     2         + pcf(9) * exp (-3.*f) / (3.*f+y)
     3         + pcf(10) * exp (-4.*f) / (4.*f+y)
          rnor = const * exp(-y) * y *  (ra + rb)
          rres = 0.0
c
        else
            kermsg =
     1       ' invalid integer for fit type in fnagex - must be 1 or 2'
          return
        endif
c
      else
c
c---    rate coefficient with resonances
c
        if (itype .eq. 1) then
c
c---    rate coefficient represented by power-log type of fit
c
          x1 = pcf(13)
          y1 = y * x1
          ra = pcf(6)/y + pcf(8)/x1
     1         + s * pcf(9) * (1./(x1 * x1)  - y/x1)
     2         + pcf(10) * log(x1) / y
          rb = pcf(7) - pcf(8)*y + s * pcf(9) * y * y + pcf(10)/y
          rnor = const * y * exp(-y1) * (ra + eiexp(y1) * rb)
          rres = const * exp (-y) * (pcf(11) * (1. + 1./y)
     1         * (1. - exp ((1. - x1) * y) * (x1 + 1./y)
     2         / (1. + 1./y))  + pcf(12) * (1. - exp ((1. - x1) * y )))
c
        else if (itype .eq. 2) then
c
c---    rate coefficient represented by exponential type of fit
c
          x1 = pcf(14)
          y1 = y * x1
          ra = pcf(6) * (1. / x1 - eiexp(y1) * y)
          f = pcf(11)
          rb = pcf(7) * exp(-f*x1) / (f+y)
     1        + pcf(8) * exp (-2.*f*x1) / (2.*f+y)
     2        + pcf(9) * exp (-3.*f*x1) / (3.*f+y)
     3        + pcf(10) * exp (-4.*f*x1) / (4.*f+y)
          rnor = const * y * exp(-y1) * (ra + rb)
          rres = const * exp (-y) * (pcf(12) * (1. + 1./y)
     1         * (1. - exp ((1. - x1) * y) * (x1 + 1./y)
     2         / (1. + 1./y))  + pcf(13) * (1. - exp ((1. - x1) * y )))
        else
            kermsg =
     1       ' invalid integer for fit type in fnagex - must be 1 or 2'
          return
        endif
c
      endif
c
      pfit = rnor + rres
c
      if (pfit .le. 0.0d0) then
            kermsg =
     1  ' error reaction rate is negative check data and temp. range '
          return
        endif
c
      end