Merge pull request #153 from grzegorzbor/master

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

.github/workflows/test.yml

@@ -12,7 +12,7 @@ jobs:
     runs-on: ubuntu-latest
     strategy:
       matrix:
-        python-version: [3.8, 3.9, 3.10.10]
+        python-version: [3.9, 3.10.10]
 
     steps:
     - uses: actions/checkout@v3

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 ..radiative_process import RadiativeProcess
 from ..targets import (
     CMB,
     PointSourceBehindJet,
@@ -21,7 +22,7 @@
 __all__ = ["ExternalCompton"]
 
 
-class ExternalCompton:
+class ExternalCompton(RadiativeProcess):
     """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 ..radiative_process import RadiativeProcess
 
 __all__ = ["SynchrotronSelfCompton"]
 
 
-class SynchrotronSelfCompton:
+class SynchrotronSelfCompton(RadiativeProcess):
     """class for Synchrotron Self Compton radiation computation
 
     Parameters

agnpy/radiative_process/__init__.py

@@ -0,0 +1 @@
+from .core import *

agnpy/radiative_process/core.py

@@ -0,0 +1,59 @@
+# core class describing all radiative processes
+import numpy as np
+import astropy.units as u
+
+
+def _nu_logspace(nu_min, nu_max, nu_points=50):
+    """A thin wrapper around numpy.logspace() function, capable of
+    spectral-equivalency conversion of input parameters."""
+    nu_min = nu_min.to(u.Hz, equivalencies=u.spectral())
+    nu_max = nu_max.to(u.Hz, equivalencies=u.spectral())
+    nu = np.logspace(np.log10(nu_min.value), np.log10(nu_max.value), nu_points) * u.Hz
+    return nu
+
+
+class RadiativeProcess:
+    """Base class for all the radiative processes.
+    It contains functionalities to handle all the basic flux calculations, e.g.
+    conversion to differential energy flux, integration, etc.
+    """
+
+    def diff_energy_flux(self, nu):
+        """Similar to sed_flux(), but returns the differential energy flux [eV-1 cm-2 s-1]."""
+        energy = nu.to("eV", equivalencies=u.spectral())
+        nuF_nu = self.sed_flux(nu)
+        dphidE = nuF_nu / energy**2
+        return dphidE.to("eV-1 cm-2 s-1")
+
+    def energy_flux_integral(self, nu_min, nu_max, nu_points=50):
+        """Evaluates the integral energy flux [erg cm-2 s-1] over a range of frequencies.
+
+        Parameters
+        ----------
+        nu_min : `~astropy.units.Quantity`
+            start frequency (in Hz or equivalent)
+        nu_max : `~astropy.units.Quantity`
+            end frequency (in Hz or equivalent)
+        nu_points: int
+            number of points (between nu_min and nu_max) in log scale
+        """
+        nu = _nu_logspace(nu_min, nu_max, nu_points)
+        Fnu = self.sed_flux(nu) / nu
+        return np.trapz(Fnu, nu).to("erg cm-2 s-1")
+
+    def flux_integral(self, nu_min, nu_max, nu_points=50):
+        """Evaluates the integral flux [cm-2 s-1] over a range of frequencies.
+
+        Parameters
+        ----------
+        nu_min : `~astropy.units.Quantity`
+            start frequency (in Hz or equivalent)
+        nu_max : `~astropy.units.Quantity`
+            end frequency (in Hz or equivalent)
+        nu_points: int
+            number of points (between nu_min and nu_max) in log scale
+        """
+        nu = _nu_logspace(nu_min, nu_max, nu_points)
+        energy = nu.to("eV", equivalencies=u.spectral())
+        dphidE = self.diff_energy_flux(nu)
+        return np.trapz(dphidE, energy).to("cm-2 s-1")

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 ..radiative_process import RadiativeProcess
 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(RadiativeProcess):
     """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 ..radiative_process import RadiativeProcess
 
 
 __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(RadiativeProcess):
     """Class for synchrotron radiation computation
 
     Parameters

agnpy/tests/test_radiative_process.py

@@ -0,0 +1,78 @@
+import numpy as np
+import astropy.units as u
+from astropy.constants import m_e
+from astropy.coordinates import Distance
+from agnpy.spectra import PowerLaw
+from agnpy.emission_regions import Blob
+from agnpy.synchrotron import Synchrotron
+from agnpy.radiative_process import RadiativeProcess
+
+
+def estimate_powerlaw_integral(x, y):
+    """Fit a power law to the data, then estimate its integral."""
+    index, log10_norm = np.polyfit(np.log10(x), np.log10(y), 1)
+    norm = 10**log10_norm
+    integral_index = index + 1
+    integral = (
+        norm / (integral_index) * (x[-1] ** integral_index - x[0] ** integral_index)
+    )
+    return integral
+
+
+def create_synchrotron_model():
+    """A simple synchrotron model to be used for test purpoises.
+    Based on the synchrotron example from the documentation."""
+    R_b = 1e16 * u.cm
+    V_b = 4 / 3 * np.pi * R_b**3
+    z = Distance(1e27, unit=u.cm).z
+    B = 1 * u.G
+    delta_D = 10
+    n_e = PowerLaw.from_total_energy(
+        1e48 * u.erg, V_b, p=2.8, gamma_min=1e2, gamma_max=1e7, mass=m_e
+    )
+    blob = Blob(R_b=R_b, z=z, delta_D=delta_D, Gamma=delta_D, B=B, n_e=n_e)
+    synch = Synchrotron(blob)
+    return synch
+
+
+class FlatFNuSedGenerator(RadiativeProcess):
+    """A dummy generator returning a flat (constant) SED in F_nu."""
+
+    def __init__(self, F_nu):
+        # the constant F_nu value
+        self.F_nu = F_nu
+
+    def sed_flux(self, nu):
+        F_nu = np.ones_like(nu).value * self.F_nu
+        return (nu * F_nu).to("erg cm-2 s-1")
+
+
+class TestSedIntegration:
+    def test_flat_energy_flux_integral(self):
+        """Integrate over flat SED (Fnu equal to 1.0 for all frequencies)."""
+        nu_min = 10 * u.Hz
+        nu_max = 1e9 * u.Hz
+        F_nu = 1 * u.Unit("erg cm-2 s-1 Hz-1")
+        energy_flux_integral = FlatFNuSedGenerator(F_nu).energy_flux_integral(
+            nu_min, nu_max
+        )
+        assert u.isclose(energy_flux_integral, 1e9 * u.Unit("erg cm-2 s-1"))
+
+    def test_synchrotron_integrals(self):
+        """Create an example synchrotron SED, fit the flux with a power-law part
+        and compare its integral."""
+        synch = create_synchrotron_model()
+        nu = np.logspace(13, 20, 50) * u.Hz
+        energy = nu.to("eV", equivalencies=u.spectral())
+
+        energy_flux_integral = synch.energy_flux_integral(nu[0], nu[-1]).value
+        flux_integral = synch.flux_integral(nu[0], nu[-1]).value
+
+        F_nu = synch.sed_flux(nu) / nu
+        dphidE = synch.diff_energy_flux(nu)
+
+        energy_flux_integral_from_fit = estimate_powerlaw_integral(nu.value, F_nu.value)
+        flux_integral_from_fit = estimate_powerlaw_integral(energy.value, dphidE.value)
+
+        assert u.isclose(energy_flux_integral, energy_flux_integral_from_fit, rtol=5e-2)
+        assert u.isclose(flux_integral, flux_integral_from_fit, rtol=5e-2)

environment.yml

@@ -5,13 +5,13 @@ channels:
   - sherpa
 
 dependencies:
-  - python
+  - python=3.9
   - pip
-  - numpy
-  - scipy
-  - astropy
+  - numpy>1.20
+  - scipy!=1.10
+  - astropy>=5.0, <6.0
   - pyyaml # needed to read astropy ecsv file
-  - matplotlib
+  - matplotlib>=3.4
   - sherpa
   - pre-commit
   - pip:

setup.py

@@ -20,6 +20,6 @@
         "License :: OSI Approved :: BSD License",
         "Operating System :: OS Independent",
     ],
-    install_requires=["astropy>=4.0", "numpy>=1.17", "scipy>=1.2", "matplotlib"],
+    install_requires=["astropy>=5.0, <6.0", "numpy>=1.21", "scipy>=1.5,!=1.10", "matplotlib>=3.4"],
     python_requires=">=3.8",
 )