This is a fixed-text formatted version of a Jupyter notebook.

Try online Binder You can also contribute with your own notebooks in this GitHub repository. Source files: source_population_model.ipynb |

Astrophysical source population modeling with Gammapy


The gammapy.astro.population package contains some simple Galactic source population models.

Here we provide some Python code to compute observable parameter distributions for Galactic gamma-ray source populations.

  • Observables: Flux, GLON, GLAT
  • Source classes: Pulsar (PSR), Supernova remnant (SNR), pulsar wind nebula (PWN)



In [1]:
%matplotlib inline
import matplotlib.pyplot as plt
In [2]:
import numpy as np
import astropy.units as u
from gammapy.utils.random import sample_powerlaw
from gammapy.astro import population

Simulate positions

In [3]:
# Spatial distribution using Lorimer (2006) model
n_sources = int(1e5)

table = population.make_base_catalog_galactic(
    max_age=1e6 * u.yr,

Simulate luminosities

Several source population models, e.g. the 1FHL paper or Strong (2007), use power-law luminosity functions.

Here we implement the “reference model” from the 1FHL catalog paper section 6.2.

In [4]:
# Source luminosity (ph s^-1)

luminosity = sample_powerlaw(
table['luminosity'] = luminosity

Compute observable parameters

In [5]:
table = population.add_observed_parameters(table)
<Table length=100000>
   name     dtype     unit                description
---------- ------- --------- --------------------------------------
       age float64        yr                      Age of the source
     n_ISM float64   1 / cm3            Interstellar medium density
 spiralarm   str18                                 Which spiralarm?
   x_birth float64       kpc   Galactocentric x coordinate at birth
   y_birth float64       kpc   Galactocentric y coordinate at birth
   z_birth float64       kpc   Galactocentric z coordinate at birth
         x float64       kpc            Galactocentric x coordinate
         y float64       kpc            Galactocentric y coordinate
         z float64       kpc            Galactocentric z coordinate
        vx float64    km / s Galactocentric velocity in x direction
        vy float64    km / s Galactocentric velocity in y direction
        vz float64    km / s Galactocentric velocity in z direction
     v_abs float64    km / s     Galactocentric velocity (absolute)
luminosity float64
  distance float64        pc     Distance observer to source center
      GLON float64       deg                     Galactic longitude
      GLAT float64       deg                      Galactic latitude
     VGLON float64 deg / Myr         Velocity in Galactic longitude
     VGLAT float64 deg / Myr          Velocity in Galactic latitude
        RA float64       deg                        Right ascension
       DEC float64       deg                            Declination
      flux float64  1 / kpc2                            Source flux

Check output

The simulation is done, you could save the simulated catalog to a file.

Here we just plot a few distributions to check if the results look OK.

In [6]:
plt.scatter(table['x'][:1000], table['y'][:1000]);
In [7]:
plt.hist(table['GLON'], bins=100);
In [8]:
plt.hist(table['GLAT'], bins=100, log=True);
In [9]:
plt.scatter(table['GLON'][:1000], table['GLAT'][:1000]);
In [10]:
plt.hist(table['distance'], bins=100, log=True);
In [11]:
plt.hist(np.log10(table['luminosity']), bins=100, log=True);
In [12]:
plt.hist(np.log10(table['flux']), bins=100, log=True);
In [13]:
# TODO: plot GLON, GLAT, FLUX distribution



In [14]:
# Start exercises here

What next?

TODO: summarise what was done here briefly.

TODO: add some pointers to other documentation.