Source code for gammapy.astro.darkmatter.utils

# Licensed under a 3-clause BSD style license - see LICENSE.rst
"""Utilities to compute J-factor maps."""
import html
import numpy as np
import astropy.units as u

__all__ = ["JFactory"]


[docs]class JFactory: """Compute J-Factor or D-Factor maps. J-Factors are computed for annihilation and D-Factors for decay. Set the argument `annihilation` to `False` to compute D-Factors. The assumed dark matter profiles will be centered on the center of the map. Parameters ---------- geom : `~gammapy.maps.WcsGeom` Reference geometry. profile : `~gammapy.astro.darkmatter.profiles.DMProfile` Dark matter profile. distance : `~astropy.units.Quantity` Distance to convert angular scale of the map. annihilation: bool, optional Decay or annihilation. Default is True. """ def __init__(self, geom, profile, distance, annihilation=True): self.geom = geom self.profile = profile self.distance = distance self.annihilation = annihilation def _repr_html_(self): try: return self.to_html() except AttributeError: return f"<pre>{html.escape(str(self))}</pre>"
[docs] def compute_differential_jfactor(self, ndecade=1e4): r"""Compute differential J-Factor. .. math:: \frac{\mathrm d J_\text{ann}}{\mathrm d \Omega} = \int_{\mathrm{LoS}} \mathrm d l \rho(l)^2 .. math:: \frac{\mathrm d J_\text{decay}}{\mathrm d \Omega} = \int_{\mathrm{LoS}} \mathrm d l \rho(l) """ separation = self.geom.separation(self.geom.center_skydir).rad rmin = u.Quantity( value=np.tan(separation) * self.distance, unit=self.distance.unit ) rmax = self.distance val = [ ( 2 * self.profile.integral( _.value * u.kpc, rmax, np.arctan(_.value / self.distance.value), ndecade, self.annihilation, ) + self.profile.integral( self.distance, 4 * rmax, np.arctan(_.value / self.distance.value), ndecade, self.annihilation, ) ) for _ in rmin.ravel() ] integral_unit = u.Unit("GeV2 cm-5") if self.annihilation else u.Unit("GeV cm-2") jfact = u.Quantity(val).to(integral_unit).reshape(rmin.shape) return jfact / u.steradian
[docs] def compute_jfactor(self, ndecade=1e4): r"""Compute astrophysical J-Factor. .. math:: J(\Delta\Omega) = \int_{\Delta\Omega} \mathrm d \Omega^{\prime} \frac{\mathrm d J}{\mathrm d \Omega^{\prime}} """ diff_jfact = self.compute_differential_jfactor(ndecade) return diff_jfact * self.geom.to_image().solid_angle()