diff --git a/agnpy/spectra/spectra.py b/agnpy/spectra/spectra.py index 6fb6f4d..8e2a9eb 100644 --- a/agnpy/spectra/spectra.py +++ b/agnpy/spectra/spectra.py @@ -1,7 +1,8 @@ # module containing the particle spectra import numpy as np import astropy.units as u -from astropy.constants import m_e, m_p +from agnpy.utils.conversion import mec2 +from astropy.constants import m_e, m_p, c from scipy.interpolate import CubicSpline import matplotlib.pyplot as plt @@ -198,6 +199,56 @@ def plot(self, gamma=None, gamma_power=0, ax=None, **kwargs): return ax + def evaluate_time(self, time, energy_loss_function, subintervals_count=10): + """Performs the time evaluation of energy change for this particle distribution. + + Parameters + ---------- + time : `~astropy.units.Quantity` + total time for the calculation + energy_loss_function : function + thr function to be used for calculation of energy loss; + the function accepts the gamma value and produces "energy per time" quantity + subintervals_count : int + optional number defining how many equal-length subintervals the total time will be split into + + Returns + ------- + InterpolatedDistribution + a new distribution + """ + unit_time_interval = time / subintervals_count + + def gamma_recalculated_after_loss(gamma): + old_energy = (gamma * mec2).to("erg") + energy_loss_per_time = energy_loss_function(gamma).to("erg s-1") + energy_loss = (energy_loss_per_time * unit_time_interval).to("erg") + new_energy = old_energy - energy_loss + if np.any(new_energy < 0): + raise ValueError( + "Energy loss formula returned value higher then original energy. Use shorter time ranges.") + new_gamma = (new_energy / mec2).value + return new_gamma + + gammas = np.logspace(np.log10(self.gamma_min), np.log10(self.gamma_max), 200) + n_array = self.__call__(gammas) + + # for each gamma point create a narrow bin, calculate the energy loss for start and end of the bin, + # and scale up density by the bin narrowing factor + for r in range(subintervals_count): + bin_size_factor = 0.0001 + bins_width = gammas * bin_size_factor + bins_end_recalc = gamma_recalculated_after_loss(gammas + bins_width) + gammas = gamma_recalculated_after_loss(gammas) + narrowed_bins = bins_end_recalc - gammas + if np.any(narrowed_bins <= 0): + raise ValueError( + "Energy loss formula returned too big value. Use shorter time ranges.") + density_increase = bins_width / narrowed_bins + n_array = n_array * density_increase + + return InterpolatedDistribution(gammas, n_array) + class PowerLaw(ParticleDistribution): r"""Class describing a power-law particle distribution. diff --git a/agnpy/spectra/tests/test_spectra.py b/agnpy/spectra/tests/test_spectra.py index 8c3332f..dc1f08d 100644 --- a/agnpy/spectra/tests/test_spectra.py +++ b/agnpy/spectra/tests/test_spectra.py @@ -2,8 +2,11 @@ from pathlib import Path import numpy as np import astropy.units as u -from astropy.constants import m_e, m_p +from astropy.constants import m_e, m_p, c import pytest +from astropy.coordinates import Distance + +from agnpy import Blob, Synchrotron from agnpy.spectra import ( PowerLaw, BrokenPowerLaw, @@ -740,3 +743,78 @@ def test_interpolation_physical(self): assert u.allclose( n_e(gamma_init), 2 * n_init, atol=0 * u.Unit("cm-3"), rtol=1e-3 ) + + +class TestSpectraTimeEvolution: + + @pytest.mark.parametrize("p", [0.5, 1.2, 2.4]) + @pytest.mark.parametrize("time_and_steps", [(1, 10), (60, 60), (200, 100)]) + def test_compare_numerical_results_with_analytical_calculation(self, p, time_and_steps): + """Test time evolution of spectral electron density for synchrotron energy losses. + Use a simple power law spectrum for easy calculation of analytical results.""" + time = time_and_steps[0] * u.s + steps = time_and_steps[1] + gamma_min = 1e2 + gamma_max = 1e7 + k = 0.1 * u.Unit("cm-3") + initial_n_e = PowerLaw(k, p, gamma_min=gamma_min, gamma_max=gamma_max, mass=m_e) + blob = Blob(1e16 * u.cm, Distance(1e27, unit=u.cm).z, 0, 10, 1 * u.G, n_e=initial_n_e) + synch = Synchrotron(blob) + evaluated_n_e = initial_n_e.evaluate_time(time, synch.electron_energy_loss_per_time, subintervals_count=steps) + + def gamma_before(gamma_after_time, time): + """Reverse-time calculation of the gamma value before the synchrotron energy loss, + using formula -dE/dt ~ (E ** 2)""" + coef = synch._electron_energy_loss_formula_prefix() / (m_e * c ** 2) + return 1 / ((1 / gamma_after_time) - time * coef) + + def integral_analytical(gamma_min, gamma_max): + """Integral for the power-law distribution""" + return k * (gamma_max ** (1 - p) - gamma_min ** (1 - p)) / (1 - p) + + assert u.isclose( + gamma_before(evaluated_n_e.gamma_max, time), + gamma_max, + rtol=0.05 + ) + assert u.isclose( + evaluated_n_e.integrate(evaluated_n_e.gamma_min, evaluated_n_e.gamma_max), + integral_analytical(gamma_before(evaluated_n_e.gamma_min, time), gamma_before(evaluated_n_e.gamma_max, time)), + rtol=0.05 + ) + assert u.isclose( + # synchrotron losses are highest at highest energy, so test the highest energy range, as the most affected + evaluated_n_e.integrate(evaluated_n_e.gamma_max / 10, evaluated_n_e.gamma_max), + integral_analytical(gamma_before(evaluated_n_e.gamma_max / 10, time), gamma_before(evaluated_n_e.gamma_max, time)), + rtol=0.05 + ) + + def test_compare_calculation_with_calculation_split_into_two(self): + initial_n_e = BrokenPowerLaw( + k=1e-8 * u.Unit("cm-3"), + p1=1.9, + p2=2.6, + gamma_b=1e4, + gamma_min=10, + gamma_max=1e6, + mass=m_e, + ) + blob = Blob(1e16 * u.cm, Distance(1e27, unit=u.cm).z, 0, 10, 1 * u.G, n_e=initial_n_e) + synch = Synchrotron(blob) + + # iterate over 60 s in 20 steps + eval_1 = initial_n_e.evaluate_time(60*u.s, synch.electron_energy_loss_per_time, subintervals_count=20) + # iterate first over 30 s, and then, starting with interpolated distribution, over the remaining 30 s, + # with slightly different number of subintervals + eval_2 = initial_n_e.evaluate_time(30 * u.s, synch.electron_energy_loss_per_time, subintervals_count=10)\ + .evaluate_time(30 * u.s, synch.electron_energy_loss_per_time, subintervals_count=8) + + gamma_min = eval_1.gamma_min + gamma_max = eval_1.gamma_max + gammas = np.logspace(np.log10(gamma_min), np.log10(gamma_max)) + assert u.allclose( + eval_1.evaluate(gammas, 1, gamma_min, gamma_max), + eval_2.evaluate(gammas, 1, gamma_min, gamma_max), + 0.001) + + diff --git a/agnpy/synchrotron/synchrotron.py b/agnpy/synchrotron/synchrotron.py index 9289030..cccab94 100644 --- a/agnpy/synchrotron/synchrotron.py +++ b/agnpy/synchrotron/synchrotron.py @@ -1,7 +1,7 @@ # module containing the synchrotron radiative process import numpy as np import astropy.units as u -from astropy.constants import e, h, c, m_e, sigma_T +from astropy.constants import e, h, c, m_e, sigma_T, mu0 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 @@ -90,6 +90,19 @@ def __init__(self, blob, ssa=False, integrator=np.trapz): self.ssa = ssa self.integrator = integrator + def electron_energy_loss_per_time(self, gamma): + # For synchrotron, energy loss dE/dt formula is: + # (4/3 * sigma_T * c * magn_energy_dens) * gamma^2 + # In case of constant B, the part in brackets is fixed, so as an improvement we could calculate it once + # and cache, and later we could use the formula (fixed * gamma^2) + fixed = self._electron_energy_loss_formula_prefix() + return fixed * gamma ** 2 + + def _electron_energy_loss_formula_prefix(self): + # using SI units here because of the ambiguity of the different CGS systems in terms of mu0 definition + magn_energy_dens = (self.blob.B.to("T") ** 2) / (2 * mu0) + return ((4 / 3) * sigma_T * c * magn_energy_dens).to("erg/s") + @staticmethod def evaluate_tau_ssa( nu, diff --git a/environment.yml b/environment.yml index f2bfd43..6fc9f68 100644 --- a/environment.yml +++ b/environment.yml @@ -5,15 +5,12 @@ channels: - sherpa dependencies: - - python=3.9 - - pip - - numpy>1.20 - - scipy!=1.10 - - astropy>=5.0, <6.0 + - python>=3.9,<3.12 + - astropy>=5.0,<6.0 + - numpy>=1.21 + - scipy>=1.5,<1.10 - pyyaml # needed to read astropy ecsv file - matplotlib>=3.4 - sherpa - pre-commit - - pip: - - agnpy - - gammapy + - gammapy<1.2 diff --git a/setup.py b/setup.py index 683badf..eb7115c 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, <6.0", "numpy>=1.21", "scipy>=1.5,!=1.10", "matplotlib>=3.4"], - python_requires=">=3.8", + install_requires=["astropy>=5.0,<6.0", "numpy>=1.21", "scipy>=1.5,<1.10", "pyyaml", "matplotlib>=3.4", "sherpa", "pre-commit", "gammapy<1.2"], + python_requires=">=3.9,<3.12", )