Source code for gammapy.astro.population.simulate

# Licensed under a 3-clause BSD style license - see LICENSE.rst
"""Simulate source catalogs."""
import numpy as np
from astropy.table import Table, Column
from astropy.units import Quantity
from astropy.coordinates import SkyCoord, spherical_to_cartesian
from ...utils import coordinates as astrometry
from ...utils.coordinates import D_SUN_TO_GALACTIC_CENTER
from ...utils.distributions import draw, pdf
from ...utils.random import sample_sphere, sample_sphere_distance, get_random_state
from ..source import SNR, SNRTrueloveMcKee, PWN, Pulsar
from ..population.spatial import (
    Exponential,
    FaucherSpiral,
    RMIN,
    RMAX,
    ZMIN,
    ZMAX,
    radial_distributions,
)
from ..population.velocity import VMIN, VMAX, velocity_distributions

__all__ = [
    "make_catalog_random_positions_cube",
    "make_catalog_random_positions_sphere",
    "make_base_catalog_galactic",
    "add_snr_parameters",
    "add_pulsar_parameters",
    "add_pwn_parameters",
    "add_observed_source_parameters",
    "add_observed_parameters",
]


[docs]def make_catalog_random_positions_cube( size=100, dimension=3, dmax=10, random_state="random-seed" ): """Make a catalog of sources randomly distributed on a line, square or cube. TODO: is this useful enough for general use or should we hide it as an internal method to generate test datasets? Parameters ---------- size : int, optional Number of sources dimension : int, optional Number of dimensions dmax : int, optional Maximum distance in pc. random_state : {int, 'random-seed', 'global-rng', `~numpy.random.RandomState`} Defines random number generator initialisation. Passed to `~gammapy.utils.random.get_random_state`. Returns ------- catalog : `~astropy.table.Table` Source catalog with columns: """ random_state = get_random_state(random_state) # Generate positions 1D, 2D, or 3D if dimension == 3: x = random_state.uniform(-dmax, dmax, size) y = random_state.uniform(-dmax, dmax, size) z = random_state.uniform(-dmax, dmax, size) elif dimension == 2: x = random_state.uniform(-dmax, dmax, size) y = random_state.uniform(-dmax, dmax, size) z = np.zeros_like(x) else: x = random_state.uniform(-dmax, dmax, size) y = np.zeros_like(x) z = np.zeros_like(x) table = Table() table["x"] = Column(x, unit="pc", description="Galactic cartesian coordinate") table["y"] = Column(y, unit="pc", description="Galactic cartesian coordinate") table["z"] = Column(z, unit="pc", description="Galactic cartesian coordinate") return table
[docs]def make_catalog_random_positions_sphere( size, center="Earth", distance=Quantity([0, 1], "Mpc"), random_state="random-seed" ): """Sample random source locations in a sphere. This can be used to generate an isotropic source population to represent extra-galactic sources. Parameters ---------- size : int Number of sources center : {'Earth', 'Milky Way'} Sphere center distance : `~astropy.units.Quantity` tuple Distance min / max range. random_state : {int, 'random-seed', 'global-rng', `~numpy.random.RandomState`} Defines random number generator initialisation. Passed to `~gammapy.utils.random.get_random_state`. Returns ------- catalog : `~astropy.table.Table` Source catalog with columns: - RAJ2000, DEJ2000 (deg) - GLON, GLAT (deg) - Distance (Mpc) """ random_state = get_random_state(random_state) lon, lat = sample_sphere(size, random_state=random_state) radius = sample_sphere_distance( distance[0], distance[1], size, random_state=random_state ) pos = SkyCoord(lon, lat, distance=radius, frame="galactic") if center == "Milky Way": pass elif center == "Earth": # TODO: add shift Galactic center -> Earth raise NotImplementedError else: msg = "Invalid center: {}\n".format(center) msg += "Choose one of: Earth, Milky Way" raise ValueError(msg) table = Table() table.meta["center"] = center icrs = pos.transform_to("icrs") table["RAJ2000"] = icrs.ra.to("deg") table["DEJ2000"] = icrs.dec.to("deg") galactic = icrs.transform_to("galactic") table["GLON"] = galactic.l.to("deg") table["GLAT"] = galactic.b.to("deg") table["Distance"] = icrs.distance.to("Mpc") return table
[docs]def make_base_catalog_galactic( n_sources, rad_dis="YK04", vel_dis="H05", max_age=Quantity(1e6, "yr"), spiralarms=True, n_ISM=Quantity(1, "cm-3"), random_state="random-seed", ): """Make a catalog of Galactic sources, with basic source parameters. Choose a radial distribution, a velocity distribution, the number of pulsars n_pulsars, the maximal age max_age[years] and the fraction of the individual morphtypes. There's an option spiralarms. If set on True a spiralarm modelling after Faucher&Kaspi is included. max_age and n_sources effectively correspond to s SN rate: SN_rate = n_sources / max_age Parameters ---------- n_sources : int Number of sources to simulate. rad_dis : callable Radial surface density distribution of sources. vel_dis : callable Proper motion velocity distribution of sources. max_age : `~astropy.units.Quantity` Maximal age of the source spiralarms : bool Include a spiralarm model in the catalog. n_ISM : `~astropy.units.Quantity` Density of the interstellar medium. random_state : {int, 'random-seed', 'global-rng', `~numpy.random.RandomState`} Defines random number generator initialisation. Passed to `~gammapy.utils.random.get_random_state`. Returns ------- table : `~astropy.table.Table` Catalog of simulated source positions and proper velocities. """ random_state = get_random_state(random_state) if isinstance(rad_dis, str): rad_dis = radial_distributions[rad_dis] if isinstance(vel_dis, str): vel_dis = velocity_distributions[vel_dis] # Draw random values for the age age = random_state.uniform(0, max_age.to_value("yr"), n_sources) age = Quantity(age, "yr") # Draw r and z values from the given distribution r = draw( RMIN.to_value("kpc"), RMAX.to_value("kpc"), n_sources, pdf(rad_dis()), random_state=random_state, ) r = Quantity(r, "kpc") z = draw( ZMIN.to_value("kpc"), ZMAX.to_value("kpc"), n_sources, Exponential(), random_state=random_state, ) z = Quantity(z, "kpc") # Apply spiralarm modelling or not if spiralarms: r, theta, spiralarm = FaucherSpiral()(r, random_state=random_state) else: theta = Quantity(random_state.uniform(0, 2 * np.pi, n_sources), "rad") spiralarm = None # Compute cartesian coordinates x, y = astrometry.cartesian(r, theta) # Draw values from velocity distribution v = draw( VMIN.to_value("km/s"), VMAX.to_value("km/s"), n_sources, vel_dis(), random_state=random_state, ) v = Quantity(v, "km/s") # Draw random direction of initial velocity theta = Quantity(random_state.uniform(0, np.pi, x.size), "rad") phi = Quantity(random_state.uniform(0, 2 * np.pi, x.size), "rad") # Compute new position dx, dy, dz, vx, vy, vz = astrometry.motion_since_birth(v, age, theta, phi) # Add displacement to birth position x_moved = x + dx y_moved = y + dy z_moved = z + dz # Set environment interstellar density n_ISM = n_ISM * np.ones(n_sources) table = Table() table["age"] = Column(age, unit="yr", description="Age of the source") table["n_ISM"] = Column( n_ISM, unit="cm-3", description="Interstellar medium density" ) if spiralarms: table["spiralarm"] = Column(spiralarm, description="Which spiralarm?") table["x_birth"] = Column( x, unit="kpc", description="Galactocentric x coordinate at birth" ) table["y_birth"] = Column( y, unit="kpc", description="Galactocentric y coordinate at birth" ) table["z_birth"] = Column( z, unit="kpc", description="Galactocentric z coordinate at birth" ) table["x"] = Column( x_moved.to("kpc"), unit="kpc", description="Galactocentric x coordinate" ) table["y"] = Column( y_moved.to("kpc"), unit="kpc", description="Galactocentric y coordinate" ) table["z"] = Column( z_moved.to("kpc"), unit="kpc", description="Galactocentric z coordinate" ) table["vx"] = Column( vx.to("km/s"), unit="km/s", description="Galactocentric velocity in x direction" ) table["vy"] = Column( vy.to("km/s"), unit="km/s", description="Galactocentric velocity in y direction" ) table["vz"] = Column( vz.to("km/s"), unit="km/s", description="Galactocentric velocity in z direction" ) table["v_abs"] = Column( v, unit="km/s", description="Galactocentric velocity (absolute)" ) return table
[docs]def add_snr_parameters(table): """Add SNR parameters to the table.""" # Read relevant columns age = table["age"].quantity n_ISM = table["n_ISM"].quantity # Compute properties snr = SNR(n_ISM=n_ISM) E_SN = snr.e_sn * np.ones(len(table)) r_out = snr.radius(age) r_in = snr.radius_inner(age) L_SNR = snr.luminosity_tev(age) # Add columns to table table["E_SN"] = Column(E_SN, unit="erg", description="SNR kinetic energy") table["r_out"] = Column(r_out, unit="pc", description="SNR outer radius") table["r_in"] = Column(r_in, unit="pc", description="SNR inner radius") table["L_SNR"] = Column(L_SNR, unit="s-1", description="SNR luminosity") return table
[docs]def add_pulsar_parameters( table, B_mean=12.05, B_stdv=0.55, P_mean=0.3, P_stdv=0.15, random_state="random-seed", ): """Add pulsar parameters to the table. For the initial normal distribution of period and logB can exist the following Parameters: B_mean=12.05[log Gauss], B_stdv=0.55, P_mean=0.3[s], P_stdv=0.15 Parameters ---------- 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) # Read relevant columns age = table["age"].quantity # Draw the initial values for the period and magnetic field def p_dist(x): return np.exp(-0.5 * ((x - P_mean) / P_stdv) ** 2) p0_birth = draw(0, 2, len(table), p_dist, random_state=random_state) p0_birth = Quantity(p0_birth, "s") logB = random_state.normal(B_mean, B_stdv, len(table)) # Compute pulsar parameters psr = Pulsar(p0_birth, logB) p0 = psr.period(age) p1 = psr.period_dot(age) p1_birth = psr.P_dot_0 tau = psr.tau(age) tau_0 = psr.tau_0 l_psr = psr.luminosity_spindown(age) l0_psr = psr.L_0 # Add columns to table table["P0"] = Column(p0, unit="s", description="Pulsar period") table["P1"] = Column(p1, unit="", description="Pulsar period derivative") table["P0_birth"] = Column(p0_birth, unit="s", description="Pulsar birth period") table["P1_birth"] = Column( p1_birth, unit="", description="Pulsar birth period derivative" ) table["CharAge"] = Column(tau, unit="yr", description="Pulsar characteristic age") table["Tau0"] = Column(tau_0, unit="yr") table["L_PSR"] = Column(l_psr, unit="erg s-1") table["L0_PSR"] = Column(l0_psr, unit="erg s-1") table["logB"] = Column(logB, unit="Gauss") return table
[docs]def add_pwn_parameters(table): """Add PWN parameters to the table.""" # Some of the computations (specifically `pwn.radius`) aren't vectorised # across all parameters; so here we loop over source parameters explicitly results = [] for idx in range(len(table)): age = table["age"].quantity[idx] E_SN = table["E_SN"].quantity[idx] n_ISM = table["n_ISM"].quantity[idx] P0_birth = table["P0_birth"].quantity[idx] logB = table["logB"][idx] # Compute properties pulsar = Pulsar(P0_birth, logB) snr = SNRTrueloveMcKee(e_sn=E_SN, n_ISM=n_ISM) pwn = PWN(pulsar, snr) r_out_pwn = pwn.radius(age).to_value("pc") L_PWN = pwn.luminosity_tev(age).to_value("erg") results.append(dict(r_out_pwn=r_out_pwn, L_PWN=L_PWN)) # Add columns to table table["r_out_PWN"] = Column( [_["r_out_pwn"] for _ in results], unit="pc", description="PWN outer radius" ) table["L_PWN"] = Column( [_["L_PWN"] for _ in results], unit="erg", description="PWN luminosity above 1 TeV", ) return table
[docs]def add_observed_source_parameters(table): """Add observed source parameters to the table.""" # Read relevant columns distance = table["distance"] r_in = table["r_in"] r_out = table["r_out"] r_out_PWN = table["r_out_PWN"] L_SNR = table["L_SNR"] L_PSR = table["L_PSR"] L_PWN = table["L_PWN"] # Compute properties ext_in_SNR = astrometry.radius_to_angle(r_in, distance) ext_out_SNR = astrometry.radius_to_angle(r_out, distance) ext_out_PWN = astrometry.radius_to_angle(r_out_PWN, distance) # Ellipse parameters not used for now theta = np.pi / 2 * np.ones(len(table)) # Position angle? epsilon = np.zeros(len(table)) # Ellipticity? S_SNR = astrometry.luminosity_to_flux(L_SNR, distance) # Ld2_PSR = astrometry.luminosity_to_flux(L_PSR, distance) Ld2_PSR = L_PSR / distance ** 2 S_PWN = astrometry.luminosity_to_flux(L_PWN, distance) # Add columns table["ext_in_SNR"] = Column(ext_in_SNR, unit="deg") table["ext_out_SNR"] = Column(ext_out_SNR, unit="deg") table["ext_out_PWN"] = Column(ext_out_PWN, unit="deg") table["theta"] = Column(theta, unit="rad") table["epsilon"] = Column(epsilon, unit="") table["S_SNR"] = Column(S_SNR, unit="cm-2 s-1") table["Ld2_PSR"] = Column(Ld2_PSR, unit="erg s-1 kpc-2") table["S_PWN"] = Column(S_PWN, unit="cm-2 s-1") return table
[docs]def add_observed_parameters(table, obs_pos=None): """Add observable parameters (such as sky position or distance). Input table columns: x, y, z, extension, luminosity Output table columns: distance, glon, glat, flux, angular_extension Position of observer in cartesian coordinates. Center of galaxy as origin, x-axis goes trough sun. Parameters ---------- table : `~astropy.table.Table` Input table obs_pos : tuple or None Observation position (X, Y, Z) in Galactocentric coordinates (default: Earth) Returns ------- table : `~astropy.table.Table` Modified input table with columns added """ obs_pos = obs_pos or [D_SUN_TO_GALACTIC_CENTER, 0, 0] # Get data x, y, z = table["x"].quantity, table["y"].quantity, table["z"].quantity vx, vy, vz = table["vx"].quantity, table["vy"].quantity, table["vz"].quantity distance, glon, glat = astrometry.galactic(x, y, z, obs_pos=obs_pos) # Compute projected velocity v_glon, v_glat = astrometry.velocity_glon_glat(x, y, z, vx, vy, vz) coordinate = SkyCoord(glon, glat, unit="deg", frame="galactic").transform_to("icrs") ra, dec = coordinate.ra.deg, coordinate.dec.deg # Add columns to table table["distance"] = Column( distance, unit="pc", description="Distance observer to source center" ) table["GLON"] = Column(glon, unit="deg", description="Galactic longitude") table["GLAT"] = Column(glat, unit="deg", description="Galactic latitude") table["VGLON"] = Column( v_glon.to("deg/Myr"), unit="deg/Myr", description="Velocity in Galactic longitude", ) table["VGLAT"] = Column( v_glat.to("deg/Myr"), unit="deg/Myr", description="Velocity in Galactic latitude", ) table["RA"] = Column(ra, unit="deg", description="Right ascension") table["DEC"] = Column(dec, unit="deg", description="Declination") try: luminosity = table["luminosity"] flux = astrometry.luminosity_to_flux(luminosity, distance) table["flux"] = Column(flux.value, unit=flux.unit, description="Source flux") except KeyError: pass try: extension = table["extension"] angular_extension = np.degrees(np.arctan(extension / distance)) table["angular_extension"] = Column( angular_extension, unit="deg", description="Source angular radius (i.e. half-diameter)", ) except KeyError: pass return table