You can interact with this notebook online: Launch notebook

Gamma-ray energy deposition

This notebook provides the initial implementation of Gamma-ray energy deposition into an arbitrary ejecta. It is a WORK IN PROGRESS and should NOT be used for any science work until further notice.

Main loop

Generates a simple 1D ejecta and a list of gamma-ray objects.

Runs packets of gamma-rays through the ejecta. Handles interactions by calling the appropriate function.

Adds deposited energy and output energy to 2 different dataframes.

Model setup

[4]:

import numpy as np

import matplotlib.pyplot as plt
from tardis.model import SimulationState
from tardis.io.config_reader import Configuration

from tardis.plasma.properties import Density, Abundance, IsotopeAbundance, IsotopeNumberDensity, AtomicData, AtomicMass, IsotopeMass, NumberDensity, SelectedAtoms
from tardis.plasma.base import BasePlasma
from tardis.io.atom_data import AtomData

import astropy.units as u
from tardis.energy_input.gamma_ray_transport import main_gamma_ray_loop

# Adjust model
config = Configuration.from_yaml("../../../../tardis/io/tests/data/tardis_configv1_density_exponential_nebular.yml")
config.model.structure.velocity.start = 1 * u.km / u.s
config.model.structure.density.rho_0 = 5e2 * u.g / (u.cm ** 3)
config.supernova.time_explosion = 2.0 * u.d

config.atom_data = "kurucz_cd23_chianti_H_He.h5"

model = SimulationState.from_config(config)

# Construct plasma
input = [Density, Abundance, IsotopeAbundance, AtomicData, AtomicMass, IsotopeNumberDensity, NumberDensity, SelectedAtoms, IsotopeMass]

plasma = BasePlasma(
        plasma_properties=input,
        density = model.density,
        abundance = model.abundance,
        isotope_abundance = model.raw_isotope_abundance,
        atomic_data = AtomData.from_hdf(config.atom_data)
    )

# Set up packet count
num_packets = 500000

# Lock seed
np.random.seed(1)

Generate plasma

[5]:


from tardis.plasma.properties import Density, Abundance, IsotopeAbundance, IsotopeNumberDensity, AtomicData, AtomicMass, IsotopeMass, NumberDensity, SelectedAtoms
from tardis.plasma.base import BasePlasma
from tardis.io.atom_data import AtomData

input = [Density, Abundance, IsotopeAbundance, AtomicData, AtomicMass, IsotopeNumberDensity, NumberDensity, SelectedAtoms, IsotopeMass]

plasma = BasePlasma(
        plasma_properties=input,
        density = model.density,
        abundance = model.abundance,
        isotope_abundance = model.raw_isotope_abundance,
        atomic_data = AtomData.from_hdf(config.atom_data)
    )

[7]:


# Compute energy deposition rate
# ejecta_energy_df is the deposited energy
# ejecta_plot_energy_df is information for plotting
# escape_energy is the escaping energy
# decayed_packet_count is the number of packets created per shell
# energy_plot_positrons is the deposited energy from positrons
# estimated_deposition is the deposited energy from the Kasen (2006) estimator (currently not functional)
(
    energy_df,
    energy_plot_df,
    escape_energy,
    decayed_packet_count,
    energy_plot_positrons,
    estimated_deposition
) = main_gamma_ray_loop(
    num_packets,
    model,
    plasma,
    time_steps=50,
    time_end=50.0,
    path_to_decay_data="~/compiled_ensdf.hdf"
)

ejecta_energy = energy_plot_df["energy_input"]
ejecta_energy_r = energy_plot_df["energy_input_r"]
energy_input_time = energy_plot_df["energy_input_time"]
energy_input_type = energy_plot_df["energy_input_type"]

Total gamma-ray energy
6.389259861704378e+49
Total positron energy
8.454621397870039e+47
Total packets: 500000
Energy per packet 7.975724746893494e+52
Initializing packets
Total positron energy from packets
8.417997803471385e+47
Total CMF energy
3.7264874091697995e+58
Total RF energy
3.7279788130927574e+58

/home/afullard/tardis/tardis/energy_input/gamma_packet_loop.py:129: NumbaPerformanceWarning: np.dot() is faster on contiguous arrays, called on (array(float64, 1d, A), array(float64, 1d, C))
  doppler_factor = doppler_factor_3d(
/home/afullard/tardis/tardis/energy_input/gamma_packet_loop.py:197: NumbaPerformanceWarning: np.dot() is faster on contiguous arrays, called on (array(float64, 1d, C), array(float64, 1d, A))
  ) = distance_trace(
/home/afullard/tardis/tardis/energy_input/gamma_packet_loop.py:245: NumbaPerformanceWarning: np.dot() is faster on contiguous arrays, called on (array(float64, 1d, A), array(float64, 1d, C))
  packet, ejecta_energy_gained = process_packet_path(packet)

Entering gamma ray loop for 500000 packets
Escaped packets: 81593
Scattered packets: 23333
Final energy to test for conservation
3.720648730441454e+58

Plotting results

Energy deposited at a given radius

[8]:

fig = plt.figure(dpi=150, facecolor='w')
ax = fig.add_subplot(111)

scatter = ax.scatter(np.array(ejecta_energy_r)/np.max(model.v_outer.value), np.array(ejecta_energy), s=1, alpha=0.1)
ax.set_xlabel("r/R")
ax.set_ylabel("E (keV)")
ax.semilogy();
../../../_images/contributing_in_progress_nebular_phase_gammaray_deposition_8_0.png

Interactions binned by radius

[9]:

fig = plt.figure(dpi=150, facecolor='w')
ax = fig.add_subplot(111)
ax.hist(np.array(ejecta_energy_r), bins=200)
#ax.set_xlim(0, 1)
ax.set_xlabel("r (cm)")
ax.set_ylabel("Interaction count");
../../../_images/contributing_in_progress_nebular_phase_gammaray_deposition_10_0.png

Density Profile

[10]:

fig = plt.figure(dpi=150, facecolor='w')
plt.semilogy(model.r_middle/np.max(model.r_outer), model.density)
plt.plot(0,0)
plt.ylabel("Density g cm$^{-3}$")
plt.xlabel("r/R");
../../../_images/contributing_in_progress_nebular_phase_gammaray_deposition_12_0.png

Fraction of energy escaping from the ejecta

[11]:

np.sum(np.sum(escape_energy) / (np.sum(escape_energy) + np.sum(energy_df)))

[11]:

1.844372140556195e-53

Spectrum of escape energy at the final time step

[12]:

from tardis.energy_input.energy_source import read_artis_lines

ni56_lines = read_artis_lines("ni56", "~/Downloads/tardisnuclear/")
co56_lines = read_artis_lines("co56", "~/Downloads/tardisnuclear/")

plt.figure(figsize=(12, 6), dpi=150)
plt.step(escape_energy.index, escape_energy.iloc[:,49], label="$\gamma$-spectrum", where="post")
plt.xlabel("Energy (keV)")
plt.ylabel("keV/keV/s");
plt.loglog()
plt.xlim(100, 4000)
#plt.ylim(0.01, 100)
plt.vlines(ni56_lines.energy*1000, 0.1, 0.2, color="k", label="Ni56")
plt.vlines(co56_lines.energy*1000, 0.1, 0.2, color="r", label="Co56")
plt.legend()

/home/afullard/tardis/tardis/energy_input/energy_source.py:170: ParserWarning: Falling back to the 'python' engine because the 'c' engine does not support regex separators (separators > 1 char and different from '\s+' are interpreted as regex); you can avoid this warning by specifying engine='python'.
  return pd.read_csv(

[12]:

<matplotlib.legend.Legend at 0x7eff9e39b8e0>
../../../_images/contributing_in_progress_nebular_phase_gammaray_deposition_16_2.png

Energy deposition rate

Dataframe index is the radial grid location. Columns are time steps in seconds

[13]:

energy_df

[13]:
1.728000e+05 1.842903e+05 1.965447e+05 2.096139e+05 2.235522e+05 2.384173e+05 2.542708e+05 2.711786e+05 2.892106e+05 3.084416e+05 ... 2.269320e+06 2.420218e+06 2.581151e+06 2.752784e+06 2.935830e+06 3.131048e+06 3.339247e+06 3.561290e+06 3.798098e+06 4.050652e+06
0 1.268952e+54 1.434033e+53 1.283598e+53 1.493515e+53 1.665060e+53 1.864771e+53 1.369552e+53 1.291566e+53 1.206870e+53 1.210695e+53 ... 2.413105e+52 1.946078e+52 1.601020e+52 1.681973e+52 1.532812e+52 1.746724e+52 1.142765e+52 1.155513e+52 1.013063e+52 6.851335e+51
1 6.352248e+54 6.256791e+53 5.862575e+53 6.364949e+53 7.788701e+53 7.056560e+53 6.954693e+53 7.225777e+53 6.150372e+53 6.159685e+53 ... 1.326250e+53 1.128764e+53 9.533177e+52 9.066819e+52 8.349490e+52 6.819232e+52 7.361927e+52 6.937388e+52 5.521254e+52 4.563154e+52
2 1.256390e+55 1.224663e+54 1.455274e+54 1.393148e+54 1.294472e+54 1.200354e+54 1.191072e+54 1.325274e+54 1.247067e+54 1.228414e+54 ... 2.548472e+53 2.394402e+53 2.271846e+53 1.892756e+53 1.689675e+53 1.628341e+53 1.372483e+53 1.343997e+53 1.264984e+53 1.023617e+53
3 1.720830e+55 1.969217e+54 1.770751e+54 1.902969e+54 1.772499e+54 1.860536e+54 1.849413e+54 1.688380e+54 1.837876e+54 1.673101e+54 ... 3.721530e+53 3.252044e+53 2.997976e+53 2.889827e+53 2.407073e+53 2.254808e+53 2.174352e+53 1.838789e+53 1.819876e+53 1.590544e+53
4 2.061597e+55 2.118889e+54 2.361284e+54 2.224901e+54 1.929322e+54 2.138034e+54 2.121064e+54 1.924094e+54 1.987540e+54 1.973124e+54 ... 4.067807e+53 3.866050e+53 3.417029e+53 3.527510e+53 2.975597e+53 2.717593e+53 2.501916e+53 2.181535e+53 1.990705e+53 1.963625e+53
5 2.183037e+55 2.313466e+54 2.486110e+54 2.332249e+54 2.466149e+54 2.202191e+54 2.202279e+54 2.211636e+54 2.227595e+54 2.093595e+54 ... 4.533677e+53 4.126323e+53 3.522390e+53 3.371528e+53 3.012767e+53 2.863011e+53 2.678039e+53 2.315553e+53 1.989705e+53 1.950900e+53
6 2.314019e+55 2.541623e+54 2.458594e+54 2.224427e+54 2.101068e+54 2.393609e+54 2.401742e+54 2.069722e+54 2.161520e+54 2.233445e+54 ... 4.401336e+53 3.829897e+53 3.344355e+53 3.361478e+53 2.807300e+53 2.594261e+53 2.613329e+53 2.173414e+53 2.075456e+53 1.861470e+53
7 2.111323e+55 2.151525e+54 2.260010e+54 2.081330e+54 2.245851e+54 2.237418e+54 2.117185e+54 2.193539e+54 2.033406e+54 1.965041e+54 ... 4.041338e+53 3.574777e+53 3.243658e+53 2.860518e+53 2.836507e+53 2.387642e+53 2.345213e+53 2.193704e+53 1.705966e+53 1.731529e+53
8 1.922621e+55 2.195889e+54 1.937526e+54 1.994053e+54 2.080152e+54 2.091939e+54 2.013608e+54 1.954220e+54 1.887144e+54 1.888389e+54 ... 3.376278e+53 3.379366e+53 2.893652e+53 2.460482e+53 2.349191e+53 2.123134e+53 1.840041e+53 1.691727e+53 1.431448e+53 1.261554e+53
9 1.725510e+55 1.852217e+54 1.755411e+54 2.035626e+54 1.885113e+54 1.804745e+54 1.800571e+54 1.683805e+54 1.719869e+54 1.687133e+54 ... 2.819492e+53 2.637785e+53 2.536277e+53 2.043085e+53 1.644105e+53 1.596846e+53 1.468581e+53 1.328962e+53 1.123189e+53 9.315640e+52
10 1.478802e+55 1.811474e+54 1.466713e+54 1.580566e+54 1.612419e+54 1.622183e+54 1.479221e+54 1.564288e+54 1.383014e+54 1.467392e+54 ... 2.130141e+53 2.118125e+53 1.711647e+53 1.648757e+53 1.546999e+53 1.292042e+53 1.012572e+53 1.090591e+53 8.554179e+52 7.878659e+52
11 1.308745e+55 1.460406e+54 1.447386e+54 1.346796e+54 1.402081e+54 1.420153e+54 1.299567e+54 1.205790e+54 1.428299e+54 1.339745e+54 ... 1.757030e+53 1.548390e+53 1.376493e+53 1.125006e+53 1.040292e+53 9.341924e+52 8.638519e+52 7.397325e+52 6.592805e+52 5.976537e+52
12 1.123812e+55 1.225134e+54 1.087325e+54 1.129301e+54 1.278660e+54 1.068451e+54 1.077110e+54 1.102109e+54 1.066945e+54 1.003824e+54 ... 1.138373e+53 9.911889e+52 8.863272e+52 9.083607e+52 7.878527e+52 6.366229e+52 6.144144e+52 4.775543e+52 4.849689e+52 4.762654e+52
13 9.236778e+54 1.055157e+54 9.396187e+53 8.982871e+53 1.015297e+54 8.711438e+53 9.211775e+53 8.318177e+53 8.380983e+53 8.134140e+53 ... 9.137514e+52 8.480772e+52 7.327633e+52 5.763397e+52 5.614848e+52 4.220384e+52 4.918503e+52 3.944850e+52 3.457032e+52 2.885190e+52
14 7.707517e+54 8.334156e+53 9.577188e+53 7.962057e+53 6.930276e+53 8.204323e+53 7.828914e+53 7.164728e+53 6.386090e+53 6.956695e+53 ... 5.413666e+52 4.705275e+52 4.130852e+52 3.872048e+52 4.527747e+52 3.207596e+52 2.490193e+52 2.495629e+52 2.457263e+52 1.981165e+52
15 6.191768e+54 5.473081e+53 6.655216e+53 6.239894e+53 6.975514e+53 6.084465e+53 7.786685e+53 5.483303e+53 5.721475e+53 5.412126e+53 ... 3.969469e+52 4.340826e+52 3.682750e+52 3.111694e+52 2.165591e+52 2.104476e+52 2.120679e+52 2.012698e+52 1.412540e+52 1.609626e+52
16 4.915662e+54 4.883089e+53 5.309932e+53 5.435722e+53 5.847968e+53 5.329650e+53 4.808035e+53 3.888501e+53 5.008029e+53 4.654530e+53 ... 2.642693e+52 2.896227e+52 2.212667e+52 2.213041e+52 1.790653e+52 1.690462e+52 1.916662e+52 1.402938e+52 8.857395e+51 1.167084e+52
17 4.147753e+54 4.613959e+53 3.960488e+53 4.516572e+53 4.716706e+53 4.121339e+53 4.278267e+53 4.491534e+53 2.649101e+53 3.225867e+53 ... 2.639209e+52 1.600225e+52 1.834074e+52 1.458769e+52 1.407996e+52 1.274043e+52 9.493741e+51 8.616331e+51 8.732235e+51 8.816470e+51
18 3.245799e+54 3.905090e+53 3.410483e+53 3.592481e+53 2.625412e+53 3.660111e+53 3.573277e+53 4.006494e+53 2.275214e+53 2.788264e+53 ... 1.712318e+52 7.255677e+51 1.114020e+52 1.186242e+52 1.025083e+52 9.771694e+51 6.949339e+51 7.629940e+51 6.335910e+51 4.726639e+51
19 2.490859e+54 3.183146e+53 1.705482e+53 2.854045e+53 2.615611e+53 2.653951e+53 2.443372e+53 2.982978e+53 2.598039e+53 2.129802e+53 ... 9.081512e+51 5.451095e+51 9.925498e+51 8.614093e+51 9.119476e+51 5.419938e+51 4.528037e+51 3.783831e+51 4.377392e+51 3.675127e+51

20 rows × 50 columns

Energy deposition rate versus radius at time_explosion

[14]:

fig = plt.figure(dpi=150, facecolor='w')
plt.semilogy(model.v_middle/np.max(model.v_outer), energy_df.iloc[:, 0])
plt.xlabel("r/R")
plt.ylabel("Energy deposition rate [eV/s$]");
../../../_images/contributing_in_progress_nebular_phase_gammaray_deposition_20_0.png

Energy Deposition rate versus time compared to the analytic estimate

[15]:

energy = energy_df.T.sum(axis=1) * u.eV.to('erg')

def analytic_estimate(t):
    return 0.6 * (0.97 * (1 - np.exp(-(40 / t)**2)) + 0.03) * (6.45 * np.exp(-t/8.8) + 1.45 * np.exp(-t/111.3)) * 1e43

t_tardis = energy.index * u.s.to('d')

plt.figure(dpi=150)
plt.plot(t_tardis, energy / analytic_estimate(t_tardis), "-")
plt.xlabel("t_exp (d)")
plt.ylabel("Energy deposition / analytic")
plt.legend()
plt.xlim(0, 50)
plt.ylim(0.7, 1.1)

No artists with labels found to put in legend.  Note that artists whose label start with an underscore are ignored when legend() is called with no argument.

[15]:

(0.7, 1.1)
../../../_images/contributing_in_progress_nebular_phase_gammaray_deposition_22_2.png 
[ ]: