Quickstart#

Welcome! You’ve found the documentation for cogsworth, good for you! This is an epic crossover episode of a package in which we bring the worlds of rapid population synthesis and galactic dynamics together and see what happens.

You should know that package relies heavily on cosmic and gala and I recommend you check out their respective documentation. As a general rule, anything that you can do with either of those packages, you can also do with cogsworth!

Okay, let’s get into this!

Basic imports#

For most of the examples that I give in this documentation you’re probably going to need the following basic imports. You can always read more about numpy and astropy.units at those respective links.

[1]:

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

Here are some plotting settings as well to make sure everything looks good.

[2]:

# ensure jupyter actually uses your fancy retina display
%config InlineBackend.figure_format = 'retina'

# make pandas show *every* column
pd.set_option("display.max_columns", None)

# various adjustments to matplotlib settings
plt.rc('font', family='serif')
plt.rcParams['text.usetex'] = False
fs = 24
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)

Create your first population#

The central class in this package is Population. The class gives you access to the functionality of the entire package with one convenient interface. Let’s start by creating a small population of binaries with the default settings.

[3]:

import cogsworth
p = cogsworth.pop.Population(2000, processes=6,
                             use_default_BSE_settings=True)
p

[3]:

<Population - 2000 systems - galactic_potential=MilkyWayPotential, SFH=Wagg2022>

Though this class is initialised, we still need to actually create the population to find its present day state. We do this by running:

[4]:

p.create_population()

Run for 2000 binaries
Sampled 2642 binaries
[1e-02s] Sample initial binaries
[0.7s] Evolve binaries (run COSMIC)

Integrating orbits: 100%|██████████| 2647/2647 [00:03<00:00, 681.63it/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.

[6.1s] Integrate galactic orbits (run gala)
Overall: 6.7s

This function has just

  • Sampled the initial binaries from cosmic

  • Drawn their initial kinematics and birth times from the Galaxy model

  • Evolved the binaries from birth until present day using cosmic

  • Integrated the orbits of the binaries through the galaxy using gala

Inspect population#

Full evolution#

p now contains information about the stellar evolution of each constituent star of each binary, as well as the path of the binary through the galaxy. Let’s investigate some of those things.

Stellar Evolution#

First we can take a look at the cosmic bpp evolution table and see how the stars evolved over time.

[5]:

p.bpp

[5]:
tphys mass_1 mass_2 kstar_1 kstar_2 sep porb ecc RRLO_1 RRLO_2 evol_type aj_1 aj_2 tms_1 tms_2 massc_he_layer_1 massc_he_layer_2 massc_co_layer_1 massc_co_layer_2 rad_1 rad_2 mass0_1 mass0_2 lum_1 lum_2 teff_1 teff_2 radc_1 radc_2 menv_1 menv_2 renv_1 renv_2 omega_spin_1 omega_spin_2 B_1 B_2 bacc_1 bacc_2 tacc_1 tacc_2 epoch_1 epoch_2 bhspin_1 bhspin_2 bin_num
0 0.000000 7.568863 3.368652 1 1 112.880258 42.028107 0.187905 0.082614 0.075268 1 0.000000 0.000000 41.283771 281.760043 0.000000 0.0 0.0 0.0 3.410712 2.149079 7.568863 3.368652 2225.798869 121.550322 21563.465800 13132.048623 0.000000 0.0 1.000000e-10 1.000000e-10 1.000000e-10 1.000000e-10 3181.733547 6272.932458 0.0 0.0 0.0 0.0 0.0 0.0 0.000000 0.000000 0.0 0.0 0
0 41.561348 7.500004 3.368874 2 1 113.571698 42.548572 0.187891 0.192258 0.080393 2 42.074918 41.555171 42.074914 281.711325 1.459549 0.0 0.0 0.0 7.971621 2.314931 7.500004 3.368874 6666.427113 129.995901 18555.367844 12867.168292 0.269674 0.0 1.000000e-10 1.000000e-10 1.000000e-10 1.000000e-10 672.346521 5407.380217 0.0 0.0 0.0 0.0 0.0 0.0 -0.513570 0.006177 0.0 0.0 0
0 41.641493 7.499482 3.368889 2 1 92.236327 31.141744 0.000000 1.001606 0.080400 3 42.161168 41.634801 42.081011 281.707840 1.477410 0.0 0.0 0.0 41.530736 2.315280 7.499482 3.368889 4577.204119 130.016544 7400.053621 12866.710213 0.271994 0.0 1.000000e-10 1.000000e-10 1.000000e-10 1.000000e-10 73.691140 73.691140 0.0 0.0 0.0 0.0 0.0 0.0 -0.519675 0.006692 0.0 0.0 0
0 41.676451 7.276935 3.591246 2 1 87.313307 28.682327 0.000000 2.204458 0.085436 5 42.604472 33.085896 42.488101 239.094055 1.479350 0.0 0.0 0.0 85.033842 2.386717 7.464898 3.591246 3827.935591 164.224631 4945.553064 13434.723457 0.272246 0.0 2.676327e-04 1.000000e-10 5.400612e+00 1.000000e-10 55.137767 8136.018270 0.0 0.0 0.0 0.0 0.0 0.0 -0.928022 8.590423 0.0 0.0 0
0 41.676451 7.276935 3.591246 2 1 87.313307 28.682327 0.000000 2.204458 0.085436 7 42.604472 33.085896 42.488101 239.094055 1.479350 0.0 0.0 0.0 85.033842 2.386717 7.464898 3.591246 3827.935591 164.224631 4945.553064 13434.723457 0.272246 0.0 2.676327e-04 1.000000e-10 5.400612e+00 1.000000e-10 55.137767 8136.018270 0.0 0.0 0.0 0.0 0.0 0.0 -0.928022 8.590423 0.0 0.0 0
... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ...
2639 7377.627266 0.929110 0.376353 1 0 6.687098 1.754068 0.000000 0.299791 0.177005 10 7377.627266 7377.627266 15966.648407 257999.381374 0.000000 0.0 0.0 0.0 0.919678 0.359876 0.929110 0.376353 0.649069 0.019640 5426.556898 3618.099715 0.000000 0.0 4.056902e-02 3.077718e-01 2.590509e-01 2.763195e-01 1307.850867 451.070064 0.0 0.0 0.0 0.0 0.0 0.0 0.000000 0.000000 0.0 0.0 2639
2640 0.000000 0.353198 0.293732 0 0 4.649638 1.444682 0.691949 0.599647 0.572766 1 0.000000 0.000000 288525.963202 422473.944836 0.000000 0.0 0.0 0.0 0.339244 0.297862 0.353198 0.293732 0.017701 0.012464 3630.909140 3549.545954 0.000000 0.0 3.531984e-01 2.937316e-01 3.392440e-01 2.978616e-01 0.517788 0.139226 0.0 0.0 0.0 0.0 0.0 0.0 0.000000 0.000000 0.0 0.0 2640
2640 9296.792609 0.353198 0.293732 0 0 3.097413 0.785493 0.535276 0.600404 0.572358 10 9296.792609 9296.792609 288525.963202 422473.944836 0.000000 0.0 0.0 0.0 0.341359 0.299128 0.353198 0.293732 0.017899 0.012556 3629.701735 3548.550016 0.000000 0.0 3.531984e-01 2.937316e-01 3.413592e-01 2.991282e-01 9110.216250 8764.228392 0.0 0.0 0.0 0.0 0.0 0.0 0.000000 0.000000 0.0 0.0 2640
2641 0.000000 0.260330 0.237264 0 0 700.706127 3047.463102 0.049573 0.001059 0.001038 1 0.000000 0.000000 526134.830937 634328.740372 0.000000 0.0 0.0 0.0 0.272888 0.256502 0.260330 0.237264 0.010255 0.008525 3531.953066 3478.499708 0.000000 0.0 2.603296e-01 2.372644e-01 2.728881e-01 2.565017e-01 0.059054 0.030385 0.0 0.0 0.0 0.0 0.0 0.0 0.000000 0.000000 0.0 0.0 2641
2641 8494.843225 0.260330 0.237264 0 0 700.706127 3047.463102 0.049573 0.001062 0.001041 10 8494.843225 8494.843225 526134.830937 634328.740372 0.000000 0.0 0.0 0.0 0.273794 0.257207 0.260330 0.237264 0.010310 0.008562 3530.806463 3477.549784 0.000000 0.0 2.603296e-01 2.372644e-01 2.737936e-01 2.572074e-01 0.058664 0.030219 0.0 0.0 0.0 0.0 0.0 0.0 0.000000 0.000000 0.0 0.0 2641

8729 rows × 46 columns

There are of course many things that you can do with this table! But let’s say, maybe you want get any binaries that experience mass transfer

[6]:

# get any binaries that have an evol_type indicating Roche Lobe overflow started
mt_bin_nums = p.bpp[p.bpp["evol_type"] == 3.0]["bin_num"].unique()

Let’s take a look at the evolution of the first binary in this list

[7]:

p.bpp[p.bpp["bin_num"] == mt_bin_nums[0]]

[7]:
tphys mass_1 mass_2 kstar_1 kstar_2 sep porb ecc RRLO_1 RRLO_2 evol_type aj_1 aj_2 tms_1 tms_2 massc_he_layer_1 massc_he_layer_2 massc_co_layer_1 massc_co_layer_2 rad_1 rad_2 mass0_1 mass0_2 lum_1 lum_2 teff_1 teff_2 radc_1 radc_2 menv_1 menv_2 renv_1 renv_2 omega_spin_1 omega_spin_2 B_1 B_2 bacc_1 bacc_2 tacc_1 tacc_2 epoch_1 epoch_2 bhspin_1 bhspin_2 bin_num
0 0.000000 7.568863 3.368652 1 1 112.880258 42.028107 0.187905 0.082614 0.075268 1 0.000000 0.000000 4.128377e+01 281.760043 0.000000 0.0 0.000000 0.0 3.410712 2.149079 7.568863 3.368652 2.225799e+03 121.550322 2.156347e+04 13132.048623 0.000000 0.0 1.000000e-10 1.000000e-10 1.000000e-10 1.000000e-10 3.181734e+03 6272.932458 0.000000e+00 0.0 0.0 0.0 0.0 0.0 0.000000 0.000000 0.0 0.0 0
0 41.561348 7.500004 3.368874 2 1 113.571698 42.548572 0.187891 0.192258 0.080393 2 42.074918 41.555171 4.207491e+01 281.711325 1.459549 0.0 0.000000 0.0 7.971621 2.314931 7.500004 3.368874 6.666427e+03 129.995901 1.855537e+04 12867.168292 0.269674 0.0 1.000000e-10 1.000000e-10 1.000000e-10 1.000000e-10 6.723465e+02 5407.380217 0.000000e+00 0.0 0.0 0.0 0.0 0.0 -0.513570 0.006177 0.0 0.0 0
0 41.641493 7.499482 3.368889 2 1 92.236327 31.141744 0.000000 1.001606 0.080400 3 42.161168 41.634801 4.208101e+01 281.707840 1.477410 0.0 0.000000 0.0 41.530736 2.315280 7.499482 3.368889 4.577204e+03 130.016544 7.400054e+03 12866.710213 0.271994 0.0 1.000000e-10 1.000000e-10 1.000000e-10 1.000000e-10 7.369114e+01 73.691140 0.000000e+00 0.0 0.0 0.0 0.0 0.0 -0.519675 0.006692 0.0 0.0 0
0 41.676451 7.276935 3.591246 2 1 87.313307 28.682327 0.000000 2.204458 0.085436 5 42.604472 33.085896 4.248810e+01 239.094055 1.479350 0.0 0.000000 0.0 85.033842 2.386717 7.464898 3.591246 3.827936e+03 164.224631 4.945553e+03 13434.723457 0.272246 0.0 2.676327e-04 1.000000e-10 5.400612e+00 1.000000e-10 5.513777e+01 8136.018270 0.000000e+00 0.0 0.0 0.0 0.0 0.0 -0.928022 8.590423 0.0 0.0 0
0 41.676451 7.276935 3.591246 2 1 87.313307 28.682327 0.000000 2.204458 0.085436 7 42.604472 33.085896 4.248810e+01 239.094055 1.479350 0.0 0.000000 0.0 85.033842 2.386717 7.464898 3.591246 3.827936e+03 164.224631 4.945553e+03 13434.723457 0.272246 0.0 2.676327e-04 1.000000e-10 5.400612e+00 1.000000e-10 5.513777e+01 8136.018270 0.000000e+00 0.0 0.0 0.0 0.0 0.0 -0.928022 8.590423 0.0 0.0 0
0 41.676451 10.064979 0.000000 2 15 0.000000 0.000000 0.000000 0.990000 0.085436 8 23.749589 33.085896 4.248810e+01 239.094055 2.210050 0.0 0.000000 0.0 85.033842 2.386717 10.064979 3.591246 3.827936e+03 164.224631 4.945553e+03 13434.723457 0.272246 0.0 2.676327e-04 1.000000e-10 5.400612e+00 1.000000e-10 5.513777e+01 8136.018270 0.000000e+00 0.0 0.0 0.0 0.0 0.0 -0.928022 8.590423 0.0 0.0 0
0 41.683526 10.064378 0.000000 3 15 0.000000 0.000000 -1.000000 0.000100 -1.000000 2 23.759245 33.085896 2.369878e+01 239.094055 2.214953 0.0 0.000000 0.0 202.269279 2.386717 10.064378 3.591246 9.513949e+03 164.224631 4.026195e+03 13434.723457 0.358949 0.0 3.926174e+00 1.000000e-10 1.312680e+02 1.000000e-10 1.370865e+01 8136.018270 0.000000e+00 0.0 0.0 0.0 0.0 0.0 17.924281 8.590555 0.0 0.0 0
0 41.692590 10.062287 0.000000 4 15 0.000000 0.000000 -1.000000 0.000100 -1.000000 2 23.768309 33.085896 2.369878e+01 239.094055 2.215140 0.0 0.000000 0.0 408.816017 2.386717 10.064378 3.591246 2.554593e+04 164.224631 3.625234e+03 13434.723457 0.358969 0.0 7.847148e+00 1.000000e-10 4.084570e+02 1.000000e-10 3.350845e+00 8136.018270 0.000000e+00 0.0 0.0 0.0 0.0 0.0 17.924281 8.590555 0.0 0.0 0
0 44.400023 9.513885 0.000000 5 15 0.000000 0.000000 -1.000000 0.000100 -1.000000 2 26.475742 33.085896 2.369878e+01 239.094055 1.358975 0.0 1.561453 0.0 430.515092 2.386717 10.064378 3.591246 2.506781e+04 164.224631 3.516046e+03 13434.723457 0.430888 0.0 6.593457e+00 1.000000e-10 4.300842e+02 1.000000e-10 2.171259e+00 8136.018270 0.000000e+00 0.0 0.0 0.0 0.0 0.0 17.924281 8.590555 0.0 0.0 0
0 44.507277 9.364967 0.000000 5 15 0.000000 0.000000 -1.000000 0.000100 -1.000000 15 26.582992 33.085896 2.369878e+01 239.094055 1.012913 0.0 1.907491 0.0 864.048546 2.386717 10.064378 3.591246 6.819991e+04 164.224631 3.187473e+03 13434.723457 10.997961 0.0 4.985976e+00 1.000000e-10 8.000444e+02 1.000000e-10 4.331925e+15 8136.018270 0.000000e+00 0.0 0.0 0.0 0.0 0.0 17.924281 8.590555 0.0 0.0 0
0 44.507277 1.277584 0.000000 13 15 0.000000 0.000000 -1.000000 0.000100 -1.000000 2 0.000000 33.085896 2.369878e+01 239.094055 0.000000 0.0 0.000000 0.0 0.000014 2.386717 10.064378 3.591246 2.356735e+00 164.224631 1.919921e+06 13434.723457 0.000014 0.0 1.000000e-10 1.000000e-10 1.000000e-10 1.000000e-10 4.331925e+15 8136.018270 0.000000e+00 0.0 0.0 0.0 0.0 0.0 44.507277 8.590555 0.0 0.0 0
0 7599.107443 1.277584 0.000000 13 15 0.000000 0.000000 -1.000000 0.000100 -1.000000 10 7554.600165 33.085896 1.000000e+10 239.094055 0.000000 0.0 0.000000 0.0 0.000014 2.386717 10.064378 3.591246 4.129409e-10 164.224631 6.985179e+03 13434.723457 0.000014 0.0 1.000000e-10 1.000000e-10 1.000000e-10 1.000000e-10 1.616444e+07 8136.018270 1.262821e+11 0.0 0.0 0.0 0.0 0.0 44.507277 8.590555 0.0 0.0 0

We can also put this to use. Let’s see how the initial orbital periods are different for binaries that experience mass transfer.

[8]:

experienced_mt = p.initC["bin_num"].isin(mt_bin_nums)
plt.hist(np.log10(p.initC["porb"][experienced_mt]), bins="fd", density=True, label="Mass transfer")
plt.hist(np.log10(p.initC["porb"][~experienced_mt]), bins="fd", density=True, label="No mass transfer", alpha=0.8)
plt.legend()
plt.xlabel("Log Initial Orbital Period [days]");
../_images/pages_getting_started_16_0.png

As we can see in the plot above, systems that experience mass transfer are biased towards shorter orbital periods. This is expected as these systems are closer and therefore it is easier to initiate mass transfer.

Orbits#

We also have access to the full orbit of each binary in the population in p.orbits. I’m going to cheat and pick a fairly simple one to show from this random population (a binary that didn’t disrupt and is pretty close to us at present day).

[9]:

fig, ax = plt.subplots()
final_rho = np.sum(p.final_pos[:, :2]**2, axis=1)**(0.5)
nice_orbits = p.orbits[(final_rho > 4 * u.kpc) & (final_rho < 5 * u.kpc)]
nice_orbit = np.random.choice(nice_orbits)
nice_orbit.cylindrical.plot(["rho", "z"], axes=ax)
plt.show()
../_images/pages_getting_started_18_0.png

If this part interests you then I recommend that you go and check out the excellent gala docs about orbits and what you can do with them.

Present Day#

In many cases, you may not really care what happened over the history of the galaxy (way to let the past go!). In these cases you can instead only focus on the present day state of each binary.

For example, we could extract any binary that ended up disrupted

[10]:

p.final_bpp[p.disrupted]

[10]:
tphys mass_1 mass_2 kstar_1 kstar_2 sep porb ecc RRLO_1 RRLO_2 evol_type aj_1 aj_2 tms_1 tms_2 massc_he_layer_1 massc_he_layer_2 massc_co_layer_1 massc_co_layer_2 rad_1 rad_2 mass0_1 mass0_2 lum_1 lum_2 teff_1 teff_2 radc_1 radc_2 menv_1 menv_2 renv_1 renv_2 omega_spin_1 omega_spin_2 B_1 B_2 bacc_1 bacc_2 tacc_1 tacc_2 epoch_1 epoch_2 bhspin_1 bhspin_2 bin_num metallicity
318 4129.976422 1.277584 0.707752 13 11 -1.0 -1.0 -1.0 0.0001 0.0001 10 4105.077530 3756.497060 1.000000e+10 1.000000e+10 0.0 0.0 0.0 0.00 0.000014 0.011401 9.871907 0.707752 1.398518e-09 1.124590e-04 9475.946905 5591.721802 0.000014 0.011401 1.000000e-10 1.000000e-10 1.000000e-10 1.000000e-10 3.106852e+07 5.176179e-06 2.079877e+11 0.000000e+00 0.0 0.0 0.0 0.0 24.898892 373.479362 0.0 0.0 318 0.025573
1555 8912.425861 1.277584 1.277584 13 13 -1.0 -1.0 -1.0 0.0001 0.0001 10 8895.972934 8870.578887 1.000000e+10 1.000000e+10 0.0 0.0 0.0 0.00 0.000014 0.000014 3.126963 8.203461 2.977995e-10 2.995070e-10 6437.045428 6446.252596 0.000014 0.000014 1.000000e-10 1.000000e-10 1.000000e-10 1.000000e-10 4.147409e+06 1.393967e+07 3.420588e+11 1.021320e+11 0.0 0.0 0.0 0.0 16.452927 41.846974 0.0 0.0 1555 0.012378
1713 8847.390073 2.277846 1.490015 13 13 -1.0 -1.0 -1.0 0.0001 0.0001 10 8833.478917 8828.975463 1.000000e+10 1.000000e+10 0.0 0.0 0.0 0.00 0.000014 0.000014 15.923893 13.207132 4.449465e-10 3.351561e-10 7116.763618 6630.058677 0.000014 0.000014 1.000000e-10 1.000000e-10 1.000000e-10 1.000000e-10 9.656101e+05 1.016968e+07 1.314209e+12 1.248759e+11 0.0 0.0 0.0 0.0 13.911155 18.414610 0.0 0.0 1713 0.005709
2021 9341.476339 1.031129 1.260782 11 13 -1.0 -1.0 -1.0 0.0001 0.0001 10 9277.752297 9265.079542 1.000000e+10 1.000000e+10 0.0 0.0 0.0 1.44 0.007706 0.000014 0.982109 7.550195 2.870320e-05 2.721200e-10 4834.369446 6293.550535 0.007706 0.000014 1.000000e-10 1.000000e-10 1.000000e-10 1.000000e-10 6.998113e+06 1.889009e+06 0.000000e+00 6.127572e+11 0.0 0.0 0.0 0.0 63.724042 76.396796 0.0 0.0 2021 0.007076
2135 9439.155932 1.603335 1.469936 13 13 -1.0 -1.0 -1.0 0.0001 0.0001 10 9428.014059 9421.054397 1.000000e+10 1.000000e+10 0.0 0.0 0.0 0.00 0.000014 0.000014 4.711866 13.316248 3.087136e-10 2.916896e-10 6495.230174 6403.771637 0.000014 0.000014 1.000000e-10 1.000000e-10 1.000000e-10 1.000000e-10 2.258104e+06 3.675249e+06 5.262125e+11 3.235500e+11 0.0 0.0 0.0 0.0 11.141873 18.101535 0.0 0.0 2135 0.009665

We can also look at where each binary ended up, and colour the disrupted binaries differently.

[11]:

final_coords = p.get_final_mw_skycoord()

[12]:

fig, ax = plt.subplots()
ax.scatter(final_coords[:len(p)][~p.disrupted].icrs.ra, final_coords[:len(p)][~p.disrupted].icrs.dec, s=0.5, c="yellow", label="Bound")
ax.scatter(final_coords[:len(p)][p.disrupted].icrs.ra, final_coords[:len(p)][p.disrupted].icrs.dec, s=15, c="C1")
ax.scatter(final_coords[len(p):].icrs.ra, final_coords[len(p):].icrs.dec, s=15, c="C1", label="Disrupted")
ax.set(xlabel="Right Ascension [deg]", ylabel="Declination [deg]", facecolor="black")
ax.legend(markerscale=6, handletextpad=.0, framealpha=0.5)
plt.show()
../_images/pages_getting_started_24_0.png

And we see that often (but not always) they are found off the galactic plane (due to the same supernova kick that disrupted binary)

Classify population#

The Population can also classify the present-day state of each binary. You can autogenerate these classes by accessing p.classes

[13]:

p.classes

[13]:
dco co-1 co-2 xrb walkaway-t-1 walkaway-t-2 runaway-t-1 runaway-t-2 walkaway-o-1 walkaway-o-2 runaway-o-1 runaway-o-2 widow-1 widow-2 stellar-merger-co-1 stellar-merger-co-2
0 False False False False False False False False False False False False False False True False
1 False False False False False False False False False False False False False False False False
2 False False False False False False False False False False False False False False False False
3 False False False False False False False False False False False False False False False False
4 False False False False False False False False False False False False False False False False
... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ...
2637 False False False False False False False False False False False False False False False False
2638 False False False False False False False False False False False False False False False False
2639 False False False False False False False False False False False False False False False False
2640 False False False False False False False False False False False False False False False False
2641 False False False False False False False False False False False False False False False False

2641 rows × 16 columns

This table has a flag for each class and for each binary. We can do some pandas DataFrame manipulation to get some statistics for the different classes

[14]:

p.classes.astype(int).sum()

[14]:

dco                     0
co-1                    4
co-2                    4
xrb                     0
walkaway-t-1            0
walkaway-t-2            0
runaway-t-1             0
runaway-t-2             0
walkaway-o-1            0
walkaway-o-2            0
runaway-o-1             0
runaway-o-2             0
widow-1                 0
widow-2                 0
stellar-merger-co-1    28
stellar-merger-co-2     8
dtype: int64

Inspect specific binaries#

cogsworth has simple plotting routines to understand the evolution of binaries in your population. Here’s an example for a binary that experiences a disruption.

These plots show the galactic orbit of the binary (indicating supernova events and different orbits for the primary and secondary after the disruption) and a cartoon timeline of its evolution.

[15]:

random_num = np.random.choice(p.bin_nums[p.disrupted])
p.plot_orbit(random_num, t_max=50 * u.Myr)
p.plot_cartoon_binary(random_num)
../_images/pages_getting_started_31_0.png
../_images/pages_getting_started_31_1.png 
[15]:

(<Figure size 1200x1800 with 1 Axes>, <Axes: >)

Predict observables#

You can also use the Population class to predict what we might actually observe for this population in different filters. This uses dust maps to correct for extinction, applies bolometric corrections and blends the stars if necessary.

Tip

Note that you need some of the optional dependencies to run this last section - this is described more on the main Install page where you can find out how to install the optional dependencies and download the data for the dust maps and Gaia selection function.

[16]:

p.get_observables(assume_mw_galactocentric=True,
                  filters=["Gaia_G_EDR3", "Gaia_BP_EDR3", "Gaia_RP_EDR3"],
                  silence_bounds_warning=True)

[16]:
Av_1 Av_2 M_abs_1 m_app_1 M_abs_2 m_app_2 Gaia_G_EDR3_app_1 Gaia_G_EDR3_app_2 teff_obs log_g_obs secondary_brighter Gaia_G_EDR3_abs_1 Gaia_G_EDR3_abs_2 Gaia_BP_EDR3_app_1 Gaia_BP_EDR3_app_2 Gaia_BP_EDR3_abs_1 Gaia_BP_EDR3_abs_2 Gaia_RP_EDR3_app_1 Gaia_RP_EDR3_app_2 Gaia_RP_EDR3_abs_1 Gaia_RP_EDR3_abs_2
0 3.465 6.0 inf inf inf inf inf inf 6985.179305 14.252201 False inf inf inf inf inf inf inf inf inf inf
1 6.000 6.0 8.910142 23.692758 11.266718 26.049334 27.730483 inf 3635.028843 4.889894 False 9.528676 inf 30.996971 inf 10.696162 inf 26.228268 inf 8.456425 inf
2 6.000 6.0 9.192547 24.965618 9.841292 25.614363 28.571187 inf 3698.549394 4.949105 False 9.268730 inf 31.536737 inf 10.267519 inf 27.104772 inf 8.281164 inf
3 3.267 6.0 8.798432 23.290875 9.309439 23.801882 25.787157 inf 3650.238608 4.875373 False 8.979955 inf 28.117533 inf 10.136782 inf 24.416741 inf 7.914339 inf
4 4.356 6.0 7.926710 23.699495 10.033889 25.806675 26.775817 inf 3947.403375 4.772906 False 8.226711 inf 29.015353 inf 9.103774 inf 25.421915 inf 7.304179 inf
... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ...
2636 3.168 6.0 4.650803 19.140386 6.163616 20.653199 21.300855 inf 5710.276655 4.364321 False 4.328715 inf 22.527301 inf 4.694581 inf 20.222424 inf 3.793110 inf
2637 6.000 6.0 9.064636 24.120723 11.115075 26.171162 28.089145 inf 3667.325440 4.923124 False 9.559735 inf 31.198630 inf 10.637055 inf 26.604218 inf 8.529256 inf
2638 3.201 6.0 5.209269 19.633776 9.007125 23.431632 21.849265 inf 5426.556898 4.478863 False 5.104472 inf 23.109769 inf 5.496716 inf 20.760303 inf 4.542251 inf
2639 3.663 6.0 9.107925 23.489469 9.492896 23.874440 26.137370 inf 3629.701735 4.919663 False 9.207942 inf 28.584215 inf 10.340070 inf 24.746045 inf 8.155353 inf
2640 4.125 6.0 9.706837 24.787505 9.908541 24.989210 27.628021 inf 3530.806463 4.978745 False 9.775730 inf 30.336678 inf 10.966758 inf 26.197000 inf 8.701930 inf

2641 rows × 21 columns

[17]:

cogsworth.plot.plot_cmd(p)
../_images/pages_getting_started_35_0.png 
[17]:

(<Figure size 700x1000 with 2 Axes>,
 <Axes: xlabel='Gaia_BP_EDR3 - Gaia_RP_EDR3', ylabel='Gaia_G_EDR3'>)

And more…!#

That’s it for your quickstart introduction to cogsworth but that’s only the beginning of what you can do with it! You can also learn how to use hydrodynamical zoom-in simulations or predict the gravitational-wave signal in LISA for your binaries. I’ve got pages of documentation ready for you to learn all of about it.

If you’re interested in learning more you’ve got a couple of options - you could:

  • Check out the other tutorials for more step-by-step guides of using cogsworth

  • Browse the gallery of examples to see how you could put cogsworth to use

  • Search through the user guide if you want to know how a specific class or function works

Hope you enjoy using cogsworth!

Note

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