This notebook is available at https://github.com/tardis-sn/tardis/tree/master/docs/io/optional/custom_source.ipynb

Running TARDIS with a Custom Packet Source

By default, TARDIS generates energy packets using its BasePacketSource class, which models the photosphere of the supernova as a perfect blackbody (see Energy Packet Initialization). However, users may create their own packet source, as will be shown in this notebook. In order to do this, a user must create a class (that can but does not have to inherit from BasePacketSource) which contains a method create_packets that takes in (in this order): - the photospheric temperature - the number of packets - a random number generator - the photospheric radius

and returns arrays of the radii, frequencies, propogation directions, and energies of the ensemble of packets (once again see Energy Packet Initialization for more information). To use your packet source in a run of TARDIS, you must pass an instance of your class into the run_tardis function under the packet_source keyword argument.

Note

In both the BasePacketSource class and in the example here, all packets are generated at the same radius. This need not be true in general (although the create_packets method should only take in a single radius as its argument).

We show an example of how a custom packet source is used:

[1]:

# Import necessary packages
import numpy as np
from tardis import constants as const
from astropy import units as u
from tardis.montecarlo.packet_source import BasePacketSource
from tardis import run_tardis
import matplotlib.pyplot as plt
from tardis.io.atom_data import download_atom_data

/usr/share/miniconda3/envs/tardis/lib/python3.7/importlib/_bootstrap.py:219: QAWarning: pyne.data is not yet QA compliant.
  return f(*args, **kwds)
/usr/share/miniconda3/envs/tardis/lib/python3.7/site-packages/traitlets/traitlets.py:3050: FutureWarning: --rc={'figure.dpi': 96} for dict-traits is deprecated in traitlets 5.0. You can pass --rc <key=value> ... multiple times to add items to a dict.
  FutureWarning,

[2]:

# Download the atomic data used for a run of TARDIS
download_atom_data('kurucz_cd23_chianti_H_He')

[3]:

# Create a packet source class that contains a create_packets method
class TruncBlackbodySource(BasePacketSource):
    """
        Custom inner boundary source class to replace the Blackbody source
        with a truncated Blackbody source.
    """

    def __init__(self, seed, truncation_wavelength):
        super().__init__(seed)
        self.rng = np.random.default_rng(seed=seed)
        self.truncation_wavelength = truncation_wavelength

    def create_packets(self, T, no_of_packets, rng, radius,
                       drawing_sample_size=None):
        """
        Packet source that generates a truncated Blackbody source.

        Parameters
        ----------
        T : float
            Blackbody temperature
        no_of_packets : int
            number of packets to be created
        truncation_wavelength : float
            truncation wavelength in Angstrom.
            Only wavelengths higher than the truncation wavelength
            will be sampled.
        """

        # Makes uniform array of packet radii
        radii = np.ones(no_of_packets) * radius

        # Use mus and energies from normal blackbody source.
        mus = self.create_zero_limb_darkening_packet_mus(no_of_packets, self.rng)
        energies = self.create_uniform_packet_energies(no_of_packets, self.rng)

        # If not specified, draw 2 times as many packets and reject any beyond no_of_packets.
        if drawing_sample_size is None:
            drawing_sample_size = 2 * no_of_packets

        # Blackbody will be truncated below truncation_wavelength / above truncation_frequency.
        truncation_frequency = u.Quantity(self.truncation_wavelength, u.Angstrom).to(
                                          u.Hz, equivalencies=u.spectral()).value

        # Draw nus from blackbody distribution and reject based on truncation_frequency.
        # If more nus.shape[0] > no_of_packets use only the first no_of_packets.
        nus = self.create_blackbody_packet_nus(T, drawing_sample_size, self.rng)
        nus = nus[nus<truncation_frequency][:no_of_packets]


        # Only required if the truncation wavelength is too big compared to the maximum
        # of the blackbody distribution. Keep sampling until nus.shape[0] > no_of_packets.
        while nus.shape[0] < no_of_packets:
            additional_nus = self.create_blackbody_packet_nus(
                T, drawing_sample_size, self.rng
            )
            mask = additional_nus < truncation_frequency
            additional_nus = additional_nus[mask][:no_of_packets]
            nus = np.hstack([nus, additional_nus])[:no_of_packets]

        return radii, nus, mus, energies

[4]:

# Call an instance of the packet source class
packet_source = TruncBlackbodySource(
    53253, truncation_wavelength=2000
)

We now run TARDIS both with and without our custom packet source, and we compare the results:

[5]:

mdl = run_tardis('tardis_example.yml',
                 packet_source=packet_source)
mdl_norm = run_tardis('tardis_example.yml')

[py.warnings         ][WARNING]  /usr/share/miniconda3/envs/tardis/lib/python3.7/site-packages/astropy/units/equivalencies.py:124: RuntimeWarning: divide by zero encountered in double_scalars
  (si.m, si.Hz, lambda x: _si.c.value / x),
 (warnings.py:110)

Shell No.	t_rad	next_t_rad	w	next_w
0	9.93e+03	9.42e+03	0.4	0.674
5	9.85e+03	9.74e+03	0.211	0.24
10	9.78e+03	9.63e+03	0.143	0.145
15	9.71e+03	9.36e+03	0.105	0.109

[py.warnings         ][WARNING]  /usr/share/miniconda3/envs/tardis/lib/python3.7/site-packages/astropy/units/equivalencies.py:124: RuntimeWarning: divide by zero encountered in double_scalars
  (si.m, si.Hz, lambda x: _si.c.value / x),
 (warnings.py:110)

Shell No.	t_rad	next_t_rad	w	next_w
0	9.42e+03	9.98e+03	0.674	0.941
5	9.74e+03	1.03e+04	0.24	0.338
10	9.63e+03	1.02e+04	0.145	0.199
15	9.36e+03	9.96e+03	0.109	0.146

Shell No.	t_rad	next_t_rad	w	next_w
0	9.98e+03	9.34e+03	0.941	0.715
5	1.03e+04	9.74e+03	0.338	0.246
10	1.02e+04	9.63e+03	0.199	0.15
15	9.96e+03	9.39e+03	0.146	0.111

Shell No.	t_rad	next_t_rad	w	next_w
0	9.34e+03	9.89e+03	0.715	0.936
5	9.74e+03	1.02e+04	0.246	0.338
10	9.63e+03	1.01e+04	0.15	0.201
15	9.39e+03	9.87e+03	0.111	0.146

Shell No.	t_rad	next_t_rad	w	next_w
0	9.89e+03	9.38e+03	0.936	0.713
5	1.02e+04	9.63e+03	0.338	0.261
10	1.01e+04	9.5e+03	0.201	0.16
15	9.87e+03	9.32e+03	0.146	0.116

Shell No.	t_rad	next_t_rad	w	next_w
0	9.38e+03	9.88e+03	0.713	0.925
5	9.63e+03	1.02e+04	0.261	0.338
10	9.5e+03	1.01e+04	0.16	0.195
15	9.32e+03	9.96e+03	0.116	0.139

Shell No.	t_rad	next_t_rad	w	next_w
0	9.88e+03	9.41e+03	0.925	0.719
5	1.02e+04	9.71e+03	0.338	0.26
10	1.01e+04	9.6e+03	0.195	0.158
15	9.96e+03	9.35e+03	0.139	0.117

Shell No.	t_rad	next_t_rad	w	next_w
0	9.41e+03	9.89e+03	0.719	0.906
5	9.71e+03	1.02e+04	0.26	0.33
10	9.6e+03	1.01e+04	0.158	0.194
15	9.35e+03	9.88e+03	0.117	0.141

Shell No.	t_rad	next_t_rad	w	next_w
0	9.89e+03	9.45e+03	0.906	0.724
5	1.02e+04	9.74e+03	0.33	0.26
10	1.01e+04	9.63e+03	0.194	0.157
15	9.88e+03	9.45e+03	0.141	0.114

Shell No.	t_rad	next_t_rad	w	next_w
0	9.45e+03	9.86e+03	0.724	0.901
5	9.74e+03	1.02e+04	0.26	0.33
10	9.63e+03	1.01e+04	0.157	0.192
15	9.45e+03	9.84e+03	0.114	0.141

Shell No.	t_rad	next_t_rad	w	next_w
0	9.86e+03	9.43e+03	0.901	0.736
5	1.02e+04	9.72e+03	0.33	0.264
10	1.01e+04	9.64e+03	0.192	0.159
15	9.84e+03	9.4e+03	0.141	0.117

Shell No.	t_rad	next_t_rad	w	next_w
0	9.43e+03	9.89e+03	0.736	0.892
5	9.72e+03	1.02e+04	0.264	0.321
10	9.64e+03	1.01e+04	0.159	0.189
15	9.4e+03	9.88e+03	0.117	0.138

Shell No.	t_rad	next_t_rad	w	next_w
0	9.89e+03	9.48e+03	0.892	0.726
5	1.02e+04	9.77e+03	0.321	0.265
10	1.01e+04	9.68e+03	0.189	0.16
15	9.88e+03	9.42e+03	0.138	0.12

Shell No.	t_rad	next_t_rad	w	next_w
0	9.48e+03	9.83e+03	0.726	0.871
5	9.77e+03	1.01e+04	0.265	0.321
10	9.68e+03	1.01e+04	0.16	0.186
15	9.42e+03	9.82e+03	0.12	0.136

Shell No.	t_rad	next_t_rad	w	next_w
0	9.83e+03	9.51e+03	0.871	0.754
5	1.01e+04	9.77e+03	0.321	0.275
10	1.01e+04	9.71e+03	0.186	0.163
15	9.82e+03	9.51e+03	0.136	0.118

Shell No.	t_rad	next_t_rad	w	next_w
0	9.51e+03	9.79e+03	0.754	0.873
5	9.77e+03	1e+04	0.275	0.319
10	9.71e+03	9.98e+03	0.163	0.189
15	9.51e+03	9.8e+03	0.118	0.136

Shell No.	t_rad	next_t_rad	w	next_w
0	9.79e+03	9.54e+03	0.873	0.752
5	1e+04	9.81e+03	0.319	0.273
10	9.98e+03	9.74e+03	0.189	0.161
15	9.8e+03	9.56e+03	0.136	0.117

Shell No.	t_rad	next_t_rad	w	next_w
0	9.54e+03	9.8e+03	0.752	0.865
5	9.81e+03	1.01e+04	0.273	0.312
10	9.74e+03	1e+04	0.161	0.185
15	9.56e+03	9.83e+03	0.117	0.134

Shell No.	t_rad	next_t_rad	w	next_w
0	9.8e+03	9.51e+03	0.865	0.761
5	1.01e+04	9.78e+03	0.312	0.276
10	1e+04	9.73e+03	0.185	0.163
15	9.83e+03	9.53e+03	0.134	0.119

Shell No.	t_rad	next_t_rad	w	next_w
0	9.93e+03	1.01e+04	0.4	0.507
5	9.85e+03	1.02e+04	0.211	0.197
10	9.78e+03	1.01e+04	0.143	0.117
15	9.71e+03	9.87e+03	0.105	0.0869

Shell No.	t_rad	next_t_rad	w	next_w
0	1.01e+04	1.15e+04	0.507	0.546
5	1.02e+04	1.15e+04	0.197	0.223
10	1.01e+04	1.13e+04	0.117	0.135
15	9.87e+03	1.1e+04	0.0869	0.101

Shell No.	t_rad	next_t_rad	w	next_w
0	1.15e+04	1.05e+04	0.546	0.442
5	1.15e+04	1.08e+04	0.223	0.165
10	1.13e+04	1.07e+04	0.135	0.101
15	1.1e+04	1.02e+04	0.101	0.0781

Shell No.	t_rad	next_t_rad	w	next_w
0	1.05e+04	1.16e+04	0.442	0.498
5	1.08e+04	1.17e+04	0.165	0.201
10	1.07e+04	1.14e+04	0.101	0.127
15	1.02e+04	1.1e+04	0.0781	0.0953

Shell No.	t_rad	next_t_rad	w	next_w
0	1.16e+04	1.06e+04	0.498	0.444
5	1.17e+04	1.09e+04	0.201	0.162
10	1.14e+04	1.07e+04	0.127	0.101
15	1.1e+04	1.03e+04	0.0953	0.0769

Shell No.	t_rad	next_t_rad	w	next_w
0	1.06e+04	1.15e+04	0.444	0.51
5	1.09e+04	1.17e+04	0.162	0.201
10	1.07e+04	1.14e+04	0.101	0.127
15	1.03e+04	1.1e+04	0.0769	0.0955

Shell No.	t_rad	next_t_rad	w	next_w
0	1.15e+04	1.06e+04	0.51	0.443
5	1.17e+04	1.1e+04	0.201	0.158
10	1.14e+04	1.07e+04	0.127	0.102
15	1.1e+04	1.03e+04	0.0955	0.0768

Shell No.	t_rad	next_t_rad	w	next_w
0	1.06e+04	1.15e+04	0.443	0.502
5	1.1e+04	1.17e+04	0.158	0.197
10	1.07e+04	1.14e+04	0.102	0.124
15	1.03e+04	1.1e+04	0.0768	0.0941

Shell No.	t_rad	next_t_rad	w	next_w
0	1.15e+04	1.06e+04	0.502	0.443
5	1.17e+04	1.1e+04	0.197	0.161
10	1.14e+04	1.07e+04	0.124	0.102
15	1.1e+04	1.04e+04	0.0941	0.0753

Shell No.	t_rad	next_t_rad	w	next_w
0	1.06e+04	1.15e+04	0.443	0.505
5	1.1e+04	1.17e+04	0.161	0.197
10	1.07e+04	1.15e+04	0.102	0.122
15	1.04e+04	1.11e+04	0.0753	0.0933

Shell No.	t_rad	next_t_rad	w	next_w
0	1.15e+04	1.06e+04	0.505	0.444
5	1.17e+04	1.1e+04	0.197	0.162
10	1.15e+04	1.09e+04	0.122	0.0958
15	1.11e+04	1.05e+04	0.0933	0.0744

Shell No.	t_rad	next_t_rad	w	next_w
0	1.06e+04	1.16e+04	0.444	0.486
5	1.1e+04	1.18e+04	0.162	0.19
10	1.09e+04	1.14e+04	0.0958	0.121
15	1.05e+04	1.11e+04	0.0744	0.089

Shell No.	t_rad	next_t_rad	w	next_w
0	1.16e+04	1.07e+04	0.486	0.438
5	1.18e+04	1.1e+04	0.19	0.162
10	1.14e+04	1.08e+04	0.121	0.101
15	1.11e+04	1.04e+04	0.089	0.0767

Shell No.	t_rad	next_t_rad	w	next_w
0	1.07e+04	1.15e+04	0.438	0.495
5	1.1e+04	1.17e+04	0.162	0.191
10	1.08e+04	1.14e+04	0.101	0.12
15	1.04e+04	1.1e+04	0.0767	0.0918

Shell No.	t_rad	next_t_rad	w	next_w
0	1.15e+04	1.07e+04	0.495	0.449
5	1.17e+04	1.1e+04	0.191	0.163
10	1.14e+04	1.09e+04	0.12	0.0997
15	1.1e+04	1.05e+04	0.0918	0.0768

Shell No.	t_rad	next_t_rad	w	next_w
0	1.07e+04	1.15e+04	0.449	0.495
5	1.1e+04	1.18e+04	0.163	0.184
10	1.09e+04	1.14e+04	0.0997	0.122
15	1.05e+04	1.11e+04	0.0768	0.0898

Shell No.	t_rad	next_t_rad	w	next_w
0	1.15e+04	1.07e+04	0.495	0.439
5	1.18e+04	1.11e+04	0.184	0.16
10	1.14e+04	1.1e+04	0.122	0.0962
15	1.11e+04	1.05e+04	0.0898	0.0763

Shell No.	t_rad	next_t_rad	w	next_w
0	1.07e+04	1.14e+04	0.439	0.492
5	1.11e+04	1.17e+04	0.16	0.186
10	1.1e+04	1.15e+04	0.0962	0.117
15	1.05e+04	1.11e+04	0.0763	0.0876

Shell No.	t_rad	next_t_rad	w	next_w
0	1.14e+04	1.07e+04	0.492	0.441
5	1.17e+04	1.11e+04	0.186	0.161
10	1.15e+04	1.09e+04	0.117	0.1
15	1.11e+04	1.05e+04	0.0876	0.0765

[6]:

%matplotlib inline
plt.plot(mdl.runner.spectrum_virtual.wavelength,
         mdl.runner.spectrum_virtual.luminosity_density_lambda,
         color='red', label='truncated blackbody (custom packet source)')
plt.plot(mdl_norm.runner.spectrum_virtual.wavelength,
         mdl_norm.runner.spectrum_virtual.luminosity_density_lambda,
         color='blue', label='normal blackbody (default packet source)')
plt.xlabel('$\lambda [\AA]$')
plt.ylabel('$L_\lambda$ [erg/s/$\AA$]')
plt.xlim(500, 10000)
plt.legend()

[6]:

<matplotlib.legend.Legend at 0x7f2041fc1cd0>

../../_images/io_optional_custom_source_8_1.svg