Fit Function: JRREC2

Fit Function \[\sigma (pt) = \text{pcf}(2) \times 10^{-14} \left(\frac{\text{pcf}(1)}{pt}\right)^{1.5} \exp\left(\frac{\text{pcf}(1)}{pt}\right) \text{exint}\left(\frac{\text{pcf}(1)}{pt}\right)\]
Comments Python code requires NumPy imported as `np`.

Fortran

Arguments
namedescriptionunitstype(s)
pt electron temperature eV real, dimension(:)
pcf coefficient data array real, dimension(9)
kncf number of coefficients in the data array integer
pfit rate coefficient cm3 s-1 real, dimension(:)
kermsg error message character
Return values
namedescriptionunitstype(s)
pfit rate coefficient cm3 s-1 real, dimension(:)
Code 
c
c###################################################################
c
      function exint(x)
c
c    exint = exponential integral
c    handbook of mathematical functions, page 231
c    by m. abramowitz and i. a. stegun
c
c     written by j. j. smith , iaea atomic and molecular data unit
c
c------------------------------------------------------------------------
c
      x2=x*x
      x3=x2*x
      x4=x3*x
      x5=x4*x
      if (x.gt.1.0) go to 10
      a0 = -0.57721566
      a1 =  0.99999193
      a2 = -0.24991055
      a3 =  0.05519968
      a4 = -0.00976004
      a5 =  0.00107857
      exint=a0+a1*x+a2*x2+a3*x3+a4*x4+a5*x5 - log(x)
      return
 10   a1 =  8.5733287401
      a2 = 18.0590169730
      a3 =  8.6347608925
      a4 =   .2677737343
      b1=   9.5733223454
      b2=  25.6329561486
      b3=  21.0996530827
      b4=   3.9584969228
      exint=(x4+a1*x3+a2*x2+a3*x+a4)/(x4+b1*x3+b2*x2+b3*x+b4)/(x*exp(x))
      return
      end

c###################################################################
c
      subroutine jrrec2(pt, pcf, kncf, pfit, kermsg)
c
c     this is a subroutine to calculate a reaction rate coefficient
c     in cm[3]/s as function of electron tempreature in ev for
c     radiative recombination for specific transitions to excited
c     states of neutral helium.
c
c     pcf is the coefficient data array, where
c
c     pcf(1)  =   threshold energy (ev)
c
c     pcf(2)  =   coefficient a (nl)
c
c     pt = tepreature (ev)
c
c     kermsg = blank if no errors
c
c     pfit = reaction rate coefficient in cm[3]/s
c
c     written by j. j. smith , iaea atomic and molecular data unit
c
c------------------------------------------------------------------------
c
      double precision pt, pcf, pfit
      dimension pcf(2)
      character*(*) kermsg
c
      kermsg = ' '
      eth=pcf(1)
      a=pcf(2)
c
      beta = eth/pt
c
      pfit = a * 1.0e-14 * (beta**1.5) * exp(beta) * exint(beta)
      return
c
      end

Python

Arguments
namedescriptionunitstype(s)
pt electron temperature eV float, np.ndarray
pcf coefficient data array float, np.ndarray
Return values
namedescriptionunitstype(s)
pfit rate coefficient cm3 s-1 float, np.ndarray
Code 
from scipy.special import expi

def jrrec2(pt, pcf):
    """
    This function calculates a reaction rate coefficient in cm3/s as a function of
    electron temperature in eV for radiative recombination for specific transitions to
    excited states of neutral helium.

        pt: Electron temperature (eV)
        pcf: Coefficient data array
            pcf[0]: Threshold energy (eV)
            pcf[1]: Fitting coefficient
    """

    eth = pcf[0]
    beta = eth / pt

    pfit = pcf[1] * 1.0e-14 * (beta**1.5) * np.exp(beta) * expi(beta)

    return pfit