You can interact with this notebook online: Launch notebook

Analysing Spectral Energy Distribution and Emission/Absorption Visualization

Note:

This notebook is only a sample demonstrating some of the features of the sdecplotter class. If you are interested in using additional features, you should directly access the sdecplotter class. You can see the rest of the features of the sdecplotter class here.

A notebook for analyzing and visualizing the spectral energy distribution, emission and absorption patterns in supernova simulations using TARDIS.

Every simulation run requires atomic data and a configuration file.

Atomic Data

We recommend using the kurucz_cd23_chianti_H_He.h5 dataset.

from tardis.io.atom_data import download_atom_data

# We download the atomic data needed to run the simulation
download_atom_data('kurucz_cd23_chianti_H_He')

[tardis.io.atom_data.atom_web_download][WARNING]
        Atomic Data kurucz_cd23_chianti_H_He already exists in /home/abhinav/Downloads/tardis-data/kurucz_cd23_chianti_H_He.h5. Will not download - override with force_download=True. (atom_web_download.py:50)

Example Configuration File

!wget -q -nc https://raw.githubusercontent.com/tardis-sn/tardis/master/docs/tardis_example.yml

!cat tardis_example.yml

# Example YAML configuration for TARDIS
tardis_config_version: v1.0

supernova:
  luminosity_requested: 9.44 log_lsun
  time_explosion: 13 day

atom_data: kurucz_cd23_chianti_H_He.h5

model:
  structure:
    type: specific
    velocity:
      start: 1.1e4 km/s
      stop: 20000 km/s
      num: 20
    density:
      type: branch85_w7

  abundances:
    type: uniform
    O: 0.19
    Mg: 0.03
    Si: 0.52
    S: 0.19
    Ar: 0.04
    Ca: 0.03

plasma:
  disable_electron_scattering: no
  ionization: lte
  excitation: lte
  radiative_rates_type: dilute-blackbody
  line_interaction_type: macroatom

montecarlo:
  seed: 23111963
  no_of_packets: 4.0e+4
  iterations: 20
  nthreads: 1

  last_no_of_packets: 1.e+5
  no_of_virtual_packets: 10

  convergence_strategy:
    type: damped
    damping_constant: 1.0
    threshold: 0.05
    fraction: 0.8
    hold_iterations: 3
    t_inner:
      damping_constant: 0.5

spectrum:
  start: 500 angstrom
  stop: 20000 angstrom
  num: 10000

Loading Simulation Data

Running simulation

To run the simulation, import the run_tardis function and create the sim object.

Note:

Get more information about the progress bars, logging configuration, and convergence plots.

from tardis import run_tardis

simulation = run_tardis(
    "tardis_example.yml",
    virtual_packet_logging=True,
    show_convergence_plots=True,
    export_convergence_plots=True,
    log_level="INFO",
)

[tardis.io.model.parse_atom_data][INFO   ]

        Reading Atomic Data from kurucz_cd23_chianti_H_He.h5 (parse_atom_data.py:40)
[tardis.io.atom_data.util][INFO   ]

        Atom Data kurucz_cd23_chianti_H_He.h5 not found in local path.
        Exists in TARDIS Data repo /home/abhinav/Downloads/tardis-data/kurucz_cd23_chianti_H_He.h5 (util.py:34)
[tardis.io.atom_data.base][INFO   ]
        Reading Atom Data with: UUID = 6f7b09e887a311e7a06b246e96350010 MD5  = 864f1753714343c41f99cb065710cace  (base.py:258)
[tardis.io.atom_data.base][INFO   ]
        Non provided Atomic Data: synpp_refs, photoionization_data, yg_data, two_photon_data, linelist_atoms, linelist_molecules (base.py:262)
[tardis.io.model.parse_density_configuration][WARNING]
        Number of density points larger than number of shells. Assuming inner point irrelevant (parse_density_configuration.py:114)
[tardis.model.matter.decay][INFO   ]
        Decaying abundances for 1123200.0 seconds (decay.py:101)
[tardis.simulation.base][INFO   ]

        Starting iteration 1 of 20 (base.py:444)

[tardis.simulation.base][INFO   ]

        Luminosity emitted   = 7.942e+42 erg / s
        Luminosity absorbed  = 2.659e+42 erg / s
        Luminosity requested = 1.059e+43 erg / s
 (base.py:657)
[tardis.simulation.base][INFO   ]

        Plasma stratification: (base.py:625)

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

[tardis.simulation.base][INFO   ]

        Current t_inner = 9933.952 K
        Expected t_inner for next iteration = 10703.212 K
 (base.py:652)

[tardis.simulation.base][INFO   ]

        Starting iteration 2 of 20 (base.py:444)
[tardis.simulation.base][INFO   ]

        Luminosity emitted   = 1.071e+43 erg / s
        Luminosity absorbed  = 3.576e+42 erg / s
        Luminosity requested = 1.059e+43 erg / s
 (base.py:657)
[tardis.simulation.base][INFO   ]

        Plasma stratification: (base.py:625)

Shell No.	t_rad	next_t_rad	w	next_w
0	1.01e+04 K	1.08e+04 K	0.507	0.525
5	1.02e+04 K	1.1e+04 K	0.197	0.203
10	1.01e+04 K	1.08e+04 K	0.117	0.125
15	9.87e+03 K	1.05e+04 K	0.0869	0.0933

[tardis.simulation.base][INFO   ]

        Current t_inner = 10703.212 K
        Expected t_inner for next iteration = 10673.712 K
 (base.py:652)
[tardis.simulation.base][INFO   ]

        Starting iteration 3 of 20 (base.py:444)
[tardis.simulation.base][INFO   ]

        Luminosity emitted   = 1.074e+43 erg / s
        Luminosity absorbed  = 3.391e+42 erg / s
        Luminosity requested = 1.059e+43 erg / s
 (base.py:657)
[tardis.simulation.base][INFO   ]
        Iteration converged 1/4 consecutive times. (base.py:260)
[tardis.simulation.base][INFO   ]

        Plasma stratification: (base.py:625)

Shell No.	t_rad	next_t_rad	w	next_w
0	1.08e+04 K	1.1e+04 K	0.525	0.483
5	1.1e+04 K	1.12e+04 K	0.203	0.189
10	1.08e+04 K	1.1e+04 K	0.125	0.118
15	1.05e+04 K	1.06e+04 K	0.0933	0.0895

[tardis.simulation.base][INFO   ]

        Current t_inner = 10673.712 K
        Expected t_inner for next iteration = 10635.953 K
 (base.py:652)
[tardis.simulation.base][INFO   ]

        Starting iteration 4 of 20 (base.py:444)
[tardis.simulation.base][INFO   ]

        Luminosity emitted   = 1.058e+43 erg / s
        Luminosity absorbed  = 3.352e+42 erg / s
        Luminosity requested = 1.059e+43 erg / s
 (base.py:657)
[tardis.simulation.base][INFO   ]
        Iteration converged 2/4 consecutive times. (base.py:260)
[tardis.simulation.base][INFO   ]

        Plasma stratification: (base.py:625)

Shell No.	t_rad	next_t_rad	w	next_w
0	1.1e+04 K	1.1e+04 K	0.483	0.469
5	1.12e+04 K	1.12e+04 K	0.189	0.182
10	1.1e+04 K	1.1e+04 K	0.118	0.113
15	1.06e+04 K	1.07e+04 K	0.0895	0.0861

[tardis.simulation.base][INFO   ]

        Current t_inner = 10635.953 K
        Expected t_inner for next iteration = 10638.407 K
 (base.py:652)
[tardis.simulation.base][INFO   ]

        Starting iteration 5 of 20 (base.py:444)
[tardis.simulation.base][INFO   ]

        Luminosity emitted   = 1.055e+43 erg / s
        Luminosity absorbed  = 3.399e+42 erg / s
        Luminosity requested = 1.059e+43 erg / s
 (base.py:657)
[tardis.simulation.base][INFO   ]
        Iteration converged 3/4 consecutive times. (base.py:260)
[tardis.simulation.base][INFO   ]

        Plasma stratification: (base.py:625)

Shell No.	t_rad	next_t_rad	w	next_w
0	1.1e+04 K	1.1e+04 K	0.469	0.479
5	1.12e+04 K	1.13e+04 K	0.182	0.178
10	1.1e+04 K	1.1e+04 K	0.113	0.113
15	1.07e+04 K	1.07e+04 K	0.0861	0.0839

[tardis.simulation.base][INFO   ]

        Current t_inner = 10638.407 K
        Expected t_inner for next iteration = 10650.202 K
 (base.py:652)
[tardis.simulation.base][INFO   ]

        Starting iteration 6 of 20 (base.py:444)
[tardis.simulation.base][INFO   ]

        Luminosity emitted   = 1.061e+43 erg / s
        Luminosity absorbed  = 3.398e+42 erg / s
        Luminosity requested = 1.059e+43 erg / s
 (base.py:657)
[tardis.simulation.base][INFO   ]
        Iteration converged 4/4 consecutive times. (base.py:260)
[tardis.simulation.base][INFO   ]

        Plasma stratification: (base.py:625)

Shell No.	t_rad	next_t_rad	w	next_w
0	1.1e+04 K	1.1e+04 K	0.479	0.47
5	1.13e+04 K	1.12e+04 K	0.178	0.185
10	1.1e+04 K	1.11e+04 K	0.113	0.112
15	1.07e+04 K	1.07e+04 K	0.0839	0.0856

[tardis.simulation.base][INFO   ]

        Current t_inner = 10650.202 K
        Expected t_inner for next iteration = 10645.955 K
 (base.py:652)
[tardis.simulation.base][INFO   ]

        Starting iteration 7 of 20 (base.py:444)
[tardis.simulation.base][INFO   ]

        Luminosity emitted   = 1.061e+43 erg / s
        Luminosity absorbed  = 3.382e+42 erg / s
        Luminosity requested = 1.059e+43 erg / s
 (base.py:657)
[tardis.simulation.base][INFO   ]
        Iteration converged 5/4 consecutive times. (base.py:260)
[tardis.simulation.base][INFO   ]

        Plasma stratification: (base.py:625)

Shell No.	t_rad	next_t_rad	w	next_w
0	1.1e+04 K	1.1e+04 K	0.47	0.47
5	1.12e+04 K	1.13e+04 K	0.185	0.178
10	1.11e+04 K	1.11e+04 K	0.112	0.112
15	1.07e+04 K	1.07e+04 K	0.0856	0.086

[tardis.simulation.base][INFO   ]

        Current t_inner = 10645.955 K
        Expected t_inner for next iteration = 10642.050 K
 (base.py:652)
[tardis.simulation.base][INFO   ]

        Starting iteration 8 of 20 (base.py:444)
[tardis.simulation.base][INFO   ]

        Luminosity emitted   = 1.062e+43 erg / s
        Luminosity absorbed  = 3.350e+42 erg / s
        Luminosity requested = 1.059e+43 erg / s
 (base.py:657)
[tardis.simulation.base][INFO   ]
        Iteration converged 6/4 consecutive times. (base.py:260)
[tardis.simulation.base][INFO   ]

        Plasma stratification: (base.py:625)

Shell No.	t_rad	next_t_rad	w	next_w
0	1.1e+04 K	1.11e+04 K	0.47	0.472
5	1.13e+04 K	1.14e+04 K	0.178	0.175
10	1.11e+04 K	1.11e+04 K	0.112	0.111
15	1.07e+04 K	1.07e+04 K	0.086	0.084

[tardis.simulation.base][INFO   ]

        Current t_inner = 10642.050 K
        Expected t_inner for next iteration = 10636.106 K
 (base.py:652)
[tardis.simulation.base][INFO   ]

        Starting iteration 9 of 20 (base.py:444)
[tardis.simulation.base][INFO   ]

        Luminosity emitted   = 1.052e+43 erg / s
        Luminosity absorbed  = 3.411e+42 erg / s
        Luminosity requested = 1.059e+43 erg / s
 (base.py:657)
[tardis.simulation.base][INFO   ]
        Iteration converged 7/4 consecutive times. (base.py:260)
[tardis.simulation.base][INFO   ]

        Plasma stratification: (base.py:625)

Shell No.	t_rad	next_t_rad	w	next_w
0	1.11e+04 K	1.11e+04 K	0.472	0.469
5	1.14e+04 K	1.15e+04 K	0.175	0.17
10	1.11e+04 K	1.11e+04 K	0.111	0.109
15	1.07e+04 K	1.08e+04 K	0.084	0.0822

[tardis.simulation.base][INFO   ]

        Current t_inner = 10636.106 K
        Expected t_inner for next iteration = 10654.313 K
 (base.py:652)
[tardis.simulation.base][INFO   ]

        Starting iteration 10 of 20 (base.py:444)
[tardis.simulation.base][INFO   ]

        Luminosity emitted   = 1.070e+43 erg / s
        Luminosity absorbed  = 3.335e+42 erg / s
        Luminosity requested = 1.059e+43 erg / s
 (base.py:657)
[tardis.simulation.base][INFO   ]
        Iteration converged 8/4 consecutive times. (base.py:260)
[tardis.simulation.base][INFO   ]

        Plasma stratification: (base.py:625)

Shell No.	t_rad	next_t_rad	w	next_w
0	1.11e+04 K	1.1e+04 K	0.469	0.475
5	1.15e+04 K	1.14e+04 K	0.17	0.177
10	1.11e+04 K	1.11e+04 K	0.109	0.112
15	1.08e+04 K	1.06e+04 K	0.0822	0.0878

[tardis.simulation.base][INFO   ]

        Current t_inner = 10654.313 K
        Expected t_inner for next iteration = 10628.190 K
 (base.py:652)
[tardis.simulation.base][INFO   ]

        Starting iteration 11 of 20 (base.py:444)
[tardis.simulation.base][INFO   ]

        Luminosity emitted   = 1.053e+43 erg / s
        Luminosity absorbed  = 3.363e+42 erg / s
        Luminosity requested = 1.059e+43 erg / s
 (base.py:657)
[tardis.simulation.base][INFO   ]
        Iteration converged 9/4 consecutive times. (base.py:260)
[tardis.simulation.base][INFO   ]

        Plasma stratification: (base.py:625)

Shell No.	t_rad	next_t_rad	w	next_w
0	1.1e+04 K	1.1e+04 K	0.475	0.472
5	1.14e+04 K	1.12e+04 K	0.177	0.184
10	1.11e+04 K	1.1e+04 K	0.112	0.114
15	1.06e+04 K	1.06e+04 K	0.0878	0.0859

[tardis.simulation.base][INFO   ]

        Current t_inner = 10628.190 K
        Expected t_inner for next iteration = 10644.054 K
 (base.py:652)
[tardis.simulation.base][INFO   ]

        Starting iteration 12 of 20 (base.py:444)
[tardis.simulation.base][INFO   ]

        Luminosity emitted   = 1.056e+43 erg / s
        Luminosity absorbed  = 3.420e+42 erg / s
        Luminosity requested = 1.059e+43 erg / s
 (base.py:657)
[tardis.simulation.base][INFO   ]
        Iteration converged 10/4 consecutive times. (base.py:260)
[tardis.simulation.base][INFO   ]

        Plasma stratification: (base.py:625)

Shell No.	t_rad	next_t_rad	w	next_w
0	1.1e+04 K	1.11e+04 K	0.472	0.467
5	1.12e+04 K	1.13e+04 K	0.184	0.176
10	1.1e+04 K	1.11e+04 K	0.114	0.11
15	1.06e+04 K	1.08e+04 K	0.0859	0.0821

[tardis.simulation.base][INFO   ]

        Current t_inner = 10644.054 K
        Expected t_inner for next iteration = 10653.543 K
 (base.py:652)
[tardis.simulation.base][INFO   ]

        Starting iteration 13 of 20 (base.py:444)
[tardis.simulation.base][INFO   ]

        Luminosity emitted   = 1.062e+43 erg / s
        Luminosity absorbed  = 3.406e+42 erg / s
        Luminosity requested = 1.059e+43 erg / s
 (base.py:657)
[tardis.simulation.base][INFO   ]
        Iteration converged 11/4 consecutive times. (base.py:260)
[tardis.simulation.base][INFO   ]

        Plasma stratification: (base.py:625)

Shell No.	t_rad	next_t_rad	w	next_w
0	1.11e+04 K	1.11e+04 K	0.467	0.466
5	1.13e+04 K	1.13e+04 K	0.176	0.18
10	1.11e+04 K	1.11e+04 K	0.11	0.111
15	1.08e+04 K	1.08e+04 K	0.0821	0.0841

[tardis.simulation.base][INFO   ]

        Current t_inner = 10653.543 K
        Expected t_inner for next iteration = 10647.277 K
 (base.py:652)
[tardis.simulation.base][INFO   ]

        Starting iteration 14 of 20 (base.py:444)
[tardis.simulation.base][INFO   ]

        Luminosity emitted   = 1.063e+43 erg / s
        Luminosity absorbed  = 3.369e+42 erg / s
        Luminosity requested = 1.059e+43 erg / s
 (base.py:657)
[tardis.simulation.base][INFO   ]
        Iteration converged 12/4 consecutive times. (base.py:260)
[tardis.simulation.base][INFO   ]

        Plasma stratification: (base.py:625)

Shell No.	t_rad	next_t_rad	w	next_w
0	1.11e+04 K	1.11e+04 K	0.466	0.469
5	1.13e+04 K	1.13e+04 K	0.18	0.182
10	1.11e+04 K	1.1e+04 K	0.111	0.113
15	1.08e+04 K	1.07e+04 K	0.0841	0.0854

[tardis.simulation.base][INFO   ]

        Current t_inner = 10647.277 K
        Expected t_inner for next iteration = 10638.875 K
 (base.py:652)
[tardis.simulation.base][INFO   ]

        Starting iteration 15 of 20 (base.py:444)
[tardis.simulation.base][INFO   ]

        Luminosity emitted   = 1.053e+43 erg / s
        Luminosity absorbed  = 3.417e+42 erg / s
        Luminosity requested = 1.059e+43 erg / s
 (base.py:657)
[tardis.simulation.base][INFO   ]
        Iteration converged 13/4 consecutive times. (base.py:260)
[tardis.simulation.base][INFO   ]

        Plasma stratification: (base.py:625)

Shell No.	t_rad	next_t_rad	w	next_w
0	1.11e+04 K	1.1e+04 K	0.469	0.484
5	1.13e+04 K	1.13e+04 K	0.182	0.181
10	1.1e+04 K	1.1e+04 K	0.113	0.113
15	1.07e+04 K	1.07e+04 K	0.0854	0.0858

[tardis.simulation.base][INFO   ]

        Current t_inner = 10638.875 K
        Expected t_inner for next iteration = 10655.125 K
 (base.py:652)
[tardis.simulation.base][INFO   ]

        Starting iteration 16 of 20 (base.py:444)
[tardis.simulation.base][INFO   ]

        Luminosity emitted   = 1.059e+43 erg / s
        Luminosity absorbed  = 3.445e+42 erg / s
        Luminosity requested = 1.059e+43 erg / s
 (base.py:657)
[tardis.simulation.base][INFO   ]
        Iteration converged 14/4 consecutive times. (base.py:260)
[tardis.simulation.base][INFO   ]

        Plasma stratification: (base.py:625)

Shell No.	t_rad	next_t_rad	w	next_w
0	1.1e+04 K	1.1e+04 K	0.484	0.472
5	1.13e+04 K	1.13e+04 K	0.181	0.177
10	1.1e+04 K	1.1e+04 K	0.113	0.113
15	1.07e+04 K	1.06e+04 K	0.0858	0.0858

[tardis.simulation.base][INFO   ]

        Current t_inner = 10655.125 K
        Expected t_inner for next iteration = 10655.561 K
 (base.py:652)
[tardis.simulation.base][INFO   ]

        Starting iteration 17 of 20 (base.py:444)
[tardis.simulation.base][INFO   ]

        Luminosity emitted   = 1.067e+43 erg / s
        Luminosity absorbed  = 3.372e+42 erg / s
        Luminosity requested = 1.059e+43 erg / s
 (base.py:657)
[tardis.simulation.base][INFO   ]
        Iteration converged 15/4 consecutive times. (base.py:260)
[tardis.simulation.base][INFO   ]

        Plasma stratification: (base.py:625)

Shell No.	t_rad	next_t_rad	w	next_w
0	1.1e+04 K	1.11e+04 K	0.472	0.468
5	1.13e+04 K	1.14e+04 K	0.177	0.175
10	1.1e+04 K	1.11e+04 K	0.113	0.11
15	1.06e+04 K	1.08e+04 K	0.0858	0.0816

[tardis.simulation.base][INFO   ]

        Current t_inner = 10655.561 K
        Expected t_inner for next iteration = 10636.536 K
 (base.py:652)
[tardis.simulation.base][INFO   ]

        Starting iteration 18 of 20 (base.py:444)
[tardis.simulation.base][INFO   ]

        Luminosity emitted   = 1.057e+43 erg / s
        Luminosity absorbed  = 3.365e+42 erg / s
        Luminosity requested = 1.059e+43 erg / s
 (base.py:657)
[tardis.simulation.base][INFO   ]
        Iteration converged 16/4 consecutive times. (base.py:260)
[tardis.simulation.base][INFO   ]

        Plasma stratification: (base.py:625)

Shell No.	t_rad	next_t_rad	w	next_w
0	1.11e+04 K	1.11e+04 K	0.468	0.464
5	1.14e+04 K	1.13e+04 K	0.175	0.177
10	1.11e+04 K	1.1e+04 K	0.11	0.113
15	1.08e+04 K	1.07e+04 K	0.0816	0.0848

[tardis.simulation.base][INFO   ]

        Current t_inner = 10636.536 K
        Expected t_inner for next iteration = 10641.692 K
 (base.py:652)
[tardis.simulation.base][INFO   ]

        Starting iteration 19 of 20 (base.py:444)
[tardis.simulation.base][INFO   ]

        Luminosity emitted   = 1.056e+43 erg / s
        Luminosity absorbed  = 3.405e+42 erg / s
        Luminosity requested = 1.059e+43 erg / s
 (base.py:657)
[tardis.simulation.base][INFO   ]
        Iteration converged 17/4 consecutive times. (base.py:260)
[tardis.simulation.base][INFO   ]

        Plasma stratification: (base.py:625)

Shell No.	t_rad	next_t_rad	w	next_w
0	1.11e+04 K	1.11e+04 K	0.464	0.466
5	1.13e+04 K	1.13e+04 K	0.177	0.177
10	1.1e+04 K	1.11e+04 K	0.113	0.111
15	1.07e+04 K	1.07e+04 K	0.0848	0.0853

[tardis.simulation.base][INFO   ]

        Current t_inner = 10641.692 K
        Expected t_inner for next iteration = 10650.463 K
 (base.py:652)
[tardis.simulation.base][INFO   ]

        Simulation finished in 19 iterations
        Simulation took 23.11 s
 (base.py:542)
[tardis.simulation.base][INFO   ]

        Starting iteration 20 of 20 (base.py:444)
[tardis.simulation.base][INFO   ]

        Luminosity emitted   = 1.061e+43 erg / s
        Luminosity absorbed  = 3.401e+42 erg / s
        Luminosity requested = 1.059e+43 erg / s
 (base.py:657)

HDF

TARDIS can save simulation data to HDF files for later analysis. The code below shows how to load a simulation from an HDF file. This is useful when you want to analyze simulation results without re-running the simulation.

# import astropy.units as u
# import pandas as pd

# hdf_fpath = "add_file_path_here"
# with pd.HDFStore(hdf_fpath, "r") as hdf:
#     sim = u.Quantity(hdf["/simulation"])

Calculating Spectral Quantities

This section demonstrates how to analyze the spectral decomposition (SDEC) output from TARDIS. SDEC helps understand how different atomic species contribute to the formation of spectral features through their emission and absorption processes.

Wavelength and Frequency Grid Setup

First, we set up the wavelength and frequency grids for analysis and plotting

# Get frequency bins from the simulation's real packets spectrum
plot_frequency_bins = (
    simulation.spectrum_solver.spectrum_real_packets._frequency
)

# Get wavelength array from the simulation's real packets spectrum
plot_wavelength = simulation.spectrum_solver.spectrum_real_packets.wavelength

# Create frequency array by excluding the last bin edge
plot_frequency = plot_frequency_bins[:-1]

# Create a boolean mask array of ones with same size as wavelength array
# This will be used to filter wavelength ranges later
packet_wvl_range_mask = np.ones(plot_wavelength.size, dtype=bool)

Luminosity Calculations

Then we calculate emission and absorption luminosities for each atomic species to understand their contributions to the spectrum

# Calculate emission luminosities for each species
# Returns a DataFrame with species-wise emission contributions and list of species
(emission_luminosities_df, emission_species) = calculate_emission_luminosities(
    simulation,
    "real",  # Use real packets mode
    packet_wvl_range=None,  # Use full wavelength range
)

# Calculate absorption luminosities for each species
# Returns a DataFrame with species-wise absorption contributions and list of species
(
    absorption_luminosities_df,
    absorption_species,
) = calculate_absorption_luminosities(
    simulation,
    packets_mode="real",  # Use real packets mode
    packet_wvl_range=None,  # Use full wavelength range
)

# Calculate total luminosities by adding absorption and emission
# Drop 'no interaction' and 'electron scattering' columns from emission before adding
total_luminosities_df = (
    absorption_luminosities_df
    + emission_luminosities_df.drop(["noint", "escatter"], axis=1)
)

# Convert atomic numbers to element symbols for plotting
species = np.array(list(total_luminosities_df.keys()))
species_name = [
    atomic_number2element_symbol(atomic_num)
    for atomic_num in species
]
species_length = len(species_name)

# Get the modeled spectrum luminosity for the selected wavelength range
modeled_spectrum_luminosity = (
    simulation.spectrum_solver.spectrum_real_packets.luminosity_density_lambda[
        packet_wvl_range_mask
    ]
)

Plotting Species-wise Emission and Absorption Contributions

Now we will plot the species-wise emission and absorption contributions that we calculated in the previous step using calculate_emission_luminosities() and calculate_absorption_luminosities() respectively. We’ll start with packets that had no interaction with the ejecta, followed by electron scattering events, and then contributions from different atomic species. The colorbar will help you to identify the atomic species and their respective contributions to the total luminosity at each wavelength.

Matplotlib

We will now plot the emission contributions from different species.

# Import required matplotlib modules for visualization
import matplotlib.cm as cm  # For colormaps
import matplotlib.colors as clr  # For color mapping and normalization
import matplotlib.pyplot as plt  # Main plotting interface

# Create figure and axis
ax = plt.figure(figsize=(12, 7)).add_subplot(111)

# Generate colors for each species using 'jet' colormap
# The colors will be evenly distributed across the colormap based on number of species
cmap = plt.get_cmap("jet", species_length)
color_list = []
color_values = []
for species_counter in range(species_length):
    # Calculate normalized color value and store in both lists
    color = cmap(species_counter / species_length)
    color_list.append(color)
    color_values.append(color)


# Create a custom colormap for species
# This will be used in the colorbar to show which color corresponds to which species
custcmap = clr.ListedColormap(color_values) # Normalize color range
norm = clr.Normalize(vmin=0, vmax=species_length)
mappable = cm.ScalarMappable(norm=norm, cmap=custcmap)
mappable.set_array(np.linspace(1, species_length + 1, 256))


# Add colorbar for species representation
cbar = plt.colorbar(mappable, ax=ax)
bounds = np.arange(species_length) + 0.5
cbar.set_ticks(bounds)
cbar.set_ticklabels(species_name) # Label ticks with species names

# Plot the emission contributions starting with 'no interaction' packets
lower_level = np.zeros(emission_luminosities_df.shape[0])  # Start at zero
upper_level = lower_level + emission_luminosities_df.noint.to_numpy()
ax.fill_between(
    plot_wavelength.value,
    lower_level,
    upper_level,
    color="#4C4C4C", #  Dark grey color
    label="No interaction",
)

# Plot electron scattering contribution
# This is stacked on top of the 'no interaction' contribution
lower_level = upper_level # Start where 'no interaction' ended
upper_level = lower_level + emission_luminosities_df.escatter.to_numpy()
ax.fill_between(
    plot_wavelength.value,
    lower_level,
    upper_level,
    color="#8F8F8F",
    label="Electron Scatter Only",
)

# Plot the emission contributions for each species by stacking them on top of each other.
# The lower level for each species starts where the previous species ended (upper_level),
# and the upper level adds that species' emission contribution.
# Each species gets a unique color from the colormap defined above.
for species_counter, identifier in enumerate(species):
    lower_level = upper_level
    upper_level = lower_level + emission_luminosities_df[identifier].to_numpy()

    ax.fill_between(
        plot_wavelength.value,
        lower_level,
        upper_level,
        color=color_list[species_counter],
        cmap=cmap,
        linewidth=0,
    )


# The absorption contributions are plotted below zero, showing how each species
# absorbs light and reduces the total luminosity. The contributions are stacked
# downward from zero, with each species getting the same color as its emission above.
lower_level = np.zeros(absorption_luminosities_df.shape[0])
for species_counter, identifier in enumerate(species):
    upper_level = lower_level
    lower_level = (
        upper_level - absorption_luminosities_df[identifier].to_numpy()
    )

    ax.fill_between(
        plot_wavelength.value,
        upper_level,
        lower_level,
        color=color_list[species_counter],
        cmap=cmap,
        linewidth=0,
    )


# Plot the modeled spectrum
ax.plot(
    plot_wavelength.value,
    modeled_spectrum_luminosity.value,
    "--b",
    label="Real Spectrum",
    linewidth=1,
)

# Add labels, legend, and formatting
ax.legend(fontsize=12)
ax.set_xlabel(r"Wavelength $[\mathrm{\AA}]$", fontsize=12)
ax.set_ylabel(
    r"$L_{\lambda}$ [erg $\mathrm{s^{-1}}$ $\mathrm{\AA^{-1}}$]",
    fontsize=12,
)
plt.gca()

<Axes: xlabel='Wavelength $[\\mathrm{\\AA}]$', ylabel='$L_{\\lambda}$ [erg $\\mathrm{s^{-1}}$ $\\mathrm{\\AA^{-1}}$]'>

Plotly

import plotly.graph_objects as go

# Create figure
fig = go.Figure()

# By specifying a common stackgroup, plotly will itself add up luminosities,
# in order, to created stacked area chart
fig.add_trace(
    go.Scatter(
        x=emission_luminosities_df.index,
        y=emission_luminosities_df.noint,
        mode="none",
        name="No interaction",
        fillcolor="#4C4C4C",
        stackgroup="emission",
        hovertemplate="(%{x:.2f}, %{y:.3g})",
    )
)

fig.add_trace(
    go.Scatter(
        x=emission_luminosities_df.index,
        y=emission_luminosities_df.escatter,
        mode="none",
        name="Electron Scatter Only",
        fillcolor="#8F8F8F",
        stackgroup="emission",
        hoverlabel={"namelength": -1},
        hovertemplate="(%{x:.2f}, %{y:.3g})",
    )
)

# The species data comes from the emission_luminosities_df and absorption_luminosities_df DataFrames
# We plot emissions as positive values stacked above zero
# And absorptions as negative values stacked below zero
# This creates a visualization showing how different species contribute to the spectrum
for (species_counter, identifier), name_of_spec in zip(
    enumerate(species), species_name
):
    fig.add_trace(
        go.Scatter(
            x=emission_luminosities_df.index,
            y=emission_luminosities_df[identifier],
            mode="none",
            name=name_of_spec + " Emission",
            hovertemplate=f"<b>{name_of_spec:s} Emission<br>"  # noqa: ISC003
            + "(%{x:.2f}, %{y:.3g})<extra></extra>",
            fillcolor=pu.to_rgb255_string(color_list[species_counter]),
            stackgroup="emission",
            showlegend=False,
            hoverlabel={"namelength": -1},
        )
    )
    # Plot absorption part
    fig.add_trace(
        go.Scatter(
            x=absorption_luminosities_df.index,
            # to plot absorption luminosities along negative y-axis
            y=absorption_luminosities_df[identifier] * -1,
            mode="none",
            name=name_of_spec + " Absorption",
            hovertemplate=f"<b>{name_of_spec:s} Absorption<br>"  # noqa: ISC003
            + "(%{x:.2f}, %{y:.3g})<extra></extra>",
            fillcolor=pu.to_rgb255_string(color_list[species_counter]),
            stackgroup="absorption",
            showlegend=False,
            hoverlabel={"namelength": -1},
        )
    )

# Plot modeled spectrum
fig.add_trace(
    go.Scatter(
        x=plot_wavelength.value,
        y=modeled_spectrum_luminosity.value,
        mode="lines",
        line={
            "color": "blue",
            "width": 1,
        },
        name="Real Spectrum",
        hovertemplate="(%{x:.2f}, %{y:.3g})",
        hoverlabel={"namelength": -1},
    )
)

# Interpolate [0, 1] range to create bins equal to number of elements
colorscale_bins = np.linspace(0, 1, num=len(species_name) + 1)

# Create a categorical colorscale [a list of (reference point, color)]
# by mapping same reference points (excluding 1st and last bin edge)
# twice in a row (https://plotly.com/python/colorscales/#constructing-a-discrete-or-discontinuous-color-scale)
categorical_colorscale = []
for species_counter in range(len(species_name)):
    color = pu.to_rgb255_string(cmap(colorscale_bins[species_counter]))
    categorical_colorscale.append((colorscale_bins[species_counter], color))
    categorical_colorscale.append((colorscale_bins[species_counter + 1], color))

# Create a categorical colorscale for the elements by mapping each species to a color
coloraxis_options = {
    "colorscale": categorical_colorscale,
    "showscale": True,
    "cmin": 0,
    "cmax": len(species_name),
    "colorbar": {
        "title": "Elements",
        "tickvals": np.arange(0, len(species_name)) + 0.5,
        "ticktext": species_name,
        # to change length and position of colorbar
        "len": 0.75,
        "yanchor": "top",
        "y": 0.75,
    },
}

# Add an invisible scatter point to make the colorbar show up in the plot
# The point is placed at the middle of the wavelength range with y=0
# The marker color is set to 0 (first color in colorscale) with opacity=0 to hide it
# coloraxis_options contains the categorical colorscale mapping species to colors
scatter_point_idx = pu.get_mid_point_idx(plot_wavelength)
fig.add_trace(
    go.Scatter(
        x=[plot_wavelength[scatter_point_idx].value],
        y=[0],
        mode="markers",
        name="Colorbar",
        showlegend=False,
        hoverinfo="skip",
        marker=dict(color=[0], opacity=0, **coloraxis_options),
    )
)

# Set label and other layout options
xlabel = pu.axis_label_in_latex("Wavelength", u.AA)
ylabel = pu.axis_label_in_latex(
    "L_{\\lambda}", u.Unit("erg/(s AA)"), only_text=False
)
fig.update_layout(
    xaxis={
        "title": xlabel,
        "exponentformat": "none",
    },
    yaxis={"title": ylabel, "exponentformat": "e"},
    height=600,
)