Tutorial: Basic functionality

Thank you for your interest in using kilopop.

This is a tutorial showing the basic functionality of how to use the module and produce kilonovae lightcurves.

For an interactive version of the same material, see the Jupyter notebook in the tutorial folder of the source repository.

%matplotlib notebook
import matplotlib.pyplot as plt

from kilopop.kilonovae import bns_kilonova as saeev
from kilopop.kilonovae import bns_kilonovae_population_distribution as s22p
import sncosmo

In order to create an instance of the model, we actually do not need to specify any inputs. In that case, all parameters will be drawn from the distributions used in the accompanying paper.

test_inst = saeev()

Now we can see the values that each of the parameters of the model have taken.

# List out the parameters of the BNS mergers kilonova and binary inspiral
test_inst.print_parameters()

Parameter 1: mass1: 1.35
Parameter 2: mass2: 1.35
Parameter 3: compactness1: 0.16726128755607012
Parameter 4: compactness2: 0.16726128755607012
Parameter 5: viewing_angle: 0.8711099188079074
Parameter 6: electron_fraction: 0.25562517175591914
Parameter 7: dynamical_ejecta_mass: 0.0023929897684419003
Parameter 8: median_ejecta_velocity: 0.2332274868584348
Parameter 9: grey_opacity: 2.1586924243196335
Parameter 10: secular_ejecta_mass: 0.023624468109337757
Parameter 11: total_ejecta_mass: 0.026017457877779655
Parameter 12: disk_unbinding_efficiency: 0.13436192902029392

Of course, we can always specify any of these parameters ourselves. For example:

test_inst = saeev(mass1=1.35, mass2=1.35)

# List out the parameters of the BNS mergers kilonova and binary inspiral
test_inst.print_parameters()

Parameter 1: mass1: 1.35
Parameter 2: mass2: 1.35
Parameter 3: compactness1: 0.16726128755607012
Parameter 4: compactness2: 0.16726128755607012
Parameter 5: viewing_angle: 0.8711099188079074
Parameter 6: electron_fraction: 0.25562517175591914
Parameter 7: dynamical_ejecta_mass: 0.0023929897684419003
Parameter 8: median_ejecta_velocity: 0.2332274868584348
Parameter 9: grey_opacity: 2.1586924243196335
Parameter 10: secular_ejecta_mass: 0.023624468109337757
Parameter 11: total_ejecta_mass: 0.026017457877779655
Parameter 12: disk_unbinding_efficiency: 0.13436192902029392

We can see that each intstance of the kilonova has different attributes:

test_inst.__dict__.keys()

dict_keys(['id', 'number_of_parameters', 'min_wave', 'max_wave', 'transient_duration', 'mass_ratio_threshold', 'EOS_path', 'emulator_path', 'opacity_data_path', 'param1_name', 'param2_name', 'param3_name', 'param4_name', 'param5_name', 'param6_name', 'param7_name', 'param8_name', 'param9_name', 'param10_name', 'param11_name', 'param12_name', 'param1', 'param2', 'param3', 'param4', 'param5', 'param6', 'param7', 'param8', 'param9', 'param10', 'param11', 'param12', 'phase', 'wave', 'flux', 'model', 'redshift', 'dist_mpc', 't0', 'tmax'])

Underneath this interface we have used the framework of sncosmo.model to handle the SED-timeseries for each realisation. So you can use most of the functionality of sncosmo with the model realisations that are generated. NB: As the model itself is quite complicated and not itself integrated into sncosmo it is not possible to change the parameters of a given realisation using the usual model.set(mass1=x) API, for example.

Note sncosmo naturally plots lightcurves only on one day increments. So the early time behaviour is low resolution in these figures, but the model is more finely detailed.

fig = sncosmo.plot_lc(model=test_inst.model, bands=['lsstg', 'lsstr', 'lssti','lsstz'])
plt.show()

We can also generate the distributions of parameters for a population of kilonovae. If we want to also compute the lightcurve properties, this can be a computationally intensive task for a personal computer or laptop. That being said, in this example we will only generate a population of 1000 kilonovae, as opposed to 50000 that was used in the paper. NB: the flag only_draw_parameters specifies whether or not the user only wants the parameter distributions or the parameter distributions and the derived light curve properties (which is more computationally intensive).

test_dist = s22p(population_size=1000, only_draw_parameters=False, chunk_size=100)

  0%|                                                                                             | 0/1000 [00:00<?, ?it/s]/Users/cnsetzer/anaconda3/envs/kne_dev/lib/python3.9/site-packages/sncosmo/models.py:189: RuntimeWarning: invalid value encountered in log10
  result[i] = -2.5 * np.log10(f / zpf)
/Users/cnsetzer/anaconda3/envs/kne_dev/lib/python3.9/site-packages/sncosmo/models.py:189: RuntimeWarning: invalid value encountered in log10
  result[i] = -2.5 * np.log10(f / zpf)
/Users/cnsetzer/anaconda3/envs/kne_dev/lib/python3.9/site-packages/sncosmo/models.py:189: RuntimeWarning: invalid value encountered in log10
  result[i] = -2.5 * np.log10(f / zpf)
/Users/cnsetzer/anaconda3/envs/kne_dev/lib/python3.9/site-packages/sncosmo/models.py:189: RuntimeWarning: invalid value encountered in log10
  result[i] = -2.5 * np.log10(f / zpf)
100%|██████████████████████████████████████████████████████████████████████████████████| 1000/1000 [07:03<00:00,  2.36it/s]

With this distribution generated, we can plot the parameter and lightcurve property distributions.

for i in range(12):
    plt.figure()
    plt.hist(getattr(test_dist, f'param{i+1}'))
    plt.xlabel(getattr(test_dist, f'param{i+1}_name'))
    plt.show()

plt.figure()
plt.hist(test_dist.peak_absmag_lssti, bins=20, density=True)
plt.xlabel(r'$M_\mathrm{abs, {\it lssti}}$')
plt.show()

plt.figure()
plt.hist(test_dist.one_mag_peak_time_lssti, bins=20, density=True)
plt.xlabel(r'Time within 1-mag. of $M_\mathrm{abs, {\it lssti}}$')
plt.xlim(0,5)
plt.show()

Note, a rare few high-opacity kilonovae have nearly flat evolution which can be seen above as those that spend significantly long amounts of time, >~5 days, within 1 magnitude of the peak mag. These will all be unobservable unless exceptionally nearby.

We can also make scatter plots like in the paper of the distributed parameters.

# ejecta parameter scatterplot
plt.figure()
plt.scatter(test_dist.param11, test_dist.param8)
plt.xlabel(r'$m_\mathrm{ej,total}$ [$M_\odot$]')
plt.ylabel(r'$v_\mathrm{ej}$ [$c$]')
plt.show()