docs: def _nuclear_factor(self, Tp):
          """
          Compute nuclear enhancement factor
          """
          sigmaRpp = 10 * np.pi * 1e-27
          sigmainel = self._sigma_inel(Tp)
          sigmainel0 = self._sigma_inel(1e3)  # at 1e3 GeV
          f = sigmainel / sigmainel0
          f2 = np.where(f > 1, f, 1.0)
          G = 1.0 + np.log(f2)
          # epsilon factors computed from Eqs 21 to 23 with local ISM abundances
          epsC = 1.37
          eps1 = 0.29
          eps2 = 0.1
  
          epstotal = np.where(
              Tp > self._Tth,
              epsC + (eps1 + eps2) * sigmaRpp * G / sigmainel,
              0.0,
          )
  
          if np.any(Tp < 1.0):
              # nuclear enhancement factor diverges towards Tp = Tth, fix Tp<1 to
              # eps(1.0) = 1.91
              loE = np.where((Tp > self._Tth) * (Tp < 1.0))
              epstotal[loE] = 1.9141
  
          return epstotal