From cca593da35557e46aeaae515ce59a8f18e7adf23 Mon Sep 17 00:00:00 2001 From: cosimoNigro Date: Tue, 16 May 2023 10:57:22 +0200 Subject: [PATCH 1/3] started to add physical particle distributions --- agnpy/spectra/__init__.py | 3 +- agnpy/spectra/spectra.py | 582 -------------------------------------- 2 files changed, 2 insertions(+), 583 deletions(-) diff --git a/agnpy/spectra/__init__.py b/agnpy/spectra/__init__.py index 2d3689f2..262de182 100644 --- a/agnpy/spectra/__init__.py +++ b/agnpy/spectra/__init__.py @@ -1 +1,2 @@ -from .spectra import * +from .analytical import * +from .physical import * diff --git a/agnpy/spectra/spectra.py b/agnpy/spectra/spectra.py index fc77db64..5ee724b9 100644 --- a/agnpy/spectra/spectra.py +++ b/agnpy/spectra/spectra.py @@ -2,21 +2,9 @@ import numpy as np import astropy.units as u from astropy.constants import m_e, m_p -from scipy.interpolate import CubicSpline import matplotlib.pyplot as plt -__all__ = [ - "ParticleDistribution", - "PowerLaw", - "ExpCutoffPowerLaw", - "BrokenPowerLaw", - "LogParabola", - "ExpCutoffBrokenPowerLaw", - "InterpolatedDistribution", -] - - class ParticleDistribution: """Base class grouping common functionalities to be used by all particles distributions. @@ -196,573 +184,3 @@ def plot(self, gamma=None, gamma_power=0, ax=None, **kwargs): ) return ax - - -class PowerLaw(ParticleDistribution): - r"""Class describing a power-law particle distribution. - When called, the particle density :math:`n(\gamma)` in :math:`\mathrm{cm}^{-3}` is returned. - - .. math:: - n(\gamma') = k \, \gamma'^{-p} \, H(\gamma'; \gamma'_{\rm min}, \gamma'_{\rm max}) - - Parameters - ---------- - k : :class:`~astropy.units.Quantity` - spectral normalisation - p : float - spectral index, note it is positive by definition, will change sign in the function - gamma_min : float - minimum Lorentz factor of the particle distribution - gamma_max : float - maximum Lorentz factor of the particle distribution - mass : :class:`~astropy.units.Quantity` - particle mass, default is the electron mass - integrator : func - function to be used for integration, default is :class:`~numpy.trapz` - """ - - def __init__( - self, - k=1e-13 * u.Unit("cm-3"), - p=2.1, - gamma_min=10, - gamma_max=1e5, - mass=m_e, - integrator=np.trapz, - ): - super().__init__(mass, integrator) - self.k = k - self.p = p - self.gamma_min = gamma_min - self.gamma_max = gamma_max - - @property - def parameters(self): - return [self.k, self.p, self.gamma_min, self.gamma_max] - - @staticmethod - def evaluate(gamma, k, p, gamma_min, gamma_max): - return np.where( - (gamma_min <= gamma) * (gamma <= gamma_max), k * gamma ** (-p), 0 - ) - - def __call__(self, gamma): - return self.evaluate(gamma, self.k, self.p, self.gamma_min, self.gamma_max) - - @staticmethod - def evaluate_SSA_integrand(gamma, k, p, gamma_min, gamma_max): - r"""Analytical integrand for the synchrotron self-absorption: - :math:`\gamma'^2 \frac{d}{d \gamma'} \left(\frac{n(\gamma)}{\gamma'^2}\right)`.""" - return k * np.where( - (gamma_min <= gamma) * (gamma <= gamma_max), - -(p + 2) * np.power(gamma, -p - 1), - 0, - ) - - def SSA_integrand(self, gamma): - return self.evaluate_SSA_integrand( - gamma, self.k, self.p, self.gamma_min, self.gamma_max - ) - - def __str__(self): - return ( - f"* {self.particle} energy distribution\n" - + f" - power law\n" - + f" - k: {self.k:.2e}\n" - + f" - p: {self.p:.2f}\n" - + f" - gamma_min: {self.gamma_min:.2e}\n" - + f" - gamma_max: {self.gamma_max:.2e}\n" - ) - - -class BrokenPowerLaw(ParticleDistribution): - r"""Class describing a broken power-law particle distribution. - When called, the particle density :math:`n(\gamma)` in :math:`\mathrm{cm}^{-3}` is returned. - - .. math:: - n(\gamma') = k \left[ - \left(\frac{\gamma'}{\gamma'_b}\right)^{-p_1} \, H(\gamma'; \gamma'_{\rm min}, \gamma'_b) + - \left(\frac{\gamma'}{\gamma'_b}\right)^{-p_2} \, H(\gamma'; \gamma'_{b}, \gamma'_{\rm max}) - \right] - - Parameters - ---------- - k : :class:`~astropy.units.Quantity` - spectral normalisation - p1 : float - spectral index before the break (positive by definition) - p2 : float - spectral index after the break (positive by definition) - gamma_b : float - Lorentz factor at which the change in spectral index is occurring - gamma_min : float - minimum Lorentz factor of the particle distribution - gamma_max : float - maximum Lorentz factor of the particle distribution - mass : :class:`~astropy.units.Quantity` - particle mass, default is the electron mass - integrator : func - function to be used for integration, default is :class:`~numpy.trapz` - """ - - def __init__( - self, - k=1e-13 * u.Unit("cm-3"), - p1=2.0, - p2=3.0, - gamma_b=1e3, - gamma_min=10, - gamma_max=1e7, - mass=m_e, - integrator=np.trapz, - ): - super().__init__(mass, integrator) - self.k = k - self.p1 = p1 - self.p2 = p2 - self.gamma_b = gamma_b - self.gamma_min = gamma_min - self.gamma_max = gamma_max - - @property - def parameters(self): - return [ - self.k, - self.p1, - self.p2, - self.gamma_b, - self.gamma_min, - self.gamma_max, - ] - - @staticmethod - def evaluate(gamma, k, p1, p2, gamma_b, gamma_min, gamma_max): - index = np.where(gamma <= gamma_b, p1, p2) - return np.where( - (gamma_min <= gamma) * (gamma <= gamma_max), - k * (gamma / gamma_b) ** (-index), - 0, - ) - - def __call__(self, gamma): - return self.evaluate( - gamma, - self.k, - self.p1, - self.p2, - self.gamma_b, - self.gamma_min, - self.gamma_max, - ) - - @staticmethod - def evaluate_SSA_integrand(gamma, k, p1, p2, gamma_b, gamma_min, gamma_max): - r"""Analytical integrand for the synchrotron self-absorption: - :math:`\gamma'^2 \frac{d}{d \gamma'} \left(\frac{n_e(\gamma)}{\gamma'^2}\right)`.""" - index = np.where(gamma <= gamma_b, p1, p2) - return np.where( - (gamma_min <= gamma) * (gamma <= gamma_max), - k * -(index + 2) / gamma * (gamma / gamma_b) ** (-index), - 0, - ) - - def SSA_integrand(self, gamma): - return self.evaluate_SSA_integrand( - gamma, - self.k, - self.p1, - self.p2, - self.gamma_b, - self.gamma_min, - self.gamma_max, - ) - - def __str__(self): - return ( - f"* {self.particle} energy distribution\n" - + f" - broken power law\n" - + f" - k: {self.k:.2e}\n" - + f" - p1: {self.p1:.2f}\n" - + f" - p2: {self.p2:.2f}\n" - + f" - gamma_b: {self.gamma_b:.2e}\n" - + f" - gamma_min: {self.gamma_min:.2e}\n" - + f" - gamma_max: {self.gamma_max:.2e}\n" - ) - - -class LogParabola(ParticleDistribution): - r"""Class describing a log-parabolic particle distribution. - When called, the particle density :math:`n(\gamma)` in :math:`\mathrm{cm}^{-3}` is returned. - - .. math:: - n(\gamma') = k \, \left(\frac{\gamma'}{\gamma'_0}\right)^{-(p + q \log_{10}(\gamma' / \gamma'_0))} - - Parameters - ---------- - k : :class:`~astropy.units.Quantity` - spectral normalisation - p : float - spectral index, note it is positive by definition, will change sign in the function - q : float - spectral curvature, note it is positive by definition, will change sign in the function - gamma_0 : float - reference Lorentz factor - gamma_min : float - minimum Lorentz factor of the particle distribution - gamma_max : float - maximum Lorentz factor of the particle distribution - mass : :class:`~astropy.units.Quantity` - particle mass, default is the electron mass - integrator : func - function to be used for integration, default is :class:`~numpy.trapz` - """ - - def __init__( - self, - k=1e-13 * u.Unit("cm-3"), - p=2.0, - q=0.1, - gamma_0=1e3, - gamma_min=10, - gamma_max=1e7, - mass=m_e, - integrator=np.trapz, - ): - super().__init__(mass, integrator) - self.k = k - self.p = p - self.q = q - self.gamma_0 = gamma_0 - self.gamma_min = gamma_min - self.gamma_max = gamma_max - - @property - def parameters(self): - return [self.k, self.p, self.q, self.gamma_0, self.gamma_min, self.gamma_max] - - @staticmethod - def evaluate(gamma, k, p, q, gamma_0, gamma_min, gamma_max): - gamma_ratio = gamma / gamma_0 - index = -p - q * np.log10(gamma_ratio) - return np.where( - (gamma_min <= gamma) * (gamma <= gamma_max), k * gamma_ratio ** index, 0, - ) - - def __call__(self, gamma): - return self.evaluate( - gamma, self.k, self.p, self.q, self.gamma_0, self.gamma_min, self.gamma_max, - ) - - @staticmethod - def evaluate_SSA_integrand(gamma, k, p, q, gamma_0, gamma_min, gamma_max): - r"""Analytical integrand for the synchrotron self-absorption: - :math:`\gamma'^2 \frac{d}{d \gamma'} \left(\frac{n_e(\gamma)}{\gamma'^2}\right)`.""" - prefactor = -(p + 2 * q * np.log10(gamma / gamma_0) + 2) / gamma - return prefactor * LogParabola.evaluate( - gamma, k, p, q, gamma_0, gamma_min, gamma_max - ) - - def SSA_integrand(self, gamma): - return self.evaluate_SSA_integrand( - gamma, self.k, self.p, self.q, self.gamma_0, self.gamma_min, self.gamma_max, - ) - - def __str__(self): - return ( - f"* {self.particle} energy distribution\n" - + f" - log parabola\n" - + f" - k: {self.k:.2e}\n" - + f" - p: {self.p:.2f}\n" - + f" - q: {self.q:.2f}\n" - + f" - gamma_0: {self.gamma_0:.2e}\n" - + f" - gamma_min: {self.gamma_min:.2e}\n" - + f" - gamma_max: {self.gamma_max:.2e}\n" - ) - - -class ExpCutoffPowerLaw(ParticleDistribution): - r"""Class describing a power-law with an exponetial cutoff particle distribution. - When called, the particle density :math:`n_e(\gamma)` in :math:`\mathrm{cm}^{-3}` is returned. - - .. math:: - n(\gamma'_c) = k \, \gamma'^{-p} exp(-\gamma'/\gamma_c) \, H(\gamma'; \gamma'{\rm min}, \gamma'{\rm max}) - - Parameters - ---------- - k : :class:`~astropy.units.Quantity` - spectral normalisation - p : float - spectral index, note it is positive by definition, will change sign in the function - gamma_c : float - cutoff Lorentz factor of the particle distribution - gamma_min : float - minimum Lorentz factor of the particle distribution - gamma_max : float - maximum Lorentz factor of the particle distribution - mass : :class:`~astropy.units.Quantity` - particle mass, default is the electron mass - integrator : func - function to be used for integration, default is :class:`~numpy.trapz` - """ - - def __init__( - self, - k=1e-13 * u.Unit("cm-3"), - p=2.1, - gamma_c=1e3, - gamma_min=10, - gamma_max=1e5, - mass=m_e, - integrator=np.trapz, - ): - super().__init__(mass, integrator) - self.k = k - self.p = p - self.gamma_c = gamma_c - self.gamma_min = gamma_min - self.gamma_max = gamma_max - - @property - def parameters(self): - return [self.k, self.p, self.gamma_c, self.gamma_min, self.gamma_max] - - @staticmethod - def evaluate(gamma, k, p, gamma_c, gamma_min, gamma_max): - return np.where( - (gamma_min <= gamma) * (gamma <= gamma_max), - k * gamma ** (-p) * np.exp(-gamma / gamma_c), - 0, - ) - - def __call__(self, gamma): - return self.evaluate( - gamma, self.k, self.p, self.gamma_c, self.gamma_min, self.gamma_max - ) - - @staticmethod - def evaluate_SSA_integrand(gamma, k, p, gamma_c, gamma_min, gamma_max): - r"""(analytical) integrand for the synchrotron self-absorption: - :math:`\gamma'^2 \frac{d}{d \gamma'} \left(\frac{n(\gamma)}{\gamma'^2}\right)`""" - prefactor = -(p + 2) / gamma + (-1 / gamma_c) - - return prefactor * ExpCutoffPowerLaw.evaluate( - gamma, k, p, gamma_c, gamma_min, gamma_max - ) - - def SSA_integrand(self, gamma): - return self.evaluate_SSA_integrand( - gamma, self.k, self.p, self.gamma_c, self.gamma_min, self.gamma_max - ) - - def __str__(self): - return ( - f"* {self.particle} energy distribution\n" - + f" - exponential cut-off power law\n" - + f" - k: {self.k:.2e}\n" - + f" - p: {self.p:.2f}\n" - + f" - gamma_c: {self.gamma_c:.2f}\n" - + f" - gamma_min: {self.gamma_min:.2e}\n" - + f" - gamma_max: {self.gamma_max:.2e}\n" - ) - - -class ExpCutoffBrokenPowerLaw(ParticleDistribution): - r"""Class describing an exponential cutoff broken power-law particle distribution. - When called, the particle density :math:`n(\gamma)` in :math:`\mathrm{cm}^{-3}` is returned. - - .. math:: - n(\gamma'_c) = k \left[ - \left(\frac{\gamma'}{\gamma'_b}\right)^{-p_1} exp(-\gamma'/\gamma_c) \, H(\gamma'; \gamma'_{\rm min}, \gamma'_b) + - \left(\frac{\gamma'}{\gamma'_b}\right)^{-p_2} exp(-\gamma'/\gamma_c)\, H(\gamma'; \gamma'_{b}, \gamma'_{\rm max}) - \right] - - Parameters - ---------- - k : :class:`~astropy.units.Quantity` - spectral normalisation - p1 : float - spectral index before the break (positive by definition) - p2 : float - spectral index after the break (positive by definition) - gamma_b : float - Lorentz factor at which the change in spectral index is occurring - gamma_min : float - minimum Lorentz factor of the particle distribution - gamma_max : float - maximum Lorentz factor of the particle distribution - gamma_c : float - cutoff Lorentz factor of the particle distribution - mass : `~astropy.units.Quantity` - particle mass, default is the electron mass - integrator : func - function to be used for integration, default is :class:`~numpy.trapz` - """ - - def __init__( - self, - k=1e-13 * u.Unit("cm-3"), - p1=2.0, - p2=3.0, - gamma_c = 1e5, - gamma_b=1e3, - gamma_min=10, - gamma_max=1e7, - mass=m_e, - integrator=np.trapz, - ): - super().__init__(mass, integrator) - self.k = k - self.p1 = p1 - self.p2 = p2 - self.gamma_c = gamma_c - self.gamma_b = gamma_b - self.gamma_min = gamma_min - self.gamma_max = gamma_max - - @property - def parameters(self): - return [ - self.k, - self.p1, - self.p2, - self.gamma_c, - self.gamma_b, - self.gamma_min, - self.gamma_max, - ] - - @staticmethod - def evaluate(gamma, k, p1, p2, gamma_c, gamma_b, gamma_min, gamma_max): - index = np.where(gamma <= gamma_b, p1, p2) - return np.where( - (gamma_min <= gamma) * (gamma <= gamma_max), - k * (gamma/gamma_b) ** (-index) * np.exp(-gamma/gamma_c), - 0 - ) - - def __call__(self, gamma): - return self.evaluate( - gamma, - self.k, - self.p1, - self.p2, - self.gamma_c, - self.gamma_b, - self.gamma_min, - self.gamma_max, - ) - - @staticmethod - def evaluate_SSA_integrand(gamma, k, p1, p2, gamma_c, gamma_b, gamma_min, gamma_max): - r"""Analytical integrand for the synchrotron self-absorption: - :math:`\gamma'^2 \frac{d}{d \gamma'} \left(\frac{n_e(\gamma)}{\gamma'^2}\right)`.""" - index = np.where(gamma <= gamma_b, p1, p2) - prefactor = -(index + 2) / gamma + (-1 / gamma_c) - return np.where( - (gamma_min <= gamma) * (gamma <= gamma_max), - prefactor * ExpCutoffBrokenPowerLaw.evaluate( - gamma, k, p1, p2, gamma_c, gamma_b, gamma_min, gamma_max), - 0 - ) - - def SSA_integrand(self, gamma): - return self.evaluate_SSA_integrand( - gamma, - self.k, - self.p1, - self.p2, - self.gamma_c, - self.gamma_b, - self.gamma_min, - self.gamma_max - ) - - def __str__(self): - return ( - f"* {self.particle} energy distribution\n" - + f" - exponential cut-off broken power law\n" - + f" - k: {self.k:.2e}\n" - + f" - p1: {self.p1:.2f}\n" - + f" - p2: {self.p2:.2f}\n" - + f" - gamma_c: {self.gamma_c:.2e}\n" - + f" - gamma_b: {self.gamma_b:.2e}\n" - + f" - gamma_min: {self.gamma_min:.2e}\n" - + f" - gamma_max: {self.gamma_max:.2e}\n" - ) - - -class InterpolatedDistribution(ParticleDistribution): - """Class describing a particle distribution with an arbitrary shape. - The spectrum is interpolated from an array of Lorentz factor and densities. - - Parameters - ---------- - gamma : :class:`~numpy.ndarray` - array of Lorentz factors where the density has been computed - n : :class:`~astropy.units.Quantity` - array of densities to be interpolated - norm : float - parameter to scale the density - mass : :class:`~astropy.units.Quantity` - particle mass, default is the electron mass - integrator : func - function to be used for integration, default is :class:`~numpy.trapz` - """ - - def __init__(self, gamma, n, norm=1, mass=m_e, integrator=np.trapz): - super().__init__(mass, integrator) - self.gamma = gamma - self.gamma_max = np.max(self.gamma) - self.gamma_min = np.min(self.gamma) - self.norm = norm - if n.unit != u.Unit("cm-3"): - raise ValueError( - f"Provide a particle distribution in cm-3, instead of {n.unit}" - ) - else: - self.n = n - # call make the interpolation - self.log10_f = self.log10_interpolation() - - def log10_interpolation(self): - """Returns the function interpolating in log10 the particle spectrum. - TODO: make possible to pass arguments to CubicSpline. - """ - interpolator = CubicSpline( - np.log10(self.gamma), np.log10(self.n.to_value("cm-3")) - ) - return interpolator - - def evaluate(self, gamma, norm, gamma_min, gamma_max): - log10_gamma = np.log10(gamma) - values = np.where( - (gamma_min <= gamma) * (gamma <= gamma_max), - norm * np.power(10, self.log10_f(log10_gamma)), - 0, - ) - return values * u.Unit("cm-3") - - def __call__(self, gamma): - return self.evaluate(gamma, self.norm, self.gamma_min, self.gamma_max) - - def SSA_integrand(self, gamma): - r"""Integrand for the synchrotron self-absorption. It is - - .. math:: - \gamma^2 \frac{d}{d \gamma} (\frac{n_e(\gamma)}{\gamma^2}) = ( \frac{dn_e(\gamma)}{d\gamma}+\frac{2n_e(\gamma)}{\gamma}) - - The derivative is: - - .. math:: - \frac{dn_e(\gamma)}{d\gamma} = \frac{d 10^{f(u(\gamma))}}{d\gamma} = \frac{d10^{f(u)}}{du} \cdot \frac{du(\gamma)}{d\gamma} - - where we have :math:`\frac{d 10^{f(u(\gamma))}}{d\gamma} = \frac{d10^{f(u)}}{du} \cdot \frac{du(\gamma)}{d\gamma}`, - where :math:`u` is the :math:`log_{10}(\gamma)`. - This is equal to :math:`\frac{d 10^{f(u(\gamma))}}{d\gamma} = 10^{f(u)} \cdot \frac{df(u)}{du} \cdot \frac{1}{\gamma}` - - """ - log10_gamma = np.log10(gamma) - df_log = self.log10_f.derivative() - int_fun = self.evaluate(gamma, self.norm, self.gamma_min, self.gamma_max) - deriv = int_fun * (1 / gamma) * df_log(log10_gamma) - return deriv - 2 * int_fun / gamma From eb624576eadad2f0000590f9b5f34590c3278a86 Mon Sep 17 00:00:00 2001 From: cosimoNigro Date: Tue, 16 May 2023 10:57:53 +0200 Subject: [PATCH 2/3] separated particle distributions in analytical and physical --- agnpy/spectra/analytical.py | 586 ++++++++++++++++++++++++++++++++++++ agnpy/spectra/physical.py | 84 ++++++ 2 files changed, 670 insertions(+) create mode 100644 agnpy/spectra/analytical.py create mode 100644 agnpy/spectra/physical.py diff --git a/agnpy/spectra/analytical.py b/agnpy/spectra/analytical.py new file mode 100644 index 00000000..e993bcc9 --- /dev/null +++ b/agnpy/spectra/analytical.py @@ -0,0 +1,586 @@ +# analytical particle energy distributions +import numpy as np +import astropy.units as u +from astropy.constants import m_e +from scipy.interpolate import CubicSpline +from .spectra import ParticleDistribution + + +__all__ = [ + "PowerLaw", + "ExpCutoffPowerLaw", + "BrokenPowerLaw", + "LogParabola", + "ExpCutoffBrokenPowerLaw", + "InterpolatedDistribution", +] + + +class PowerLaw(ParticleDistribution): + r"""Class describing a power-law particle distribution. + When called, the particle density :math:`n(\gamma)` in :math:`\mathrm{cm}^{-3}` is returned. + + .. math:: + n(\gamma') = k \, \gamma'^{-p} \, H(\gamma'; \gamma'_{\rm min}, \gamma'_{\rm max}) + + Parameters + ---------- + k : :class:`~astropy.units.Quantity` + spectral normalisation + p : float + spectral index, note it is positive by definition, will change sign in the function + gamma_min : float + minimum Lorentz factor of the particle distribution + gamma_max : float + maximum Lorentz factor of the particle distribution + mass : :class:`~astropy.units.Quantity` + particle mass, default is the electron mass + integrator : func + function to be used for integration, default is :class:`~numpy.trapz` + """ + + def __init__( + self, + k=1e-13 * u.Unit("cm-3"), + p=2.1, + gamma_min=10, + gamma_max=1e5, + mass=m_e, + integrator=np.trapz, + ): + super().__init__(mass, integrator) + self.k = k + self.p = p + self.gamma_min = gamma_min + self.gamma_max = gamma_max + + @property + def parameters(self): + return [self.k, self.p, self.gamma_min, self.gamma_max] + + @staticmethod + def evaluate(gamma, k, p, gamma_min, gamma_max): + return np.where( + (gamma_min <= gamma) * (gamma <= gamma_max), k * gamma ** (-p), 0 + ) + + def __call__(self, gamma): + return self.evaluate(gamma, self.k, self.p, self.gamma_min, self.gamma_max) + + @staticmethod + def evaluate_SSA_integrand(gamma, k, p, gamma_min, gamma_max): + r"""Analytical integrand for the synchrotron self-absorption: + :math:`\gamma'^2 \frac{d}{d \gamma'} \left(\frac{n(\gamma)}{\gamma'^2}\right)`.""" + return k * np.where( + (gamma_min <= gamma) * (gamma <= gamma_max), + -(p + 2) * np.power(gamma, -p - 1), + 0, + ) + + def SSA_integrand(self, gamma): + return self.evaluate_SSA_integrand( + gamma, self.k, self.p, self.gamma_min, self.gamma_max + ) + + def __str__(self): + return ( + f"* {self.particle} energy distribution\n" + + f" - power law\n" + + f" - k: {self.k:.2e}\n" + + f" - p: {self.p:.2f}\n" + + f" - gamma_min: {self.gamma_min:.2e}\n" + + f" - gamma_max: {self.gamma_max:.2e}\n" + ) + + +class BrokenPowerLaw(ParticleDistribution): + r"""Class describing a broken power-law particle distribution. + When called, the particle density :math:`n(\gamma)` in :math:`\mathrm{cm}^{-3}` is returned. + + .. math:: + n(\gamma') = k \left[ + \left(\frac{\gamma'}{\gamma'_b}\right)^{-p_1} \, H(\gamma'; \gamma'_{\rm min}, \gamma'_b) + + \left(\frac{\gamma'}{\gamma'_b}\right)^{-p_2} \, H(\gamma'; \gamma'_{b}, \gamma'_{\rm max}) + \right] + + Parameters + ---------- + k : :class:`~astropy.units.Quantity` + spectral normalisation + p1 : float + spectral index before the break (positive by definition) + p2 : float + spectral index after the break (positive by definition) + gamma_b : float + Lorentz factor at which the change in spectral index is occurring + gamma_min : float + minimum Lorentz factor of the particle distribution + gamma_max : float + maximum Lorentz factor of the particle distribution + mass : :class:`~astropy.units.Quantity` + particle mass, default is the electron mass + integrator : func + function to be used for integration, default is :class:`~numpy.trapz` + """ + + def __init__( + self, + k=1e-13 * u.Unit("cm-3"), + p1=2.0, + p2=3.0, + gamma_b=1e3, + gamma_min=10, + gamma_max=1e7, + mass=m_e, + integrator=np.trapz, + ): + super().__init__(mass, integrator) + self.k = k + self.p1 = p1 + self.p2 = p2 + self.gamma_b = gamma_b + self.gamma_min = gamma_min + self.gamma_max = gamma_max + + @property + def parameters(self): + return [ + self.k, + self.p1, + self.p2, + self.gamma_b, + self.gamma_min, + self.gamma_max, + ] + + @staticmethod + def evaluate(gamma, k, p1, p2, gamma_b, gamma_min, gamma_max): + index = np.where(gamma <= gamma_b, p1, p2) + return np.where( + (gamma_min <= gamma) * (gamma <= gamma_max), + k * (gamma / gamma_b) ** (-index), + 0, + ) + + def __call__(self, gamma): + return self.evaluate( + gamma, + self.k, + self.p1, + self.p2, + self.gamma_b, + self.gamma_min, + self.gamma_max, + ) + + @staticmethod + def evaluate_SSA_integrand(gamma, k, p1, p2, gamma_b, gamma_min, gamma_max): + r"""Analytical integrand for the synchrotron self-absorption: + :math:`\gamma'^2 \frac{d}{d \gamma'} \left(\frac{n_e(\gamma)}{\gamma'^2}\right)`.""" + index = np.where(gamma <= gamma_b, p1, p2) + return np.where( + (gamma_min <= gamma) * (gamma <= gamma_max), + k * -(index + 2) / gamma * (gamma / gamma_b) ** (-index), + 0, + ) + + def SSA_integrand(self, gamma): + return self.evaluate_SSA_integrand( + gamma, + self.k, + self.p1, + self.p2, + self.gamma_b, + self.gamma_min, + self.gamma_max, + ) + + def __str__(self): + return ( + f"* {self.particle} energy distribution\n" + + f" - broken power law\n" + + f" - k: {self.k:.2e}\n" + + f" - p1: {self.p1:.2f}\n" + + f" - p2: {self.p2:.2f}\n" + + f" - gamma_b: {self.gamma_b:.2e}\n" + + f" - gamma_min: {self.gamma_min:.2e}\n" + + f" - gamma_max: {self.gamma_max:.2e}\n" + ) + + +class LogParabola(ParticleDistribution): + r"""Class describing a log-parabolic particle distribution. + When called, the particle density :math:`n(\gamma)` in :math:`\mathrm{cm}^{-3}` is returned. + + .. math:: + n(\gamma') = k \, \left(\frac{\gamma'}{\gamma'_0}\right)^{-(p + q \log_{10}(\gamma' / \gamma'_0))} + + Parameters + ---------- + k : :class:`~astropy.units.Quantity` + spectral normalisation + p : float + spectral index, note it is positive by definition, will change sign in the function + q : float + spectral curvature, note it is positive by definition, will change sign in the function + gamma_0 : float + reference Lorentz factor + gamma_min : float + minimum Lorentz factor of the particle distribution + gamma_max : float + maximum Lorentz factor of the particle distribution + mass : :class:`~astropy.units.Quantity` + particle mass, default is the electron mass + integrator : func + function to be used for integration, default is :class:`~numpy.trapz` + """ + + def __init__( + self, + k=1e-13 * u.Unit("cm-3"), + p=2.0, + q=0.1, + gamma_0=1e3, + gamma_min=10, + gamma_max=1e7, + mass=m_e, + integrator=np.trapz, + ): + super().__init__(mass, integrator) + self.k = k + self.p = p + self.q = q + self.gamma_0 = gamma_0 + self.gamma_min = gamma_min + self.gamma_max = gamma_max + + @property + def parameters(self): + return [self.k, self.p, self.q, self.gamma_0, self.gamma_min, self.gamma_max] + + @staticmethod + def evaluate(gamma, k, p, q, gamma_0, gamma_min, gamma_max): + gamma_ratio = gamma / gamma_0 + index = -p - q * np.log10(gamma_ratio) + return np.where( + (gamma_min <= gamma) * (gamma <= gamma_max), k * gamma_ratio ** index, 0, + ) + + def __call__(self, gamma): + return self.evaluate( + gamma, self.k, self.p, self.q, self.gamma_0, self.gamma_min, self.gamma_max, + ) + + @staticmethod + def evaluate_SSA_integrand(gamma, k, p, q, gamma_0, gamma_min, gamma_max): + r"""Analytical integrand for the synchrotron self-absorption: + :math:`\gamma'^2 \frac{d}{d \gamma'} \left(\frac{n_e(\gamma)}{\gamma'^2}\right)`.""" + prefactor = -(p + 2 * q * np.log10(gamma / gamma_0) + 2) / gamma + return prefactor * LogParabola.evaluate( + gamma, k, p, q, gamma_0, gamma_min, gamma_max + ) + + def SSA_integrand(self, gamma): + return self.evaluate_SSA_integrand( + gamma, self.k, self.p, self.q, self.gamma_0, self.gamma_min, self.gamma_max, + ) + + def __str__(self): + return ( + f"* {self.particle} energy distribution\n" + + f" - log parabola\n" + + f" - k: {self.k:.2e}\n" + + f" - p: {self.p:.2f}\n" + + f" - q: {self.q:.2f}\n" + + f" - gamma_0: {self.gamma_0:.2e}\n" + + f" - gamma_min: {self.gamma_min:.2e}\n" + + f" - gamma_max: {self.gamma_max:.2e}\n" + ) + + +class ExpCutoffPowerLaw(ParticleDistribution): + r"""Class describing a power-law with an exponetial cutoff particle distribution. + When called, the particle density :math:`n_e(\gamma)` in :math:`\mathrm{cm}^{-3}` is returned. + + .. math:: + n(\gamma'_c) = k \, \gamma'^{-p} exp(-\gamma'/\gamma_c) \, H(\gamma'; \gamma'{\rm min}, \gamma'{\rm max}) + + Parameters + ---------- + k : :class:`~astropy.units.Quantity` + spectral normalisation + p : float + spectral index, note it is positive by definition, will change sign in the function + gamma_c : float + cutoff Lorentz factor of the particle distribution + gamma_min : float + minimum Lorentz factor of the particle distribution + gamma_max : float + maximum Lorentz factor of the particle distribution + mass : :class:`~astropy.units.Quantity` + particle mass, default is the electron mass + integrator : func + function to be used for integration, default is :class:`~numpy.trapz` + """ + + def __init__( + self, + k=1e-13 * u.Unit("cm-3"), + p=2.1, + gamma_c=1e3, + gamma_min=10, + gamma_max=1e5, + mass=m_e, + integrator=np.trapz, + ): + super().__init__(mass, integrator) + self.k = k + self.p = p + self.gamma_c = gamma_c + self.gamma_min = gamma_min + self.gamma_max = gamma_max + + @property + def parameters(self): + return [self.k, self.p, self.gamma_c, self.gamma_min, self.gamma_max] + + @staticmethod + def evaluate(gamma, k, p, gamma_c, gamma_min, gamma_max): + return np.where( + (gamma_min <= gamma) * (gamma <= gamma_max), + k * gamma ** (-p) * np.exp(-gamma / gamma_c), + 0, + ) + + def __call__(self, gamma): + return self.evaluate( + gamma, self.k, self.p, self.gamma_c, self.gamma_min, self.gamma_max + ) + + @staticmethod + def evaluate_SSA_integrand(gamma, k, p, gamma_c, gamma_min, gamma_max): + r"""(analytical) integrand for the synchrotron self-absorption: + :math:`\gamma'^2 \frac{d}{d \gamma'} \left(\frac{n(\gamma)}{\gamma'^2}\right)`""" + prefactor = -(p + 2) / gamma + (-1 / gamma_c) + + return prefactor * ExpCutoffPowerLaw.evaluate( + gamma, k, p, gamma_c, gamma_min, gamma_max + ) + + def SSA_integrand(self, gamma): + return self.evaluate_SSA_integrand( + gamma, self.k, self.p, self.gamma_c, self.gamma_min, self.gamma_max + ) + + def __str__(self): + return ( + f"* {self.particle} energy distribution\n" + + f" - exponential cut-off power law\n" + + f" - k: {self.k:.2e}\n" + + f" - p: {self.p:.2f}\n" + + f" - gamma_c: {self.gamma_c:.2f}\n" + + f" - gamma_min: {self.gamma_min:.2e}\n" + + f" - gamma_max: {self.gamma_max:.2e}\n" + ) + + +class ExpCutoffBrokenPowerLaw(ParticleDistribution): + r"""Class describing an exponential cutoff broken power-law particle distribution. + When called, the particle density :math:`n(\gamma)` in :math:`\mathrm{cm}^{-3}` is returned. + + .. math:: + n(\gamma'_c) = k \left[ + \left(\frac{\gamma'}{\gamma'_b}\right)^{-p_1} exp(-\gamma'/\gamma_c) \, H(\gamma'; \gamma'_{\rm min}, \gamma'_b) + + \left(\frac{\gamma'}{\gamma'_b}\right)^{-p_2} exp(-\gamma'/\gamma_c)\, H(\gamma'; \gamma'_{b}, \gamma'_{\rm max}) + \right] + + Parameters + ---------- + k : :class:`~astropy.units.Quantity` + spectral normalisation + p1 : float + spectral index before the break (positive by definition) + p2 : float + spectral index after the break (positive by definition) + gamma_b : float + Lorentz factor at which the change in spectral index is occurring + gamma_min : float + minimum Lorentz factor of the particle distribution + gamma_max : float + maximum Lorentz factor of the particle distribution + gamma_c : float + cutoff Lorentz factor of the particle distribution + mass : `~astropy.units.Quantity` + particle mass, default is the electron mass + integrator : func + function to be used for integration, default is :class:`~numpy.trapz` + """ + + def __init__( + self, + k=1e-13 * u.Unit("cm-3"), + p1=2.0, + p2=3.0, + gamma_c = 1e5, + gamma_b=1e3, + gamma_min=10, + gamma_max=1e7, + mass=m_e, + integrator=np.trapz, + ): + super().__init__(mass, integrator) + self.k = k + self.p1 = p1 + self.p2 = p2 + self.gamma_c = gamma_c + self.gamma_b = gamma_b + self.gamma_min = gamma_min + self.gamma_max = gamma_max + + @property + def parameters(self): + return [ + self.k, + self.p1, + self.p2, + self.gamma_c, + self.gamma_b, + self.gamma_min, + self.gamma_max, + ] + + @staticmethod + def evaluate(gamma, k, p1, p2, gamma_c, gamma_b, gamma_min, gamma_max): + index = np.where(gamma <= gamma_b, p1, p2) + return np.where( + (gamma_min <= gamma) * (gamma <= gamma_max), + k * (gamma/gamma_b) ** (-index) * np.exp(-gamma/gamma_c), + 0 + ) + + def __call__(self, gamma): + return self.evaluate( + gamma, + self.k, + self.p1, + self.p2, + self.gamma_c, + self.gamma_b, + self.gamma_min, + self.gamma_max, + ) + + @staticmethod + def evaluate_SSA_integrand(gamma, k, p1, p2, gamma_c, gamma_b, gamma_min, gamma_max): + r"""Analytical integrand for the synchrotron self-absorption: + :math:`\gamma'^2 \frac{d}{d \gamma'} \left(\frac{n_e(\gamma)}{\gamma'^2}\right)`.""" + index = np.where(gamma <= gamma_b, p1, p2) + prefactor = -(index + 2) / gamma + (-1 / gamma_c) + return np.where( + (gamma_min <= gamma) * (gamma <= gamma_max), + prefactor * ExpCutoffBrokenPowerLaw.evaluate( + gamma, k, p1, p2, gamma_c, gamma_b, gamma_min, gamma_max), + 0 + ) + + def SSA_integrand(self, gamma): + return self.evaluate_SSA_integrand( + gamma, + self.k, + self.p1, + self.p2, + self.gamma_c, + self.gamma_b, + self.gamma_min, + self.gamma_max + ) + + def __str__(self): + return ( + f"* {self.particle} energy distribution\n" + + f" - exponential cut-off broken power law\n" + + f" - k: {self.k:.2e}\n" + + f" - p1: {self.p1:.2f}\n" + + f" - p2: {self.p2:.2f}\n" + + f" - gamma_c: {self.gamma_c:.2e}\n" + + f" - gamma_b: {self.gamma_b:.2e}\n" + + f" - gamma_min: {self.gamma_min:.2e}\n" + + f" - gamma_max: {self.gamma_max:.2e}\n" + ) + + +class InterpolatedDistribution(ParticleDistribution): + """Class describing a particle distribution with an arbitrary shape. + The spectrum is interpolated from an array of Lorentz factor and densities. + + Parameters + ---------- + gamma : :class:`~numpy.ndarray` + array of Lorentz factors where the density has been computed + n : :class:`~astropy.units.Quantity` + array of densities to be interpolated + norm : float + parameter to scale the density + mass : :class:`~astropy.units.Quantity` + particle mass, default is the electron mass + integrator : func + function to be used for integration, default is :class:`~numpy.trapz` + """ + + def __init__(self, gamma, n, norm=1, mass=m_e, integrator=np.trapz): + super().__init__(mass, integrator) + self.gamma = gamma + self.gamma_max = np.max(self.gamma) + self.gamma_min = np.min(self.gamma) + self.norm = norm + if n.unit != u.Unit("cm-3"): + raise ValueError( + f"Provide a particle distribution in cm-3, instead of {n.unit}" + ) + else: + self.n = n + # call make the interpolation + self.log10_f = self.log10_interpolation() + + def log10_interpolation(self): + """Returns the function interpolating in log10 the particle spectrum. + TODO: make possible to pass arguments to CubicSpline. + """ + interpolator = CubicSpline( + np.log10(self.gamma), np.log10(self.n.to_value("cm-3")) + ) + return interpolator + + def evaluate(self, gamma, norm, gamma_min, gamma_max): + log10_gamma = np.log10(gamma) + values = np.where( + (gamma_min <= gamma) * (gamma <= gamma_max), + norm * np.power(10, self.log10_f(log10_gamma)), + 0, + ) + return values * u.Unit("cm-3") + + def __call__(self, gamma): + return self.evaluate(gamma, self.norm, self.gamma_min, self.gamma_max) + + def SSA_integrand(self, gamma): + r"""Integrand for the synchrotron self-absorption. It is + + .. math:: + \gamma^2 \frac{d}{d \gamma} (\frac{n_e(\gamma)}{\gamma^2}) = ( \frac{dn_e(\gamma)}{d\gamma}+\frac{2n_e(\gamma)}{\gamma}) + + The derivative is: + + .. math:: + \frac{dn_e(\gamma)}{d\gamma} = \frac{d 10^{f(u(\gamma))}}{d\gamma} = \frac{d10^{f(u)}}{du} \cdot \frac{du(\gamma)}{d\gamma} + + where we have :math:`\frac{d 10^{f(u(\gamma))}}{d\gamma} = \frac{d10^{f(u)}}{du} \cdot \frac{du(\gamma)}{d\gamma}`, + where :math:`u` is the :math:`log_{10}(\gamma)`. + This is equal to :math:`\frac{d 10^{f(u(\gamma))}}{d\gamma} = 10^{f(u)} \cdot \frac{df(u)}{du} \cdot \frac{1}{\gamma}` + + """ + log10_gamma = np.log10(gamma) + df_log = self.log10_f.derivative() + int_fun = self.evaluate(gamma, self.norm, self.gamma_min, self.gamma_max) + deriv = int_fun * (1 / gamma) * df_log(log10_gamma) + return deriv - 2 * int_fun / gamma diff --git a/agnpy/spectra/physical.py b/agnpy/spectra/physical.py new file mode 100644 index 00000000..ef8382cf --- /dev/null +++ b/agnpy/spectra/physical.py @@ -0,0 +1,84 @@ +# physical particle energy distributions +import numpy as np +import astropy.units as u +from astropy.constants import m_e, c, sigma_T +from sympy import Symbol, uppergamma, lambdify +from .spectra import ParticleDistribution +from ..utils.conversion import B_to_cgs + + +# lambdify the uppergamma function of sympy (make it applicable to numpy arrays) +# https://docs.sympy.org/latest/modules/functions/special.html#sympy.functions.special.gamma_functions.uppergamma +s = Symbol("s") +t = Symbol("t") +expr = uppergamma(s, t) +lambda_uppergamma = lambdify((s, t), expr, modules=["numpy", "sympy"]) +vecotrised_uppergamma = np.vectorize(lambda_uppergamma) + + +def gamma_b_tau_synchrotron_tau_ad(B, R): + """Break Lorentz factor assuming equal time scales in cooling and adiabatic expansion. + Eq. 2.26 in [Inoue]_""" + B = B_to_cgs(B) + u_B = B**2 / (4 * np.pi) + return (3 * m_e * c**2 / (4 * sigma_T * u_B * R)).to_value("") + + +class SteadyStateSynchrotronCooling(ParticleDistribution): + r"""Solution of the steady state electron differential equation + + .. math:: + N(\gamma') = ... + + Parameters + ---------- + k : :class:`~astropy.units.Quantity` + injection rate + p : float + spectral index of the injected power law + B : float + magnetic field + tau_ad : :class:`~astropy.units.Quantity` + adiabatic time scale + """ + def __init__( + self, + q_inj=1e-5 * u.Unit("cm-3 s-1"), + p=2, + B=0.1 * u.G, + R_b=1e17 * u.cm, + mass=m_e, + integrator=np.trapz, + ): + super().__init__(mass, integrator) + self.q_inj = q_inj + self.p = p + self.B = B + self.R_b = R_b + + @property + def parameters(self): + return [self.q_inj, self.p, self.B, self.R_b] + + @staticmethod + def evaluate(gamma, q_inj, p, B, R_b): + """Solution, Eq. 2.26 in [Inoue]""" + tau_ad = R_b / c + gamma_b = gamma_b_tau_synchrotron_tau_ad(B, R_b) + print(gamma_b) + prefactor = ( + np.exp(gamma / gamma_b) + * tau_ad.to("s") + * q_inj.to("cm-3 s-1") + * gamma_b ** (2 - p) + / gamma**2 + ) + return prefactor * vecotrised_uppergamma(1 - p, gamma / gamma_b) + + def __call__(self, gamma): + return self.evaluate(gamma, self.q_inj, self.p, self.B, self.R_b) + + def __str__(self): + return ( + f"ciao" + ) From 1492e8c188adbb9ab098a94e89e8bfd45dc0476a Mon Sep 17 00:00:00 2001 From: cosimoNigro Date: Tue, 16 May 2023 10:58:19 +0200 Subject: [PATCH 3/3] modified blob and synchrotron classes to use physical solutions --- agnpy/emission_regions/blob.py | 20 ++++++++++++++++---- agnpy/synchrotron/synchrotron.py | 4 ++-- 2 files changed, 18 insertions(+), 6 deletions(-) diff --git a/agnpy/emission_regions/blob.py b/agnpy/emission_regions/blob.py index 246d649e..2a46f75c 100644 --- a/agnpy/emission_regions/blob.py +++ b/agnpy/emission_regions/blob.py @@ -5,7 +5,7 @@ import astropy.units as u from astropy.coordinates import Distance from astropy.constants import c, sigma_T, m_e -from ..spectra import PowerLaw +from ..spectra import PowerLaw, SteadyStateSynchrotronCooling from ..utils.conversion import mec2, mpc2, B_to_cgs @@ -63,17 +63,29 @@ def __init__( gamma_e_size=200, gamma_p_size=200, ): - self.R_b = R_b.to("cm") + # if this is a physical solutions, then we want to inherit its B and R parameters + if isinstance(n_e, SteadyStateSynchrotronCooling): + self.R_b = n_e.R_b.to("cm") + self.B = n_e.B + else: + self.R_b = R_b.to("cm") + self.B = B + self.z = z self.delta_D = delta_D self.Gamma = Gamma - self.B = B self._n_e = n_e self._n_p = n_p self.xi = xi + # set the parameters of the array of electrons Lorentz factors + # physical solutions don't have gamma_min and gamma_max + if isinstance(n_e, SteadyStateSynchrotronCooling): + self.set_gamma_e(gamma_e_size, 10, 1e6) + else: + self.set_gamma_e(gamma_e_size, self._n_e.gamma_min, self._n_e.gamma_max) + # we might want to have different array of Lorentz factors for e and p - self.set_gamma_e(gamma_e_size, self._n_e.gamma_min, self._n_e.gamma_max) if self._n_p is not None: self.set_gamma_p(gamma_p_size, self._n_p.gamma_min, self._n_p.gamma_max) diff --git a/agnpy/synchrotron/synchrotron.py b/agnpy/synchrotron/synchrotron.py index 3b4c85da..33980329 100644 --- a/agnpy/synchrotron/synchrotron.py +++ b/agnpy/synchrotron/synchrotron.py @@ -6,7 +6,7 @@ from ..utils.conversion import nu_to_epsilon_prime, B_to_cgs, lambda_c_e -__all__ = ["R", "nu_synch_peak", "Synchrotron"] +__all__ = ["R", "nu_synch_peak", "single_particle_synch_power", "Synchrotron"] e = e.gauss B_cr = 4.414e13 * u.G # critical magnetic field @@ -178,7 +178,7 @@ def evaluate_sed_flux( V_b = 4 / 3 * np.pi * np.power(R_b, 3) N_e = V_b * n_e.evaluate(_gamma, *args) # fold the electron distribution with the synchrotron power - integrand = N_e * single_electron_synch_power(B_cgs, _epsilon, _gamma) + integrand = N_e * single_particle_synch_power(B_cgs, _epsilon, _gamma) emissivity = integrator(integrand, gamma, axis=0).reshape(epsilon.shape) prefactor = np.power(delta_D, 4) / (4 * np.pi * np.power(d_L, 2)) sed = (prefactor * epsilon * emissivity).to("erg cm-2 s-1")