Merge pull request #157 from grzegorzbor/synchrotron-time-evolution

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,

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

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",
 )