Projecting your population onto the sky

In this tutorial we’ll discuss how to project your simulated population onto the sky.

Learning Goals

By the end of this tutorial you should know how to:

• Make a simple sky projection plot

• Use healpy to make a sky projection plot

import cogsworth
import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
import astropy.units as u

# this all just makes plots look nice
%config InlineBackend.figure_format = 'retina'

plt.style.use("dark_background")
plt.rc('font', family='serif')
plt.rcParams['text.usetex'] = False
fs = 24

# update various fontsizes to match
params = {'figure.figsize': (12, 8),
          'legend.fontsize': fs,
          'axes.labelsize': fs,
          'xtick.labelsize': 0.9 * fs,
          'ytick.labelsize': 0.9 * fs,
          'axes.linewidth': 1.1,
          'xtick.major.size': 7,
          'xtick.minor.size': 4,
          'ytick.major.size': 7,
          'ytick.minor.size': 4}
plt.rcParams.update(params)
pd.options.display.max_columns = 999

Sky scatter plot

Warning

All of these functions assume that the population is in the Milky Way and we are using the galacocentric frame

Before doing some a little fancier there are built-in functions in cogsworth to visualise your population in terms RA and Dec on the sky. Let’s make a simple population to show this.

p = cogsworth.pop.Population(1000, final_kstar1=[10, 11, 12, 13, 14],
                             use_default_BSE_settings=True)
p.create_population()

Run for 1000 binaries
Sampled 1147 binaries
[8e-03s] Sample initial binaries
[0.7s] Evolve binaries (run COSMIC)

Integrating orbits: 100%|██████████| 1159/1159 [00:01<00:00, 630.48it/s]
cogsworth warning: 1 orbit(s) failed numerical integration, removing them. This can occur due to NaNs in stellar evolution or extreme orbits (e.g. passing directly through the galactic centre) that Gala cannot handle. Information for these systems was saved to `./failed_integration_binaries.h5`. This includes their initial_binaries, bpp, kick_info, and initial galaxy objects.

[2.8s] Integrate galactic orbits (run gala)
Overall: 3.5s

From this population we can find the present-day astropy SkyCoord object with

final_coords = p.get_final_mw_skycoord()
final_coords

<SkyCoord (Galactocentric: galcen_coord=<ICRS Coordinate: (ra, dec) in deg
    (266.4051, -28.936175)>, galcen_distance=8.122 kpc, galcen_v_sun=(12.9, 245.6, 7.78) km / s, z_sun=20.8 pc, roll=0.0 deg): (x, y, z) in kpc
    [( 8.69623465e-04, -2.48169157e+00,  8.79814106e-02),
     (-2.05077600e+00,  6.59705615e+00,  1.11441218e-01),
     (-3.56046654e-01, -3.35589184e-01,  3.68610176e-01), ...,
     (-1.51882266e+00,  1.44614010e+00,  9.65844024e-02),
     ( 2.48292926e+00, -9.26150352e-01,  4.03096205e-01),
     (-4.92179718e+02, -2.69582788e+03, -3.74169118e+03)]
 (v_x, v_y, v_z) in km / s
    [( 238.29314276,  117.55031534,   25.88750096),
     (-203.42235805, -142.79891164,  -27.25261596),
     ( -13.60875036,  -22.01348711,  110.63308024), ...,
     (  91.70709999,  -98.54905486,   -7.78366447),
     ( -21.39116307,  157.38803711,   34.36064019),
     ( -54.67050019, -299.4494503 , -415.64074305)]>

You can transform this coordinate in all sorts of ways with astropy. Particularly relevant in this case is transforming to icrs, which lets us find the RA and Dec of each source.

final_coords.icrs

<SkyCoord (ICRS): (ra, dec, distance) in (deg, deg, kpc)
    [(253.77562462, -42.54925741,    8.49375563),
     (288.94343069,  13.14471962,    8.96598032),
     (262.26397963, -29.57172879,    7.78095172), ...,
     (272.42087055, -17.85073957,    6.76007868),
     (261.1236937 , -31.94215788,   10.65212957),
     ( 50.49556232, -48.83575608, 4637.0473851 )]
 (pm_ra_cosdec, pm_dec, radial_velocity) in (mas / yr, mas / yr, km / s)
    [(-1.19405997, -8.42952486e-01,  253.10704887),
     (-0.50355342, -2.49707110e+00, -432.61310189),
     (-6.35411731, -4.52525263e+00,  -10.31795026), ...,
     (-4.9059492 , -9.87605612e+00,    3.18160871),
     (-1.46127734, -1.18290625e+00,  -25.51557648),
     ( 0.00754942, -4.63407964e-03,  665.59293932)]>

Now we don’t actually need to grab these coordinates, we can let cogsworth do that under the hood and just plot the sky positions

p.plot_sky_locations()

(<Figure size 1200x800 with 1 Axes>,
 <Axes: xlabel='Right Ascension [deg]', ylabel='Declination [deg]'>)

We could also style this differently to colour the points by their distance from the galactic centre. The function allows you to pass any kwargs you like to the .scatter call so we can restyle the point a bit as well.

dist = (final_coords.x**2 + final_coords.y**2 + final_coords.z**2)**(0.5)

# define a wider axis
fig, ax = plt.subplots(figsize=(15, 7))

# plot the points coloured by distance (and make them a little smaller)
fig, ax = p.plot_sky_locations(s=5, c=dist.value, vmin=0, vmax=20, fig=fig, ax=ax, show=False)

# add a colourbar for the scatter
fig.colorbar(ax.collections[0], extend='max', label="Galactocentric distance [kpc]")
plt.show()

We can even highlight the disrupted objects as well!

dist = (final_coords.x**2 + final_coords.y**2 + final_coords.z**2)**(0.5)

# define a wider axis
fig, ax = plt.subplots(figsize=(15, 7))

# plot the points coloured by distance (and make them a little smaller)
fig, ax = p.plot_sky_locations(s=5, c=dist.value, vmin=0, vmax=20, fig=fig, ax=ax, show=False)

disrupted_pop = p[p.disrupted]
disrupted_pop.plot_sky_locations(s=100, edgecolor="red", facecolor="none", fig=fig, ax=ax, show=False,
                                 label="Disrupted systems", show_galactic_plane=False)

# add a colourbar for the scatter
fig.colorbar(ax.collections[0], extend='max', label="Galactocentric distance [kpc]")
plt.show()

Healpix sky plot

Now let’s assume you’ve got a population that’s way larger and can’t be plotted as simple scatter points - in this case something like a healpix density plot could be useful. Let’s make one!

p.plot_map(ra="auto", dec="auto", nside=8)

You can also change the coordinates or colourmap of this with the various parameters

p.plot_map(ra="auto", dec="auto", nside=64, norm="linear", coord="G", cmap="cividis")

Wrap-up

And now you can see all of your glorious binaries on the sky! Looking forward to seeing you next time for another cogsworth tutorial, remember to like, comment and subscribe…or since this isn’t YouTube maybe just give me a star on GitHub…or just take a moment and smile to yourself that works too (:

Note

This tutorial was generated from a Jupyter notebook that can be found here.