# Licensed under a 3-clause BSD style license - see LICENSE.rst"""Utilities to compute J-factor maps."""importhtmlimportnumpyasnpimportastropy.unitsasu__all__=["JFactory"]
[docs]classJFactory:"""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=geomself.profile=profileself.distance=distanceself.annihilation=annihilationdef_repr_html_(self):try:returnself.to_html()exceptAttributeError:returnf"<pre>{html.escape(str(self))}</pre>"
[docs]defcompute_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).radrmin=u.Quantity(value=np.tan(separation)*self.distance,unit=self.distance.unit)rmax=self.distanceval=[(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_inrmin.ravel()]integral_unit=u.Unit("GeV2 cm-5")ifself.annihilationelseu.Unit("GeV cm-2")jfact=u.Quantity(val).to(integral_unit).reshape(rmin.shape)returnjfact/u.steradian
[docs]defcompute_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)returndiff_jfact*self.geom.to_image().solid_angle()