Source code for gammapy.astro.population.spatial

# Licensed under a 3-clause BSD style license - see LICENSE.rst
"""Galactic radial source distribution probability density functions.

Attributes
----------
radial_distributions : `~collections.OrderedDict`
    Dictionary of available spatial distributions.

    Useful for automatic processing.
"""
from __future__ import absolute_import, division, print_function, unicode_literals
from collections import OrderedDict
import numpy as np
from numpy import exp, pi, log, abs, cos, sin
from astropy.units import Quantity
from astropy.modeling import Fittable1DModel, Parameter
from ...utils.coordinates import cartesian, polar, D_SUN_TO_GALACTIC_CENTER
from ...utils.random import get_random_state

__all__ = [
    'CaseBattacharya1998',
    'FaucherKaspi2006',
    'Lorimer2006',
    'Paczynski1990',
    'YusifovKucuk2004',
    'YusifovKucuk2004B',
    'Exponential',
    'LogSpiral',
    'FaucherSpiral',
    'ValleeSpiral',
    'radial_distributions',
    'RMIN', 'RMAX',
    'ZMIN', 'ZMAX',
]

# Simulation range used for random number drawing
RMIN, RMAX = Quantity([0, 20], 'kpc')
ZMIN, ZMAX = Quantity([-0.5, 0.5], 'kpc')


[docs]class Paczynski1990(Fittable1DModel): """Radial distribution of the birth surface density of neutron stars - Paczynski 1990. .. math :: f(r) = A r_{exp}^{-2} \\exp \\left(-\\frac{r}{r_{exp}} \\right) Reference: http://adsabs.harvard.edu/abs/1990ApJ...348..485P (Formula (2)) Parameters ---------- amplitude : float See formula r_exp : float See formula See Also -------- CaseBattacharya1998, YusifovKucuk2004, Lorimer2006, YusifovKucuk2004B, FaucherKaspi2006, Exponential """ amplitude = Parameter() r_exp = Parameter() evolved = False def __init__(self, amplitude=1, r_exp=4.5, **kwargs): super(Paczynski1990, self).__init__(amplitude=amplitude, r_exp=r_exp, **kwargs) @staticmethod
[docs] def evaluate(r, amplitude, r_exp): """One dimensional Paczynski 1990 model function""" return amplitude * r_exp ** -2 * np.exp(-r / r_exp)
[docs]class CaseBattacharya1998(Fittable1DModel): """Radial distribution of the surface density of supernova remnants in the galaxy - Case & Battacharya 1998. .. math :: f(r) = A \\left( \\frac{r}{r_{\\odot}} \\right) ^ \\alpha \\exp \\left[ -\\beta \\left( \\frac{ r - r_{\\odot}}{r_{\\odot}} \\right) \\right] Reference: http://adsabs.harvard.edu//abs/1998ApJ...504..761C (Formula (14)) Parameters ---------- amplitude : float See model formula alpha : float See model formula beta : float See model formula See Also -------- Paczynski1990, YusifovKucuk2004, Lorimer2006, YusifovKucuk2004B, FaucherKaspi2006, Exponential """ amplitude = Parameter() alpha = Parameter() beta = Parameter() evolved = True def __init__(self, amplitude=1., alpha=2, beta=3.53, **kwargs): super(CaseBattacharya1998, self).__init__(amplitude=amplitude, alpha=alpha, beta=beta, **kwargs) @staticmethod
[docs] def evaluate(r, amplitude, alpha, beta): """One dimensional Case & Battacharya 2006 model function""" term1 = (r / D_SUN_TO_GALACTIC_CENTER.value) ** alpha term2 = np.exp(-beta * (r - D_SUN_TO_GALACTIC_CENTER.value) / D_SUN_TO_GALACTIC_CENTER.value) return amplitude * term1 * term2
[docs]class YusifovKucuk2004(Fittable1DModel): """Radial distribution of the surface density of pulsars in the galaxy - Yusifov & Kucuk 2004. .. math :: f(r) = A \\left ( \\frac{r + r_1}{r_{\\odot} + r_1} \\right )^a \\exp \\left [-b \\left( \\frac{r - r_{\\odot}}{r_{\\odot} + r_1} \\right ) \\right ] Used by Faucher-Guigere and Kaspi. Density at ``r = 0`` is nonzero. Reference: http://adsabs.harvard.edu/abs/2004A%26A...422..545Y (Formula (15)) Parameters ---------- amplitude : float See model formula a : float See model formula b : float See model formula r_1 : float See model formula See Also -------- CaseBattacharya1998, Paczynski1990, Lorimer2006, YusifovKucuk2004B, FaucherKaspi2006, Exponential """ amplitude = Parameter() a = Parameter() b = Parameter() r_1 = Parameter() evolved = True def __init__(self, amplitude=1, a=1.64, b=4.01, r_1=0.55, **kwargs): super(YusifovKucuk2004, self).__init__(amplitude=amplitude, a=a, b=b, r_1=r_1, **kwargs) @staticmethod
[docs] def evaluate(r, amplitude, a, b, r_1): """One dimensional Yusifov & Kucuk 2004 model function""" term1 = ((r + r_1) / (D_SUN_TO_GALACTIC_CENTER.value + r_1)) ** a term2 = np.exp(-b * (r - D_SUN_TO_GALACTIC_CENTER.value) / (D_SUN_TO_GALACTIC_CENTER.value + r_1)) return amplitude * term1 * term2
[docs]class YusifovKucuk2004B(Fittable1DModel): """Radial distribution of the surface density of OB stars in the galaxy - Yusifov & Kucuk 2004. .. math :: f(r) = A \\left( \\frac{r}{r_{\\odot}} \\right) ^ a \\exp \\left[ -b \\left( \\frac{r}{r_{\\odot}} \\right) \\right] Derived empirically from OB-stars distribution. Reference: http://adsabs.harvard.edu/abs/2004A%26A...422..545Y (Formula (17)) Parameters ---------- amplitude : float See model formula a : float See model formula b : float See model formula See Also -------- CaseBattacharya1998, Paczynski1990, YusifovKucuk2004, Lorimer2006, FaucherKaspi2006, Exponential """ amplitude = Parameter() a = Parameter() b = Parameter() evolved = False def __init__(self, amplitude=1, a=4, b=6.8, **kwargs): super(YusifovKucuk2004B, self).__init__(amplitude=amplitude, a=a, b=b, **kwargs) @staticmethod
[docs] def evaluate(r, amplitude, a, b): """One dimensional Yusifov & Kucuk 2004 model function""" return amplitude * (r / D_SUN_TO_GALACTIC_CENTER.value) ** a * np.exp(-b * (r / D_SUN_TO_GALACTIC_CENTER.value))
[docs]class FaucherKaspi2006(Fittable1DModel): """Radial distribution of the birth surface density of pulsars in the galaxy - Faucher-Giguere & Kaspi 2006. .. math :: f(r) = A \\frac{1}{\\sqrt{2 \pi} \sigma} \\exp \\left(- \\frac{(r - r_0)^2}{2 \sigma ^ 2}\\right) Reference: http://adsabs.harvard.edu/abs/2006ApJ...643..332F (Appendix B) Parameters ---------- amplitude : float See model formula r_0 : float See model formula sigma : float See model formula See Also -------- CaseBattacharya1998, Paczynski1990, YusifovKucuk2004, Lorimer2006, YusifovKucuk2004B, Exponential """ amplitude = Parameter() r_0 = Parameter() sigma = Parameter() evolved = False def __init__(self, amplitude=1, r_0=7.04, sigma=1.83, **kwargs): super(FaucherKaspi2006, self).__init__(amplitude=amplitude, r_0=r_0, sigma=sigma, **kwargs) @staticmethod
[docs] def evaluate(r, amplitude, r_0, sigma): """One dimensional Faucher-Giguere & Kaspi 2006 model function""" term1 = 1. / np.sqrt(2 * pi * sigma) term2 = np.exp(-(r - r_0) ** 2 / (2 * sigma ** 2)) return amplitude * term1 * term2
[docs]class Lorimer2006(Fittable1DModel): """Radial distribution of the surface density of pulsars in the galaxy - Lorimer 2006. .. math :: f(r) = A \\left( \\frac{r}{r_{\\odot}} \\right) ^ B \\exp \\left[ -C \\left( \\frac{r - r_{\\odot}}{r_{\\odot}} \\right) \\right] Reference: http://adsabs.harvard.edu/abs/2006MNRAS.372..777L (Formula (10)) Parameters ---------- amplitude : float See model formula B : float See model formula C : float See model formula See Also -------- CaseBattacharya1998, Paczynski1990, YusifovKucuk2004, Lorimer2006, YusifovKucuk2004B, FaucherKaspi2006 """ amplitude = Parameter() B = Parameter() C = Parameter() evolved = True def __init__(self, amplitude=1, B=1.9, C=5.0, **kwargs): super(Lorimer2006, self).__init__(amplitude=amplitude, B=B, C=C, **kwargs) @staticmethod
[docs] def evaluate(r, amplitude, B, C): """One dimensional Lorimer 2006 model function""" term1 = (r / D_SUN_TO_GALACTIC_CENTER.value) ** B term2 = np.exp(-C * (r - D_SUN_TO_GALACTIC_CENTER.value) / D_SUN_TO_GALACTIC_CENTER.value) return amplitude * term1 * term2
[docs]class Exponential(Fittable1DModel): """Exponential distribution. .. math :: f(z) = A \\exp \\left(- \\frac{|z|}{z_0} \\right) Usually used for height distribution above the Galactic plane, with 0.05 kpc as a commonly used birth height distribution. Parameters ---------- amplitude : float See model formula z_0 : float Scale height of the distribution See Also -------- CaseBattacharya1998, Paczynski1990, YusifovKucuk2004, Lorimer2006, YusifovKucuk2004B, FaucherKaspi2006, Exponential """ amplitude = Parameter() z_0 = Parameter() evolved = False def __init__(self, amplitude=1, z_0=0.05, **kwargs): super(Exponential, self).__init__(amplitude=amplitude, z_0=z_0, **kwargs) @staticmethod
[docs] def evaluate(z, amplitude, z_0): """One dimensional exponential model function""" return amplitude * np.exp(-np.abs(z) / z_0)
[docs]class LogSpiral(object): """Logarithmic spiral. Reference: http://en.wikipedia.org/wiki/Logarithmic_spiral """
[docs] def xy_position(self, theta=None, radius=None, spiralarm_index=0): """Compute (x, y) position for a given angle or radius. Parameters ---------- theta : array_like Angle (deg) radius : array_like Radius (kpc) spiralarm_index : int Spiral arm index Returns ------- x, y : array_like Position (x, y) """ if (theta is None) and not (radius is None): theta = self.theta(radius, spiralarm_index=spiralarm_index) elif (radius is None) and not (theta is None): radius = self.radius(theta, spiralarm_index=spiralarm_index) else: ValueError('Specify only one of: theta, radius') theta = np.radians(theta) x = radius * cos(theta) y = radius * sin(theta) return x, y
[docs] def radius(self, theta, spiralarm_index): """Radius for a given angle. Parameters ---------- theta : array_like Angle (deg) spiralarm_index : int Spiral arm index Returns ------- radius : array_like Radius (kpc) """ k = self.k[spiralarm_index] r_0 = self.r_0[spiralarm_index] theta_0 = self.theta_0[spiralarm_index] d_theta = np.radians(theta - theta_0) radius = r_0 * exp(d_theta / k) return radius
[docs] def theta(self, radius, spiralarm_index): """Angle for a given radius. Parameters ---------- radius : array_like Radius (kpc) spiralarm_index : int Spiral arm index Returns ------- theta : array_like Angle (deg) """ k = self.k[spiralarm_index] r_0 = self.r_0[spiralarm_index] theta_0 = self.theta_0[spiralarm_index] theta_0 = np.radians(theta_0) theta = k * log(radius / r_0) + theta_0 return np.degrees(theta)
[docs]class FaucherSpiral(LogSpiral): """Milky way spiral arm used in Faucher et al (2006). Reference: http://adsabs.harvard.edu/abs/2006ApJ...643..332F """ # Parameters k = Quantity([4.25, 4.25, 4.89, 4.89], 'rad') r_0 = Quantity([3.48, 3.48, 4.9, 4.9], 'kpc') theta_0 = Quantity([1.57, 4.71, 4.09, 0.95], 'rad') spiralarms = np.array(['Norma', 'Carina Sagittarius', 'Perseus', 'Crux Scutum']) def _blur(self, radius, theta, amount=0.07, random_state='random-seed'): """Blur the positions around the centroid of the spiralarm. The given positions are blurred by drawing a displacement in radius from a normal distribution, with sigma = amount * radius. And a direction theta from a uniform distribution in the interval [0, 2 * pi]. Parameters ---------- radius : `~astropy.units.Quantity` Radius coordinate theta : `~astropy.units.Quantity` Angle coordinate amount: float, optional Amount of blurring of the position, given as a fraction of `radius`. random_state : {int, 'random-seed', 'global-rng', `~numpy.random.RandomState`} Defines random number generator initialisation. Passed to `~gammapy.utils.random.get_random_state`. """ random_state = get_random_state(random_state) dr = Quantity(abs(random_state.normal(0, amount * radius, radius.size)), 'kpc') dtheta = Quantity(random_state.uniform(0, 2 * np.pi, radius.size), 'rad') x, y = cartesian(radius, theta) dx, dy = cartesian(dr, dtheta) return polar(x + dx, y + dy) def _gc_correction(self, radius, theta, r_corr=Quantity(2.857, 'kpc'), random_state='random-seed'): """Correction of source distribution towards the galactic center. To avoid spiralarm features near the Galactic Center, the position angle theta is blurred by a certain amount towards the GC. Parameters ---------- radius : `~astropy.units.Quantity` Radius coordinate theta : `~astropy.units.Quantity` Angle coordinate r_corr : `~astropy.units.Quantity`, optional Scale of the correction towards the GC random_state : {int, 'random-seed', 'global-rng', `~numpy.random.RandomState`} Defines random number generator initialisation. Passed to `~gammapy.utils.random.get_random_state`. """ random_state = get_random_state(random_state) theta_corr = Quantity(random_state.uniform(0, 2 * pi, radius.size), 'rad') return radius, theta + theta_corr * np.exp(-radius / r_corr)
[docs] def __call__(self, radius, blur=True, random_state='random-seed'): """Draw random position from spiral arm distribution. Returns the corresponding angle theta[rad] to a given radius[kpc] and number of spiralarm. Possible numbers are: * Norma = 0, * Carina Sagittarius = 1, * Perseus = 2 * Crux Scutum = 3. Parameters ---------- random_state : {int, 'random-seed', 'global-rng', `~numpy.random.RandomState`} Defines random number generator initialisation. Passed to `~gammapy.utils.random.get_random_state`. Returns ------- Returns dx and dy, if blurring= true. """ random_state = get_random_state(random_state) N = random_state.randint(0, 4, radius.size) # Choose Spiralarm theta = self.k[N] * log(radius / self.r_0[N]) + self.theta_0[N] # Compute angle spiralarm = self.spiralarms[N] # List that contains in wich spiralarm a postion lies if blur: # Apply blurring model according to Faucher radius, theta = self._blur(radius, theta, random_state=random_state) radius, theta = self._gc_correction(radius, theta, random_state=random_state) return radius, theta, spiralarm
[docs]class ValleeSpiral(LogSpiral): """Milky way spiral arm model from Vallee (2008). Reference: http://adsabs.harvard.edu/abs/2008AJ....135.1301V """ # Model parameters p = Quantity(12.8, 'deg') # pitch angle in deg m = 4 # number of spiral arms r_sun = Quantity(7.6, 'kpc') # distance sun to Galactic center in kpc r_0 = Quantity(2.1, 'kpc') # spiral inner radius in kpc theta_0 = Quantity(-20, 'deg') # Norma spiral arm start angle bar_radius = Quantity(3.0, 'kpc') # Radius of the galactic bar (not equal r_0!) spiralarms = np.array(['Norma', 'Perseus', 'Carina Sagittarius', 'Crux Scutum']) def __init__(self): self.r_0 = self.r_0 * np.ones(4) self.theta_0 = self.theta_0 + Quantity([0, 90, 180, 270], 'deg') self.k = Quantity(1. / np.tan(np.radians(self.p.value)) * np.ones(4), 'rad') # Compute start and end point of the bar x_0, y_0 = self.xy_position(radius=self.bar_radius, spiralarm_index=0) x_1, y_1 = self.xy_position(radius=self.bar_radius, spiralarm_index=2) self.bar = dict(x=Quantity([x_0, x_1]), y=Quantity([y_0, y_1]))
# TODO: this is not picked up in the HTML docs ... don't know why. # http://sphinx-doc.org/latest/ext/example_numpy.html # For now I add it in the module-level docstring in an `Attributes` section. radial_distributions = OrderedDict() """Dictionary of available spatial distributions. Useful for automatic processing. """ radial_distributions['CB98'] = CaseBattacharya1998 radial_distributions['F06'] = FaucherKaspi2006 radial_distributions['L06'] = Lorimer2006 radial_distributions['P90'] = Paczynski1990 radial_distributions['YK04'] = YusifovKucuk2004 radial_distributions['YK04B'] = YusifovKucuk2004B