Synchrotron losses time evolution

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.

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)
+
+

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,