docs: def _sigma_1(gam, eps):
          """
          gam and eps in units of m_e c^2
          Eq. A2 of Baring et al. (1999)
          Return in units of cm2 / mec2
          """
          s1 = 4 * r0 ** 2 * alpha / eps / mec2_unit
          s2 = 1 + (1.0 / 3.0 - eps / gam) * (1 - eps / gam)
          s3 = np.log(2 * gam * (gam - eps) / eps) - 1.0 / 2.0
          s3[np.where(gam < eps)] = 0.0
          return s1 * s2 * s3