Fit Function: PHACX

Fit Function \[\begin{align*} \sigma(\text{pe}) &= \frac{\text{pcf}(3) \cdot \text{pcf}(7) \cdot \ln\left(\frac{\text{pcf}(4) \sqrt{\text{pcf}(7)}}{{\text{pek}}}\right)}{1 + \frac{\text{pcf}(5) \cdot \text{pek}^2}{\text{pcf}(7)} + \text{pcf}(6) \left(\frac{\text{pek}}{\sqrt{\text{pcf}(7)}}\right)^{4.5}}; \quad \text{pek} = 10^{-3} \; \text{pe} \end{align*}\]
Comments pcf(7) = q, the charge state of the projectile ion. This must be input into the fit coefficient array, as it is not included in the fit coefficients. Python code requires NumPy imported as `np`.

Fortran

Arguments
namedescriptionunitstype(s)
pe projectile energy eV amu-1 real, dimension(:)
pcf coefficient data array real, dimension(7)
kncf number of coefficients in the data array integer
pfit cross section cm2 real, dimension(:)
kermsg error message character
Return values
namedescriptionunitstype(s)
pfit cross section cm2 real, dimension(:)
Code
c
c###################################################################
c
      subroutine alphacx(pe, pcf, kncf, pfit, kermsg)
c
c     this is a subroutine to calculate cross section for
c     projectile energy (ev/amu).
c
c     the approximate range of valiity of the fit is given by the
c     parameters pcf(1) and pcf(2) which are the scaled energies
c     pe/dsqrt(q) (ev/u), where q is the charge state of the ion.
c
c     pe = projectile energy (ev/amu)
c
c     kermsg = blank if no errors
c
c     pfit = cross section in cm[2]
c
c     pcf(1) = lower limit of the range of validity of the scaled
c              projectile energy in ev/(u*sqrt(q)).
c     pcf(2) = upper limit of the range of validity of the scaled
c              projectile energy in ev/(u*sqrt(q)).
c     pcf(3-6) = a, b, c, d fitting parameters
c     pcf(7) = q, the charge state of the projectile ion. this must be
c              input into the coefficient array. this is not included
c              in the data entries.
c
c     written by j. j. smith , iaea atomic and molecular data unit
c     corrected for invalid definition of energy pek (30/01/90)
c
c------------------------------------------------------------------------
c
      double precision pe, pcf, pfit
      double precision pek, q, q2, a, b, c, d, f1, f2, elow, ehigh
      dimension pcf(7)
      character*(*) kermsg
      a = pcf(3)
      b = pcf(4)
      c = pcf(5)
      d = pcf(6)
      q = pcf(7)
      q2 = dsqrt(q)
c
c---  determine the projectile energy - pek in kev/amu
c
      pek = pe*1.0d-03
c
c---  check that the energy is within the valid energy ranges
c
      elow = pcf(1)*1.0d-03
      ehigh = pcf(2)*1.0d-03
c
      kermsg = ' '
c
      if(q .lt. 5 .or. q .gt. 26) then
        kermsg =
     &   'the charge state (q) must be in the range 5 <= q <= 26'
        return
      endif
c
      f1 = a * q * dlog ( b * q2 / pek )
      f2 = ( c * pek * pek / q )  +  d * ( ( pek / q2 ) ** 4.5d0 )
      pfit = f1 / ( 1.0d0 + f2 )
c
      return
c
      end

Python

Arguments
namedescriptionunitstype(s)
pe projectile energy eV amu-1 float, np.ndarray
pcf coefficient data array float, np.ndarray
Return values
namedescriptionunitstype(s)
pfit cross section cm2 float, np.ndarray
Code
def alphacx(pe, pcf):
    """
    This function calculates cross-sections in cm[2] as a function
    of projectile energy in ev/amu.

    pe: projectile energy (ev/amu)
    pcf: fit coefficient array
        pcf[0]: lower limit of the range of validity of the scaled projectile energy 
                     in ev/(u*sqrt(q))
        pcf[1]: upper limit of the range of validity of the scaled projectile energy
                    in ev/(u*sqrt(q))
        pcf[2] - pcf[6]: fit parameters 
        pcf[6]: charge state of the projectile ion (q)
        pcf[7]: unused fit parameter
        pcf[8]: unused fit parameter
    """
    q = pcf[6]
    q2 = np.sqrt(q)

    # Determine the projectile energy - pek in kev/amu
    pek = pe * 1.0e-03

    # Check that the energy is within the valid energy ranges
    elow = pcf[0] * 1.0e-03
    ehigh = pcf[1] * 1.0e-03

    if pek < elow or pek > ehigh:
        raise ValueError('Energy outside the validity range')

    if q < 5 or q > 26:
        raise ValueError('The charge state (q) must be in the range 5 <= q <= 26')

    f1 = pcf[2] * q * np.log(pcf[3] * q2 / pek)
    f2 = pcf[4] * pek * pek / q + pcf[5] * (pek / q2) ** 4.5
    pfit = f1 / (1.0 + f2)

    return pfit