From 49c504f4fa196113fa912fcbbae87a728ced2d5c Mon Sep 17 00:00:00 2001 From: Grzegorz Borkowski Date: Sun, 5 Nov 2023 19:55:05 +0100 Subject: [PATCH 01/13] Added a mixin for integrating the SED flux for different processes --- agnpy/compton/external_compton.py | 3 ++- agnpy/compton/synchrotron_self_compton.py | 4 ++-- agnpy/synchrotron/proton_synchrotron.py | 4 +++- agnpy/synchrotron/synchrotron.py | 3 ++- agnpy/tests/test_sedintegrable.py | 15 +++++++++++++++ agnpy/utils/sedintegrable.py | 16 ++++++++++++++++ 6 files changed, 40 insertions(+), 5 deletions(-) create mode 100644 agnpy/tests/test_sedintegrable.py create mode 100644 agnpy/utils/sedintegrable.py diff --git a/agnpy/compton/external_compton.py b/agnpy/compton/external_compton.py index 6f6cd8d..88d5a8d 100644 --- a/agnpy/compton/external_compton.py +++ b/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 diff --git a/agnpy/compton/synchrotron_self_compton.py b/agnpy/compton/synchrotron_self_compton.py index 640e005..8c2499a 100644 --- a/agnpy/compton/synchrotron_self_compton.py +++ b/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 diff --git a/agnpy/synchrotron/proton_synchrotron.py b/agnpy/synchrotron/proton_synchrotron.py index 56c74c0..014f1bb 100644 --- a/agnpy/synchrotron/proton_synchrotron.py +++ b/agnpy/synchrotron/proton_synchrotron.py @@ -4,6 +4,7 @@ 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"] @@ -11,7 +12,8 @@ e = e.gauss B_cr = 4.414e13 * u.G # critical magnetic field -class ProtonSynchrotron: + +class ProtonSynchrotron(SedFluxIntegrable): """Class for synchrotron radiation computation Parameters diff --git a/agnpy/synchrotron/synchrotron.py b/agnpy/synchrotron/synchrotron.py index 0469449..2c85a52 100644 --- a/agnpy/synchrotron/synchrotron.py +++ b/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 diff --git a/agnpy/tests/test_sedintegrable.py b/agnpy/tests/test_sedintegrable.py new file mode 100644 index 0000000..4100ac2 --- /dev/null +++ b/agnpy/tests/test_sedintegrable.py @@ -0,0 +1,15 @@ +import numpy as np +import astropy.units as u +from ..utils.sedintegrable import SedFluxIntegrable + + +class TestSedIntegrable: + + def test_flatintegral(self): + """Integrate over flat SED (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) + + flux = FlatSedGenerator().integrate_flux(np.logspace(1, 20) * u.Hz) + assert flux.value == 1e+20 diff --git a/agnpy/utils/sedintegrable.py b/agnpy/utils/sedintegrable.py new file mode 100644 index 0000000..0de6824 --- /dev/null +++ b/agnpy/utils/sedintegrable.py @@ -0,0 +1,16 @@ +import numpy as np +import astropy.units as u + + +class SedFluxIntegrable: + def integrate_flux(self, nu): + r""" Evaluates the SED flux integral over the span of provided frequencies + + Parameters + ---------- + nu : :class:`~astropy.units.Quantity` + array of frequencies, in Hz, to compute the SED + **note** these are observed frequencies (observer frame) + """ + sed = self.sed_flux(nu) + return np.trapz(sed, nu) # should i rather use trapz_loglog from utils.math? From f01732d4fef4414fdba8261d7ad816f99b84d3b1 Mon Sep 17 00:00:00 2001 From: Grzegorz Borkowski Date: Fri, 10 Nov 2023 13:10:11 +0100 Subject: [PATCH 02/13] Fixed the logic for integrating over dlog(nu); added a test for more realistic data (based on the synchrotron radiation sample from the agnpy documentation); added the method for photon flux SED integral; changed input parameters (allow any equivalency of Hz) --- agnpy/tests/test_sedintegrable.py | 77 +++++++++++++++++++++++++++++-- agnpy/utils/sedintegrable.py | 36 ++++++++++++--- 2 files changed, 102 insertions(+), 11 deletions(-) diff --git a/agnpy/tests/test_sedintegrable.py b/agnpy/tests/test_sedintegrable.py index 4100ac2..8219c20 100644 --- a/agnpy/tests/test_sedintegrable.py +++ b/agnpy/tests/test_sedintegrable.py @@ -1,15 +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_flatintegral(self): - """Integrate over flat SED (equal to 1.0 for all frequencies""" + 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) + 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 - flux = FlatSedGenerator().integrate_flux(np.logspace(1, 20) * u.Hz) - assert flux.value == 1e+20 + assert math.isclose(actual_integral, estimated_integral, rel_tol=0.05) diff --git a/agnpy/utils/sedintegrable.py b/agnpy/utils/sedintegrable.py index 0de6824..d422709 100644 --- a/agnpy/utils/sedintegrable.py +++ b/agnpy/utils/sedintegrable.py @@ -3,14 +3,38 @@ class SedFluxIntegrable: - def integrate_flux(self, nu): + 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 : :class:`~astropy.units.Quantity` - array of frequencies, in Hz, to compute the SED - **note** these are observed frequencies (observer frame) + 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 """ - sed = self.sed_flux(nu) - return np.trapz(sed, nu) # should i rather use trapz_loglog from utils.math? + 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 From ce999b7c7d8e6b1ab90f24c2d51662fb97735a33 Mon Sep 17 00:00:00 2001 From: Grzegorz Borkowski Date: Sat, 16 Dec 2023 13:31:30 +0100 Subject: [PATCH 03/13] FlatSedGenerator class extracted from the unit test method body --- agnpy/tests/test_sedintegrable.py | 24 +++++++++++++++--------- 1 file changed, 15 insertions(+), 9 deletions(-) diff --git a/agnpy/tests/test_sedintegrable.py b/agnpy/tests/test_sedintegrable.py index 8219c20..c3f00ce 100644 --- a/agnpy/tests/test_sedintegrable.py +++ b/agnpy/tests/test_sedintegrable.py @@ -18,7 +18,7 @@ def _estimate_powerlaw_integral(x_values, y_values): def _sample_synchrotron_model(): - # Based on the synchrotron example from agnpy documentation + """The model based on the synchrotron example from the 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) @@ -27,22 +27,28 @@ def _sample_synchrotron_model(): return synch +class FlatSedGenerator(SedFluxIntegrable): + """A dummy generator returning flat (constant) flux, of the same value across the whole spectrum.""" + def __init__(self, flat_value): + self.flat_value = flat_value + + def sed_flux(self, nu): + return np.full(fill_value=self.flat_value, shape=nu.shape, dtype=np.float64) * ( + u.erg / (u.cm ** 2 * u.s * u.Hz)) * nu + + 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 - + """Integrate over flat SED (Fnu equal to 1.0 for all frequencies).""" nu_min = 10 * u.Hz nu_max = 10 ** 20 * u.Hz - flux = FlatSedGenerator().integrate_flux(nu_min, nu_max) + flux = FlatSedGenerator(1.0).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""" + """Integrate over sample synchrotron flux and compare with the value estimated using the power law integral.""" synch = _sample_synchrotron_model() nu_min_log = 13 @@ -59,7 +65,7 @@ def test_synchrotron_energy_flux_integral(self): 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""" + """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 From 259c9b3b00441104de36afc016536721e99ffdcd Mon Sep 17 00:00:00 2001 From: Grzegorz Borkowski Date: Sat, 16 Dec 2023 14:12:30 +0100 Subject: [PATCH 04/13] SedIntegrable class renamed to RadiativeProcess --- agnpy/compton/external_compton.py | 4 ++-- agnpy/compton/synchrotron_self_compton.py | 4 ++-- agnpy/radiation/__init__.py | 0 .../sedintegrable.py => radiation/radiative_process.py} | 4 +++- agnpy/synchrotron/proton_synchrotron.py | 4 ++-- agnpy/synchrotron/synchrotron.py | 4 ++-- .../{test_sedintegrable.py => test_sed_integration.py} | 6 +++--- 7 files changed, 14 insertions(+), 12 deletions(-) create mode 100644 agnpy/radiation/__init__.py rename agnpy/{utils/sedintegrable.py => radiation/radiative_process.py} (87%) rename agnpy/tests/{test_sedintegrable.py => test_sed_integration.py} (96%) diff --git a/agnpy/compton/external_compton.py b/agnpy/compton/external_compton.py index 88d5a8d..d42aa71 100644 --- a/agnpy/compton/external_compton.py +++ b/agnpy/compton/external_compton.py @@ -9,7 +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 agnpy.radiation.radiative_process import RadiativeProcess from ..targets import ( CMB, PointSourceBehindJet, @@ -22,7 +22,7 @@ __all__ = ["ExternalCompton"] -class ExternalCompton(SedFluxIntegrable): +class ExternalCompton(RadiativeProcess): """class for External Compton radiation computation Parameters diff --git a/agnpy/compton/synchrotron_self_compton.py b/agnpy/compton/synchrotron_self_compton.py index 8c2499a..0b0c168 100644 --- a/agnpy/compton/synchrotron_self_compton.py +++ b/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 +from agnpy.radiation.radiative_process import RadiativeProcess __all__ = ["SynchrotronSelfCompton"] -class SynchrotronSelfCompton(SedFluxIntegrable): +class SynchrotronSelfCompton(RadiativeProcess): """class for Synchrotron Self Compton radiation computation Parameters diff --git a/agnpy/radiation/__init__.py b/agnpy/radiation/__init__.py new file mode 100644 index 0000000..e69de29 diff --git a/agnpy/utils/sedintegrable.py b/agnpy/radiation/radiative_process.py similarity index 87% rename from agnpy/utils/sedintegrable.py rename to agnpy/radiation/radiative_process.py index d422709..202a7f8 100644 --- a/agnpy/utils/sedintegrable.py +++ b/agnpy/radiation/radiative_process.py @@ -2,7 +2,8 @@ import astropy.units as u -class SedFluxIntegrable: +class RadiativeProcess: + """Base class for radiative processes. Contains common logic, e.g. calculating SED flux integrals.""" def integrate_flux(self, nu_min, nu_max, nu_points=50): r""" Evaluates the SED flux integral over the span of provided frequencies @@ -34,6 +35,7 @@ def integrate_photon_flux(self, nu_min, nu_max, nu_points=50): return np.trapz(n, nu) # photons / s / cm2 def _nu_logspace(self, nu_min, nu_max, nu_points=50): + """A thin wrapper around numpy.logspace() function, capable of spectral-equivalency conversion of input parameters.""" 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 diff --git a/agnpy/synchrotron/proton_synchrotron.py b/agnpy/synchrotron/proton_synchrotron.py index 014f1bb..d81e4fb 100644 --- a/agnpy/synchrotron/proton_synchrotron.py +++ b/agnpy/synchrotron/proton_synchrotron.py @@ -4,7 +4,7 @@ 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 agnpy.radiation.radiative_process import RadiativeProcess from .synchrotron import single_particle_synch_power, tau_to_attenuation __all__ = ["ProtonSynchrotron"] @@ -13,7 +13,7 @@ B_cr = 4.414e13 * u.G # critical magnetic field -class ProtonSynchrotron(SedFluxIntegrable): +class ProtonSynchrotron(RadiativeProcess): """Class for synchrotron radiation computation Parameters diff --git a/agnpy/synchrotron/synchrotron.py b/agnpy/synchrotron/synchrotron.py index 2c85a52..8391bed 100644 --- a/agnpy/synchrotron/synchrotron.py +++ b/agnpy/synchrotron/synchrotron.py @@ -4,7 +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 +from agnpy.radiation.radiative_process import RadiativeProcess __all__ = ["R", "nu_synch_peak", "Synchrotron"] @@ -68,7 +68,7 @@ def tau_to_attenuation(tau): return np.where(tau < 1e-3, 1, 3 * u / tau) -class Synchrotron(SedFluxIntegrable): +class Synchrotron(RadiativeProcess): """Class for synchrotron radiation computation Parameters diff --git a/agnpy/tests/test_sedintegrable.py b/agnpy/tests/test_sed_integration.py similarity index 96% rename from agnpy/tests/test_sedintegrable.py rename to agnpy/tests/test_sed_integration.py index c3f00ce..c7f42ab 100644 --- a/agnpy/tests/test_sedintegrable.py +++ b/agnpy/tests/test_sed_integration.py @@ -6,7 +6,7 @@ from agnpy.spectra import PowerLaw from agnpy.emission_regions import Blob from agnpy.synchrotron import Synchrotron -from ..utils.sedintegrable import SedFluxIntegrable +from agnpy.radiation.radiative_process import RadiativeProcess def _estimate_powerlaw_integral(x_values, y_values): @@ -27,7 +27,7 @@ def _sample_synchrotron_model(): return synch -class FlatSedGenerator(SedFluxIntegrable): +class FlatSedGenerator(RadiativeProcess): """A dummy generator returning flat (constant) flux, of the same value across the whole spectrum.""" def __init__(self, flat_value): self.flat_value = flat_value @@ -37,7 +37,7 @@ def sed_flux(self, nu): u.erg / (u.cm ** 2 * u.s * u.Hz)) * nu -class TestSedIntegrable: +class TestSedIntegration: def test_flat_integral(self): """Integrate over flat SED (Fnu equal to 1.0 for all frequencies).""" From 5cd33e42f00f6e049f8ed738acc5cc94a6f6ccbc Mon Sep 17 00:00:00 2001 From: Grzegorz Borkowski Date: Sat, 16 Dec 2023 20:11:51 +0100 Subject: [PATCH 05/13] extracted diff_photon_flux() function to RadiativeProcess base class --- agnpy/radiation/radiative_process.py | 21 +++++++++++++++------ agnpy/tests/test_sed_integration.py | 10 +++------- 2 files changed, 18 insertions(+), 13 deletions(-) diff --git a/agnpy/radiation/radiative_process.py b/agnpy/radiation/radiative_process.py index 202a7f8..1f2f436 100644 --- a/agnpy/radiation/radiative_process.py +++ b/agnpy/radiation/radiative_process.py @@ -1,5 +1,6 @@ import numpy as np import astropy.units as u +from astropy.constants import h class RadiativeProcess: @@ -19,7 +20,7 @@ def integrate_flux(self, nu_min, nu_max, nu_points=50): 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 + r""" Evaluates the photon flux integral (photons/s*cm2) from SED over the span of provided frequencies Parameters ---------- @@ -28,11 +29,19 @@ def integrate_photon_flux(self, nu_min, nu_max, nu_points=50): 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 + flux = self.diff_photon_flux(nu) + energy = nu.to("eV", equivalencies=u.spectral()) + return np.trapz(flux, energy) + + def diff_photon_flux(self, nu): + """Similar to sed_flux(), but returns the differential photon flux in energy (photons/s*cm2*eV).""" + photon_energy_eV = nu.to("eV", equivalencies=u.spectral()) + sed_nufnu = self.sed_flux(nu) # erg / s / cm2 / log(nu) + sed_fnu = sed_nufnu / nu # erg / s / cm2 / Hz + sed_fnu_eV = sed_fnu.to("eV / (cm2 s Hz)") # eV / s / cm2 / Hz + photons_per_Hz = sed_fnu_eV / photon_energy_eV # photons / s / cm2 / Hz + photons_per_eV = photons_per_Hz / h.to("eV/Hz") # photons / s / cm2 / eV + return photons_per_eV def _nu_logspace(self, nu_min, nu_max, nu_points=50): """A thin wrapper around numpy.logspace() function, capable of spectral-equivalency conversion of input parameters.""" diff --git a/agnpy/tests/test_sed_integration.py b/agnpy/tests/test_sed_integration.py index c7f42ab..c3e3b8b 100644 --- a/agnpy/tests/test_sed_integration.py +++ b/agnpy/tests/test_sed_integration.py @@ -74,14 +74,10 @@ def test_synchrotron_photon_flux_integral(self): 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 + flux = synch.diff_photon_flux(nu) + energies_eV = nu.to("eV", equivalencies=u.spectral()) - estimated_integral = _estimate_powerlaw_integral(energies_eV, photons_per_eV) + estimated_integral = _estimate_powerlaw_integral(energies_eV, flux) actual_integral = synch.integrate_photon_flux(nu_min, nu_max).value From 49ec93cbe8ba7b9ed6f2e57bec9a08aff7e0ac88 Mon Sep 17 00:00:00 2001 From: cosimoNigro Date: Thu, 4 Jan 2024 16:51:54 +0100 Subject: [PATCH 06/13] fixed pre-commit issue --- agnpy/radiation/__init__.py | 0 agnpy/radiation/radiative_process.py | 51 ---------------- agnpy/radiative_process/__init__.py | 1 + agnpy/radiative_process/core.py | 59 +++++++++++++++++++ ...tegration.py => test_radiative_process.py} | 0 5 files changed, 60 insertions(+), 51 deletions(-) delete mode 100644 agnpy/radiation/__init__.py delete mode 100644 agnpy/radiation/radiative_process.py create mode 100644 agnpy/radiative_process/__init__.py create mode 100644 agnpy/radiative_process/core.py rename agnpy/tests/{test_sed_integration.py => test_radiative_process.py} (100%) diff --git a/agnpy/radiation/__init__.py b/agnpy/radiation/__init__.py deleted file mode 100644 index e69de29..0000000 diff --git a/agnpy/radiation/radiative_process.py b/agnpy/radiation/radiative_process.py deleted file mode 100644 index 1f2f436..0000000 --- a/agnpy/radiation/radiative_process.py +++ /dev/null @@ -1,51 +0,0 @@ -import numpy as np -import astropy.units as u -from astropy.constants import h - - -class RadiativeProcess: - """Base class for radiative processes. Contains common logic, e.g. calculating SED flux integrals.""" - 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 (photons/s*cm2) 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) - flux = self.diff_photon_flux(nu) - energy = nu.to("eV", equivalencies=u.spectral()) - return np.trapz(flux, energy) - - def diff_photon_flux(self, nu): - """Similar to sed_flux(), but returns the differential photon flux in energy (photons/s*cm2*eV).""" - photon_energy_eV = nu.to("eV", equivalencies=u.spectral()) - sed_nufnu = self.sed_flux(nu) # erg / s / cm2 / log(nu) - sed_fnu = sed_nufnu / nu # erg / s / cm2 / Hz - sed_fnu_eV = sed_fnu.to("eV / (cm2 s Hz)") # eV / s / cm2 / Hz - photons_per_Hz = sed_fnu_eV / photon_energy_eV # photons / s / cm2 / Hz - photons_per_eV = photons_per_Hz / h.to("eV/Hz") # photons / s / cm2 / eV - return photons_per_eV - - def _nu_logspace(self, nu_min, nu_max, nu_points=50): - """A thin wrapper around numpy.logspace() function, capable of spectral-equivalency conversion of input parameters.""" - 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 diff --git a/agnpy/radiative_process/__init__.py b/agnpy/radiative_process/__init__.py new file mode 100644 index 0000000..8b13789 --- /dev/null +++ b/agnpy/radiative_process/__init__.py @@ -0,0 +1 @@ + diff --git a/agnpy/radiative_process/core.py b/agnpy/radiative_process/core.py new file mode 100644 index 0000000..6cad33c --- /dev/null +++ b/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") diff --git a/agnpy/tests/test_sed_integration.py b/agnpy/tests/test_radiative_process.py similarity index 100% rename from agnpy/tests/test_sed_integration.py rename to agnpy/tests/test_radiative_process.py From a34f12e9e85837f83fd75665234ef4c71e67ab4b Mon Sep 17 00:00:00 2001 From: cosimoNigro Date: Thu, 4 Jan 2024 16:53:50 +0100 Subject: [PATCH 07/13] improved tests on the RadiativeProcess class --- agnpy/radiative_process/__init__.py | 2 +- agnpy/tests/test_radiative_process.py | 118 ++++++++++++-------------- 2 files changed, 57 insertions(+), 63 deletions(-) diff --git a/agnpy/radiative_process/__init__.py b/agnpy/radiative_process/__init__.py index 8b13789..bb67a43 100644 --- a/agnpy/radiative_process/__init__.py +++ b/agnpy/radiative_process/__init__.py @@ -1 +1 @@ - +from .core import * diff --git a/agnpy/tests/test_radiative_process.py b/agnpy/tests/test_radiative_process.py index c3e3b8b..4b044fe 100644 --- a/agnpy/tests/test_radiative_process.py +++ b/agnpy/tests/test_radiative_process.py @@ -1,84 +1,78 @@ import numpy as np -import math import astropy.units as u -from astropy.constants import m_e, h +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.radiation.radiative_process import RadiativeProcess +from agnpy.radiative_process import RadiativeProcess -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 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 _sample_synchrotron_model(): - """The model based on the synchrotron example from the agnpy documentation.""" +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 - 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) + 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 FlatSedGenerator(RadiativeProcess): - """A dummy generator returning flat (constant) flux, of the same value across the whole spectrum.""" - def __init__(self, flat_value): - self.flat_value = flat_value +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): - return np.full(fill_value=self.flat_value, shape=nu.shape, dtype=np.float64) * ( - u.erg / (u.cm ** 2 * u.s * u.Hz)) * 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_integral(self): + 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 = 10 ** 20 * u.Hz - flux = FlatSedGenerator(1.0).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 - flux = synch.diff_photon_flux(nu) - energies_eV = nu.to("eV", equivalencies=u.spectral()) - - estimated_integral = _estimate_powerlaw_integral(energies_eV, flux) - - actual_integral = synch.integrate_photon_flux(nu_min, nu_max).value - - assert math.isclose(actual_integral, estimated_integral, rel_tol=0.05) + 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) From e830730c2b14c8186caf96a709d17c70517688d2 Mon Sep 17 00:00:00 2001 From: cosimoNigro Date: Thu, 4 Jan 2024 16:54:36 +0100 Subject: [PATCH 08/13] modified all radiative classes to inherit RadiativeProcess --- agnpy/compton/external_compton.py | 2 +- agnpy/compton/synchrotron_self_compton.py | 2 +- agnpy/synchrotron/proton_synchrotron.py | 2 +- agnpy/synchrotron/synchrotron.py | 2 +- 4 files changed, 4 insertions(+), 4 deletions(-) diff --git a/agnpy/compton/external_compton.py b/agnpy/compton/external_compton.py index d42aa71..692c24e 100644 --- a/agnpy/compton/external_compton.py +++ b/agnpy/compton/external_compton.py @@ -9,7 +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 agnpy.radiation.radiative_process import RadiativeProcess +from ..radiative_process import RadiativeProcess from ..targets import ( CMB, PointSourceBehindJet, diff --git a/agnpy/compton/synchrotron_self_compton.py b/agnpy/compton/synchrotron_self_compton.py index 0b0c168..32bbe67 100644 --- a/agnpy/compton/synchrotron_self_compton.py +++ b/agnpy/compton/synchrotron_self_compton.py @@ -10,7 +10,7 @@ nu_to_integrate, ) from ..utils.conversion import nu_to_epsilon_prime -from agnpy.radiation.radiative_process import RadiativeProcess +from ..radiative_process import RadiativeProcess __all__ = ["SynchrotronSelfCompton"] diff --git a/agnpy/synchrotron/proton_synchrotron.py b/agnpy/synchrotron/proton_synchrotron.py index d81e4fb..4c0c5b6 100644 --- a/agnpy/synchrotron/proton_synchrotron.py +++ b/agnpy/synchrotron/proton_synchrotron.py @@ -4,7 +4,7 @@ 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 agnpy.radiation.radiative_process import RadiativeProcess +from ..radiative_process import RadiativeProcess from .synchrotron import single_particle_synch_power, tau_to_attenuation __all__ = ["ProtonSynchrotron"] diff --git a/agnpy/synchrotron/synchrotron.py b/agnpy/synchrotron/synchrotron.py index 8391bed..9289030 100644 --- a/agnpy/synchrotron/synchrotron.py +++ b/agnpy/synchrotron/synchrotron.py @@ -4,7 +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 agnpy.radiation.radiative_process import RadiativeProcess +from ..radiative_process import RadiativeProcess __all__ = ["R", "nu_synch_peak", "Synchrotron"] From 31bee7da6554ac0fb42a54fdbbd6fa68a561a87e Mon Sep 17 00:00:00 2001 From: cosimoNigro Date: Thu, 4 Jan 2024 17:02:07 +0100 Subject: [PATCH 09/13] tried to fix dependenies to make tests pass --- setup.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/setup.py b/setup.py index e6958e3..0da168f 100644 --- a/setup.py +++ b/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", "numpy>=1.21", "scipy>=1.5,!=1.10", "matplotlib>=3.4"], python_requires=">=3.8", ) From ce517cf65a7a221f128086a56a88fde14e312ae9 Mon Sep 17 00:00:00 2001 From: cosimoNigro Date: Thu, 4 Jan 2024 17:04:47 +0100 Subject: [PATCH 10/13] fixed dependencies also in the environment.yml --- environment.yml | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/environment.yml b/environment.yml index 7a54a17..58ce293 100644 --- a/environment.yml +++ b/environment.yml @@ -5,13 +5,13 @@ channels: - sherpa dependencies: - - python + - python=3.9 - pip - - numpy - - scipy + - numpy>1.20 + - scipy!=1.10 - astropy - pyyaml # needed to read astropy ecsv file - - matplotlib + - matplotlib>=3.4 - sherpa - pre-commit - pip: From 1a2990d44e94734cee0125508c36c93749b9b7ca Mon Sep 17 00:00:00 2001 From: cosimoNigro Date: Thu, 4 Jan 2024 17:07:33 +0100 Subject: [PATCH 11/13] removed python 3.8 from the tests --- .github/workflows/test.yml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/.github/workflows/test.yml b/.github/workflows/test.yml index 8502b08..c56f2d9 100644 --- a/.github/workflows/test.yml +++ b/.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] steps: - uses: actions/checkout@v3 From 3fefaeff8e567aa60584bdee974f13824683332e Mon Sep 17 00:00:00 2001 From: cosimoNigro Date: Thu, 4 Jan 2024 18:42:25 +0100 Subject: [PATCH 12/13] 3.10 -> 3.10.10 --- .github/workflows/test.yml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/.github/workflows/test.yml b/.github/workflows/test.yml index c56f2d9..f49592b 100644 --- a/.github/workflows/test.yml +++ b/.github/workflows/test.yml @@ -12,7 +12,7 @@ jobs: runs-on: ubuntu-latest strategy: matrix: - python-version: [3.9, 3.10] + python-version: [3.9, 3.10.10] steps: - uses: actions/checkout@v3 From be47b3f8f55097c4ab36b0ceba557e4313acfe82 Mon Sep 17 00:00:00 2001 From: cosimoNigro Date: Thu, 4 Jan 2024 18:49:58 +0100 Subject: [PATCH 13/13] request astropy < 6.0 --- environment.yml | 2 +- setup.py | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/environment.yml b/environment.yml index 58ce293..f2bfd43 100644 --- a/environment.yml +++ b/environment.yml @@ -9,7 +9,7 @@ dependencies: - pip - numpy>1.20 - scipy!=1.10 - - astropy + - astropy>=5.0, <6.0 - pyyaml # needed to read astropy ecsv file - matplotlib>=3.4 - sherpa diff --git a/setup.py b/setup.py index 0da168f..683badf 100644 --- a/setup.py +++ b/setup.py @@ -20,6 +20,6 @@ "License :: OSI Approved :: BSD License", "Operating System :: OS Independent", ], - install_requires=["astropy>=5.0", "numpy>=1.21", "scipy>=1.5,!=1.10", "matplotlib>=3.4"], + install_requires=["astropy>=5.0, <6.0", "numpy>=1.21", "scipy>=1.5,!=1.10", "matplotlib>=3.4"], python_requires=">=3.8", )