Fit Function: NEEXH1

Fit Function \[\begin{align*} \sigma (pe) &= \begin{cases} 10^{-16} \left[ \text{pcf}(3) + \text{pcf}(4) \cdot (pe - \text{pcf}(1)) \right] & \text{if } pe \leq \text{pcf}(2) \\ \text{pcf}(6) & \text{if } \text{pcf}(2) < pe \leq \text{pcf}(5) \\ \frac{5.984 \times 10^{-16} \text{cstr}}{pe} & \text{otherwise} \\ \end{cases} \\ \end{align*}\]
Comments Python code requires NumPy imported as `np` as well as `warnings`. Please refer to the source code for the information on `cstr`.

Fortran

Arguments
namedescriptionunitstype(s)
pe energy eV real, dimension(:)
pcf coefficient data array real, dimension(12)
kncf number of coefficients in the data array integer
pxs cross section cm2 real, dimension(:)
kermsg error message character
Return values
namedescriptionunitstype(s)
pxs cross section cm2 real, dimension(:)
Code
c
c##############################################################################
c
      subroutine neexh1(pe, pcf, kncf, pxs, kermsg)
c
c     this is a subroutine to calculate cross sections (cm[2])
c     versus energy (ev) for electron impact excitation.
c
c     pe = collision energy in ev
c
c     pcf(1) = threshold energy (ev), eth.
c     pcf(2) = upper value of energy range over which the
c              cross section, sig, is approximated in a linear form
c                     sig = a + b*(e-eth)
c              where a,b and c are parameters and e is the incident
c              energy in ev
c              (lower end of the range is the threshold energy)
c
c     pcf(3) = parameter a
c     pcf(4) = parameter b
c     pcf(5) = upper value of energy range over which the cross section
c              is approximated by a constant value. (lower limit for
c              this range is pcf(2).
c     pcf(6) = constant value for second region.
c     pcf(7-12) = parameters for fit to the cross section in the third
c              region.
c
c     kermsg = blank if no errors
c
c     pxs = cross section in cm[2]
c
c------------------------------------------------------------------------
c
      double precision pe, pcf, pxs
      double precision cstr, ryd, ksq, x, alog
c
      dimension pcf(12)
      character*(*) kermsg
c
      data ryd/1.36d+01/
c
      eth = pcf(1)
      if(pe .ge. eth) then
        kermsg = ' '
      else
        pxs = 0.0d0
        return
      endif
c
      if (pe .le. pcf(2)) then
        pxs =  1.0e-16 * (pcf(3) + (pcf(4) * (pe - eth)))
c
      else if (pe .gt. pcf(2) .and. pe .le. pcf(5)) then
         pxs =  pcf(6)
c
      else
        ksq= pe/ryd
        etrans = eth/ryd
        x=ksq/etrans
        xsq = x*x
        alog = 0.0d0
        if (pcf(12) .ne. zero) alog = pcf(12) * dlog (x)
        cstr = alog + pcf(7) + pcf(8)/x + pcf(9)/xsq + pcf(10)/(x*xsq)
     &        + pcf(11)/(xsq*xsq)
        pxs =   5.984d-16 * cstr / pe
c
      endif
c
      return
      end

Python

Arguments
namedescriptionunitstype(s)
pe energy eV float, np.ndarray
pcf coefficient data array float, np.ndarray
Return values
namedescriptionunitstype(s)
pxs cross section cm2 float, np.ndarray
Code
def neexh1(pe, pcf):
    """
    This is a function to calculate cross sections (cm^2) versus energy (eV)
    for electron impact excitation.

    pe: collision energy (eV)
    pcf: fit coefficient array
        pcf[0]: threshold energy (eV), eth
        pcf[1]: upper value of energy range for linear approximation
        pcf[2], pcf[3]: fit parameters for the linear approximation 
        pcf[4]: upper value of energy range for constant approximation
        pcf[5]: constant value for the second region
        pcf[6:12]: parameters for fit to the cross section in the third region
    """

    ryd = 13.6
    eth = pcf[0]

    if pe < eth:
        warnings.warn('Energy below threshold. Cross section set to 0.')
        pxs = 0.0
        return pxs

    if pe <= pcf[1]:
        pxs = 1.0e-16 * (pcf[2] + pcf[3] * (pe - eth))

    elif pe > pcf[1] and pe <= pcf[4]:
        pxs = pcf[5]

    else:
        ksq = pe/ryd
        etrans = eth/ryd
        x = ksq/etrans
        xsq = x**2
        alog = 0.0
        if pcf[11] != 0.0:
            alog = pcf[11] * np.log(x)
        cstr = alog + pcf[6] + pcf[7]/x + pcf[8]/xsq + pcf[9]/(x * xsq) 
                   + pcf[10]/(xsq * xsq)
        pxs = 5.984e-16 * cstr/pe

    return pxs