First implementation of the photomeson interactions in photo_meson.py

agnpy/photo_meson/kernels.py

@@ -2,6 +2,7 @@
 from scipy.interpolate import CubicSpline
 import numpy as np
 from pathlib import Path
+import astropy.units as u
 
 data_dir = Path(__file__).parent.parent
 secondaries = [
@@ -40,11 +41,14 @@ def interpolate_phi_parameter(particle, parameter):
     )
 
     if parameter == "s":
-        func = CubicSpline(eta_eta0, s)
+        #func = CubicSpline(eta_eta0, s)
+        func = lambda x: np.interp(x, eta_eta0, s)
     elif parameter == "delta":
-        func = CubicSpline(eta_eta0, delta)
+        #func = CubicSpline(eta_eta0, delta)
+        func = lambda x: np.interp(x, eta_eta0, delta)
     elif parameter == "B":
-        func = CubicSpline(eta_eta0, B)
+        #func = CubicSpline(eta_eta0, B)
+        func = lambda x: np.interp(x, eta_eta0, B)
     else:
         raise ValueError(
             f"{parameter} not available among the parameters to be interpolated"
@@ -182,7 +186,7 @@ def __init__(self, particle):
             ]:
                 self.x_minus = x_minus_leptons_1
                 self.x_plus = x_plus_gamma
-            elif self.particle in ["electrons", "electron_antineutrino"]:
+            elif self.particle in ["electron", "electron_antineutrino"]:
                 self.x_minus = x_minus_leptons_2
                 self.x_plus = x_plus_leptons_2
             elif self.particle == "muon_neutrino":
@@ -201,7 +205,7 @@ def __call__(self, eta, x):
         # evaluate the interpolated parameters
         s = self.s(eta / eta_0)
         delta = self.delta(eta / eta_0)
-        B = self.B(eta / eta_0)
+        B = self.B(eta / eta_0) * u.Unit("cm3 s-1")
 
         x_minus = self.x_minus(eta)
         x_plus = self.x_plus(eta)

agnpy/photo_meson/photo_meson.py

@@ -1,12 +1,31 @@
 # class for photomeson production
 import numpy as np
+import astropy.units as u
+from astropy.constants import c, h, m_e
 from .kernels import PhiKernel
 from ..utils.math import axes_reshaper
 from ..utils.conversion import mpc2
 
 
 class PhotoMesonProduction:
-    def __init__(self, blob, target, integrator=np.trapz):
+    """Class for computation of the energetic spectra of secondaries of photomeson interactions.
+
+    Parameters
+    ----------
+    blob : :class:`~agnpy.emission_region.Blob`
+        emitting region with a proton distribution
+    target : ...TBD...
+        ...TBD...
+    integrator : func
+        function to be used for integration (default = `np.trapz`)
+    """
+
+    def __init__(
+        self,
+        blob,
+        target,
+        integrator=np.trapz
+    ):
         self.blob = blob
         # check that this blob has a proton distribution
         if self.blob._n_p is None:
@@ -15,38 +34,96 @@ def __init__(self, blob, target, integrator=np.trapz):
             )
         self.target = target
         self.integrator = integrator
-        self.phi_kernel_gamma = PhiKernel("gamma")
-        self.phi_kernel_electron = PhiKernel("electron")
-        self.phi_kernel_positron = PhiKernel("positron")
-        self.phi_kernel_electron_neutrino = PhiKernel("electron_neutrino")
-        self.phi_kernel_electron_antineutrino = PhiKernel("electron_antineutrino")
-        self.phi_kernel_muon_neutrino = PhiKernel("muon_neutrino")
-        self.phi_kernel_moun_antineutrino = PhiKernel("muon_antineutrino")
+        return
+
 
-    def H(self, eta, E, phi_kernel):
-        """Compute the H function in Eq. (70) [KelnerAharonian2008]_.
+    def H(
+        self,
+        eta,
+        E,
+        phi_kernel,
+        integrator = np.trapz
+    ):
+        """ Compute the H function in Eq. (70) [KelnerAharonian2008]_.
 
         Parameters
         ----------
+        E : :class:`~astropy.units.Quantity`
+            energy of the secondary particles
         eta : float
             kinematic variable (:math:`eta = 4 \epsilon \gamma_{\rm p}`)
-        E : float
-            energy of the secondary
         phi_kernel : `~agnpy.photo_meson.PhiKernel`
             kernel to be used for the integration (depends on the particle)
+        integrator : func
+            function to be used for integration (default = `np.trapz`)
         """
-        # we want this output to be a two-dimensional function, so let us reshape it
-        _gamma_p, _eta, _E = axes_reshaper(self.blob.gamma_p, eta, E)
-        _x = (_E / (_gamma_p * mpc2)).to_value("")
-        # obtain the energy of the target from eta
-        _E_ph = _eta * mpc2 / (4 * _gamma_p)
+        # Integral on E_p to be made from E to infinity
+        _eta, _E = axes_reshaper(eta, E)   # shape (len(eta), 1), (1, len(E))
+        _E_p = np.logspace(
+            log10(_E.to_value("eV")),
+            log10(_E.to_value("eV")) + 8,
+            200
+        ) * u.Unit("eV")                   # shape (200, 1, len(E))
+        _gamma_p = _E_p / mpc2
+        _epsilon = _eta * mpc2**2 / (4*_E_p)
+        _nu = _epsilon / h
+        _x = _E / _E_p
         _H_integrand = (
-            mpc2
-            / 4
-            * self.blob.n_p(_gamma_p)
-            * self.target.n_ph(_E_ph)
-            * phi_kernel(_eta, _x)
-        )
-        import IPython
+            mpc2**2 / 4                           # erg2
+            * (self.blob.n_p(_gamma_p) / _E_p**3) # cm-3 erg-3
+            * self.target(_nu)                    # cm-3 erg-1
+            * phi_kernel(_eta, _x)                # cm3 s-1
+        ).to("erg-2 cm-3 s-1")
 
-        IPython.embed()
+        _H = integrator(
+            _H_integrand,
+            _E_p,
+            axis = 0,
+        ).to("erg-1 cm-3 s-1")
+        return _H
+
+
+    def evaluate_spectrum (
+        self,
+        E,
+        particle,
+        integrator = np.trapz
+    ):
+        """ Evaluate the spectrum of secondaries in the emission region reference frame
+        as in Eq. (69) [KelnerAharonian2008]_.
+
+        Parameters
+        ----------
+        E : float
+            energy of the secondary particles
+        particle: str
+            name of the secondary particle, to be chosen among
+            "gamma", "electron", "positron", "electron_neutrino",
+            "electron_antineutrino", "muon_neutrino", "muon_antineutrino"
+        integrator : func
+            function to be used for integration (default = `np.trapz`)
+        """
+        if particle not in secondaries:
+            raise AttributeError(
+                f"There is no secondary particle from photomeson interactions named {particle}."
+            )
+
+        phi_kernel = PhiKernel(particle)
+        # Integral on eta to be done from eta_0 to infinity
+        eta = np.logspace(
+            log10(eta_0),
+            log10(eta_0) + 5,
+            100,
+        )
+        _H = self.H(
+            eta,
+            E,
+            phi_kernel,
+            integrator = integrator,
+        )
+        dN_dEdVdt = integrator(
+            _H,
+            eta,
+            axis = 0
+        ).to("erg-1 cm-3 s-1")
+        return dN_dEdVdt

agnpy/utils/math.py

@@ -52,6 +52,12 @@ def log(x):
     return np.log(values)
 
 
+def log10(x):
+    """clipped log to avoid RuntimeWarning: divide by zero encountered in log"""
+    values = np.clip(x, ftiny, fmax)
+    return np.log10(values)
+
+
 def trapz_loglog(y, x, axis=0):
     """
     Integrate a function approximating its discrete intervals as power-laws