Add a function for integrating SED flux, for different types of radiative processes

agnpy/compton/external_compton.py

@@ -9,6 +9,7 @@
 )
 from ..utils.conversion import nu_to_epsilon_prime, to_R_g_units
 from ..utils.geometry import x_re_shell, mu_star_shell, x_re_ring
+from ..utils.sedintegrable import SedFluxIntegrable
 from ..targets import (
     CMB,
     PointSourceBehindJet,
@@ -21,7 +22,7 @@
 __all__ = ["ExternalCompton"]
 
 
-class ExternalCompton:
+class ExternalCompton(SedFluxIntegrable):
     """class for External Compton radiation computation
 
     Parameters

agnpy/compton/synchrotron_self_compton.py

@@ -10,12 +10,12 @@
     nu_to_integrate,
 )
 from ..utils.conversion import nu_to_epsilon_prime
-
+from ..utils.sedintegrable import SedFluxIntegrable
 
 __all__ = ["SynchrotronSelfCompton"]
 
 
-class SynchrotronSelfCompton:
+class SynchrotronSelfCompton(SedFluxIntegrable):
     """class for Synchrotron Self Compton radiation computation
 
     Parameters

agnpy/synchrotron/proton_synchrotron.py

@@ -4,14 +4,16 @@
 from astropy.constants import e, c, m_p
 from ..utils.math import axes_reshaper, gamma_e_to_integrate
 from ..utils.conversion import nu_to_epsilon_prime, B_to_cgs, lambda_c_p
+from ..utils.sedintegrable import SedFluxIntegrable
 from .synchrotron import single_particle_synch_power, tau_to_attenuation
 
 __all__ = ["ProtonSynchrotron"]
 
 e = e.gauss
 B_cr = 4.414e13 * u.G  # critical magnetic field
 
-class ProtonSynchrotron:
+
+class ProtonSynchrotron(SedFluxIntegrable):
     """Class for synchrotron radiation computation
 
     Parameters

agnpy/synchrotron/synchrotron.py

@@ -4,6 +4,7 @@
 from astropy.constants import e, h, c, m_e, sigma_T
 from ..utils.math import axes_reshaper, gamma_e_to_integrate
 from ..utils.conversion import nu_to_epsilon_prime, B_to_cgs, lambda_c_e
+from ..utils.sedintegrable import SedFluxIntegrable
 
 
 __all__ = ["R", "nu_synch_peak", "Synchrotron"]
@@ -67,7 +68,7 @@ def tau_to_attenuation(tau):
     return np.where(tau < 1e-3, 1, 3 * u / tau)
 
 
-class Synchrotron:
+class Synchrotron(SedFluxIntegrable):
     """Class for synchrotron radiation computation
 
     Parameters

agnpy/tests/test_sedintegrable.py

@@ -0,0 +1,82 @@
+import numpy as np
+import math
+import astropy.units as u
+from astropy.constants import m_e, h
+from astropy.coordinates import Distance
+from agnpy.spectra import PowerLaw
+from agnpy.emission_regions import Blob
+from agnpy.synchrotron import Synchrotron
+from ..utils.sedintegrable import SedFluxIntegrable
+
+
+def _estimate_powerlaw_integral(x_values, y_values):
+    slope, intercept = np.polyfit(np.log10(x_values.value), np.log10(y_values.value), 1)
+    power = slope + 1
+    estimated_integral = (10 ** intercept) * (
+            x_values[-1].value ** power / power - x_values[0].value ** power / power)
+    return estimated_integral
+
+
+def _sample_synchrotron_model():
+    # Based on the synchrotron example from agnpy documentation
+    R_b = 1e16 * u.cm
+    n_e = PowerLaw.from_total_energy(1e48 * u.erg, 4 / 3 * np.pi * R_b ** 3, p=2.8, gamma_min=1e2,
+                                     gamma_max=1e7, mass=m_e)
+    blob = Blob(R_b, z=Distance(1e27, unit=u.cm).z, delta_D=10, Gamma=10, B=1 * u.G, n_e=n_e)
+    synch = Synchrotron(blob)
+    return synch
+
+
+class TestSedIntegrable:
+
+    def test_flat_integral(self):
+        """Integrate over flat SED (Fnu equal to 1.0 for all frequencies)"""
+        class FlatSedGenerator(SedFluxIntegrable):
+            def sed_flux(self, nu):
+                return np.full(fill_value=1.0, shape=nu.shape, dtype=np.float64) * (u.erg / (u.cm ** 2 * u.s * u.Hz)) * nu
+
+        nu_min = 10 * u.Hz
+        nu_max = 10 ** 20 * u.Hz
+        flux = FlatSedGenerator().integrate_flux(nu_min, nu_max)
+        assert flux == 1e+20 * u.erg / (u.cm ** 2 * u.s)
+
+
+    def test_synchrotron_energy_flux_integral(self):
+        """Integrate over sample synchrotron flux and compare with the value estimated using the power law integral"""
+        synch = _sample_synchrotron_model()
+
+        nu_min_log = 13
+        nu_max_log = 20
+        nu_min = 10 ** nu_min_log * u.Hz
+        nu_max = 10 ** nu_max_log * u.Hz
+
+        nu = np.logspace(nu_min_log, nu_max_log) * u.Hz
+        Fnu = synch.sed_flux(nu) / nu
+        estimated_integral = _estimate_powerlaw_integral(nu, Fnu)
+
+        actual_integral = synch.integrate_flux(nu_min, nu_max).value
+
+        assert math.isclose(actual_integral, estimated_integral, rel_tol=0.05)
+
+    def test_synchrotron_photon_flux_integral(self):
+        """Integrate over sample synchrotron photon flux and compare with the value estimated using the power law integral"""
+        synch = _sample_synchrotron_model()
+
+        nu_min_log = 13
+        nu_max_log = 20
+        nu_min = 10 ** nu_min_log * u.Hz
+        nu_max = 10 ** nu_max_log * u.Hz
+
+        nu = np.logspace(nu_min_log, nu_max_log) * u.Hz
+        Fnu = synch.sed_flux(nu) / nu
+        # convert the energy flux to photon flux in photons / s / cm2 / eV
+        photons_per_Hz = Fnu / nu.to("erg", equivalencies=u.spectral())  # Fnu is in ergs so must use the same unit
+        h_in_eV = h.to("eV/Hz")
+        energies_eV = nu * h_in_eV
+        photons_per_eV = photons_per_Hz / h_in_eV
+
+        estimated_integral = _estimate_powerlaw_integral(energies_eV, photons_per_eV)
+
+        actual_integral = synch.integrate_photon_flux(nu_min, nu_max).value
+
+        assert math.isclose(actual_integral, estimated_integral, rel_tol=0.05)

agnpy/utils/sedintegrable.py

@@ -0,0 +1,40 @@
+import numpy as np
+import astropy.units as u
+
+
+class SedFluxIntegrable:
+    def integrate_flux(self, nu_min, nu_max, nu_points=50):
+        r""" Evaluates the SED flux integral over the span of provided frequencies
+
+        Parameters
+        ----------
+        nu_min : start frequency (in Hz or equivalent)
+        nu_max : end frequency (in Hz or equivalent)
+        nu_points: number of points (between nu_min and nu_max) on the log scale
+        """
+        nu = self._nu_logspace(nu_min, nu_max, nu_points)
+        sed_nufnu = self.sed_flux(nu)  # erg / s / cm2 / log(nu)
+        sed_fnu = sed_nufnu / nu       # erg / s / cm2 / Hz
+        return np.trapz(sed_fnu, nu)   # erg / s / cm2
+
+    def integrate_photon_flux(self, nu_min, nu_max, nu_points=50):
+        r""" Evaluates the photon flux integral from SED over the span of provided frequencies
+
+        Parameters
+        ----------
+        nu_min : start frequency (in Hz or equivalent)
+        nu_max : end frequency (in Hz or equivalent)
+        nu_points: number of points (between nu_min and nu_max) on the log scale
+        """
+        nu = self._nu_logspace(nu_min, nu_max, nu_points)
+        photon_energy = nu.to("erg", equivalencies=u.spectral())
+        sed_nufnu = self.sed_flux(nu)   # erg / s / cm2 / log(nu)
+        sed_fnu = sed_nufnu / nu        # erg / s / cm2 / Hz
+        n = sed_fnu / photon_energy     # photons / s / cm2 / Hz
+        return np.trapz(n, nu)          # photons / s / cm2
+
+    def _nu_logspace(self, nu_min, nu_max, nu_points=50):
+        nu_min_hz = nu_min.to(u.Hz, equivalencies=u.spectral())
+        nu_max_hz = nu_max.to(u.Hz, equivalencies=u.spectral())
+        nu = np.logspace(start=np.log10(nu_min_hz.value), stop=np.log10(nu_max_hz.value), num=nu_points) * u.Hz
+        return nu