diff --git a/tardis/energy_input/energy_source.py b/tardis/energy_input/energy_source.py
index 1be7f6e5a74..4fa33734944 100644
--- a/tardis/energy_input/energy_source.py
+++ b/tardis/energy_input/energy_source.py
@@ -190,7 +190,7 @@ def get_nuclear_lines_database(
     pandas DataFrame
         The decay radiation lines
     """
-    decay_radiation_db = pd.read_hdf(path, "decay_data")
+    decay_radiation_db = pd.read_hdf(path, "decay_radiation_data")
     return decay_radiation_db
 
 
diff --git a/tardis/energy_input/gamma_ray_transport.py b/tardis/energy_input/gamma_ray_transport.py
index 6aea2856080..7542572ce2b 100644
--- a/tardis/energy_input/gamma_ray_transport.py
+++ b/tardis/energy_input/gamma_ray_transport.py
@@ -178,7 +178,6 @@ def initialize_packets(
 
     packet_index = 0
     for k, shell in enumerate(decays_per_shell):
-
         initial_radii = initial_packet_radius(
             shell, inner_velocities[k], outer_velocities[k]
         )
@@ -261,7 +260,6 @@ def initialize_packets(
 
 
 def calculate_shell_masses(model):
-
     """Function to calculate shell masses
     Parameters
     ----------
@@ -274,186 +272,162 @@ def calculate_shell_masses(model):
 
     """
 
-    ejecta_density = model.density.to("g/cm^3").value
-    ejecta_volume = model.volume.to("cm^3").value
-    shell_masses = ejecta_volume * ejecta_density
-
-    return shell_masses
-
-
-def main_gamma_ray_loop(
-    num_decays,
-    model,
-    plasma,
-    time_steps=10,
-    time_end=80.0,
-    grey_opacity=-1,
-    spectrum_bins=500,
-    time_space="log",
-    photoabsorption_opacity="tardis",
-    pair_creation_opacity="tardis",
-    seed=1,
-    path_to_decay_data="~/Downloads/tardisnuclear/decay_radiation.h5",
-    positronium_fraction=0.0,
-):
-    """Main loop that determines the gamma ray propagation
+    ejecta_density = model.density.to("g/cm^3")
+    ejecta_volume = model.volume.to("cm^3")
+    return (ejecta_volume * ejecta_density).to(u.g)
 
+
+def calculate_total_decays(inventories, time_delta):
+    """Function to create inventories of isotope
     Parameters
     ----------
-    num_decays : int
-        Number of decays requested
-    model : tardis.SimulationState
+    model : tardis.Radial1DModel
         The tardis model to calculate gamma ray propagation through
-    plasma : tardis.plasma.BasePlasma
-        The tardis plasma with calculated atomic number density
-    time_steps : int
-        Number of time steps requested
+
     time_end : float
         End time of simulation in days
-    grey_opacity : float
-        Grey photoabsorption opacity for gamma-rays in cm^2 g^-1, set to -1 to turn off
-    spectrum_bins : int
-        Number of energy bins for the gamma-ray spectrum
-    time_space : str
-        `'log'` for log-space time steps, otherwise time steps will be linear
-    photoabsorption_opacity : str
-        Set the photoabsorption opacity calculation.
-        Defaults to Ambwani & Sutherland (1988) approximation.
-        `'kasen'` uses the Kasen et al. 2006 method.
-    pair_creation_opacity : str
-        Set the pair creation opacity calculation. Defaults to Ambwani & Sutherland (1988) approximation.
-        `'artis'` uses the ARTIS implementation of the Ambwani & Sutherland (1988) approximation.
-    seed : int
-        Sets the seed for the random number generator. Uses deprecated methods.
-    path_to_decay_data : str
-        The path to a decay radiation file from the `nuclear` package.
+    Returns
+    -------
+        Total decay list : List
+            list of total decays for x g of isotope for time 't'
 
+    """
+
+    time_delta = u.Quantity(time_delta, u.s)
+
+    total_decays_list = []
+    for inv in inventories:
+        total_decays = inv.cumulative_decays(time_delta.value)
+        total_decays_list.append(total_decays)
+
+    return total_decays_list
+
+
+def create_isotope_dicts(raw_isotope_abundance, shell_masses):
+    """
+    Function to create a dictionary of isotopes for each shell with their masses.
+    Parameters
+    ----------
+    raw_isotope_abundance : pd.DataFrame
+        isotope abundance in mass fractions.
+    shell_masses : numpy.ndarray
+        shell masses in units of g
     Returns
     -------
-    pandas.DataFrame
-        Energy per shell per time in units of eV/s/cm^-3
-    pandas.DataFrame
-        Columns:
-        packet index,
-        Energy input per packet,
-        radius of deposition,
-        time of deposition,
-        compton opacity,
-        photoabsorption opacity,
-        pair creation opacity
-    pandas.DataFrame
-        Energy of escaping packets
-    numpy.ndarray
-        Packets emitted per shell
-    pandas.DataFrame
-        Energy from positrons
-    pandas.DataFrame
-        Estimated energy deposition in units of keV/s/cm^-3
+        isotope_dicts : Dict
+            dictionary of isotopes for each shell with their ``masses``.
+            For eg: {0: {(28, 56): {'Ni56': 0.0001}, (27, 57): {'Co56': 0.0001}}
+                    {1: {(28, 56): {'Ni56': 0.0001}, (27, 57): {'Co56': 0.0001}}} etc
+
     """
-    # Note: not best numpy practice, but works better in numba than the alternatives
-    np.random.seed(seed)
-
-    # Enforce cgs
-    outer_velocities = model.v_outer.to("cm/s").value
-    inner_velocities = model.v_inner.to("cm/s").value
-    ejecta_density = model.density.to("g/cm^3").value
-    ejecta_volume = model.volume.to("cm^3").value
-    ejecta_velocity_volume = (
-        4 * np.pi / 3 * (outer_velocities**3.0 - inner_velocities**3.0)
-    )
-    time_explosion = model.time_explosion.to("s").value
-    number_of_shells = model.no_of_shells
-    raw_isotope_abundance = model.raw_isotope_abundance.sort_values(
-        by=["atomic_number", "mass_number"], ascending=False
-    )
+    isotope_dicts = {}
+    for i in range(len(raw_isotope_abundance.columns)):
+        isotope_dicts[i] = {}
+        for (
+            atomic_number,
+            mass_number,
+        ), abundances in raw_isotope_abundance.iterrows():
+            isotope_dicts[i][atomic_number, mass_number] = {}
+            nuclear_symbol = f"{rd.utils.Z_to_elem(atomic_number)}{mass_number}"
+            isotope_dicts[i][atomic_number, mass_number][nuclear_symbol] = (
+                abundances[i] * shell_masses[i].to(u.g).value
+            )
 
-    shell_masses = calculate_shell_masses(model)
+    return isotope_dicts
 
-    time_start = time_explosion
-    time_end *= u.d.to(u.s)
 
-    assert (
-        time_start < time_end
-    ), "Error, simulation start time greater than end time!"
+def create_inventories_dict(isotope_dict):
+    """Function to create dictionary of inventories for each shell
+    Parameters
+    ----------
+    isotope_dict : Dict
+        dictionary of isotopes for each shell with their ``masses``.
+    Returns
+    -------
+        inv : Dict
+            dictionary of inventories for each shell
+            For eg: {0: {'Ni56': <radioactivedecay.Inventory at 0x7f7d2c0d0e50>,
+                         'Co56': <radioactivedecay.Inventory at 0x7f7d2c0d0e50>},
+                    {1: {'Ni56': <radioactivedecay.Inventory at 0x7f7d2c0d0e50>,
+                         'Co56': <radioactivedecay.Inventory at 0x7f7d2c0d0e50>}} etc
+    """
+    inv = {}
+    for shell, isotopes in isotope_dict.items():
+        inv[shell] = {}
+        for isotope, name in isotopes.items():
+            inv[shell][isotope] = rd.Inventory(name, "g")
 
-    if time_space == "log":
-        times = np.geomspace(time_start, time_end, time_steps + 1)
-    else:
-        times = np.linspace(time_start, time_end, time_steps + 1)
+    return inv
 
-    dt_array = np.diff(times)
-    effective_time_array = np.array(
-        [np.sqrt(times[i] * times[i + 1]) for i in range(time_steps)]
-    )
 
-    # Use isotopic number density
-    for atom_number in plasma.isotope_number_density.index.get_level_values(0):
-        values = plasma.isotope_number_density.loc[atom_number].values
-        if values.shape[1] > 1:
-            plasma.number_density.loc[atom_number] = np.sum(values, axis=0)
-        else:
-            plasma.number_density.loc[atom_number] = values
-
-    # Calculate electron number density
-    electron_number_density = (
-        plasma.number_density.mul(plasma.number_density.index, axis=0)
-    ).sum()
-
-    electron_number_density_time = np.zeros(
-        (len(ejecta_velocity_volume), len(effective_time_array))
-    )
+def calculate_total_decays(inventory_dict, time_delta):
+    """
+    Function to calculate total decays for each isotope in each shell
+    Parameters
+    ----------
+    inventory_dict : Dict
+        dictionary of inventories for each shell
+    time_delta : float
+        time interval in units of time (days/mins/secs etc)
+    Returns
+    -------
+        total_decays : Dict
+            dictionary of total decays for each isotope in each shell
 
-    mass_density_time = np.zeros(
-        (len(ejecta_velocity_volume), len(effective_time_array))
-    )
+    """
+    time_delta = u.Quantity(time_delta, u.s)
+    total_decays = {}
+    for shell, isotopes in inventory_dict.items():
+        total_decays[shell] = {}
+        for isotope, name in isotopes.items():
+            total_decays[shell][isotope] = name.cumulative_decays(
+                time_delta.value
+            )
 
-    electron_number = (electron_number_density * ejecta_volume).to_numpy()
+    return total_decays
 
-    inv_volume_time = np.zeros(
-        (len(ejecta_velocity_volume), len(effective_time_array))
-    )
 
-    # Pre-calculate quantities as they change with time
-    for i, t in enumerate(effective_time_array):
-        inv_volume_time[:, i] = (1.0 / ejecta_velocity_volume) / (t**3.0)
-        mass_density_time[:, i] = shell_masses * inv_volume_time[:, i]
-        electron_number_density_time[:, i] = (
-            electron_number * inv_volume_time[:, i]
-        )
+def calculate_average_energies(raw_isotope_abundance, gamma_ray_lines):
+    """
+    Function to calculate average energies of positrons and gamma rays
+    from a list of gamma ray lines from nndc.
+    Parameters
+    ----------
+    raw_isotope_abundance : pd.DataFrame
+        isotope abundance in mass fractions
+    gamma_ray_lines : pd.DataFrame
+        decay data
 
-    energy_df_rows = np.zeros((number_of_shells, time_steps))
+    Returns
+    -------
+    average_energies_list : List
+        list of gamma ray energies
+    average_positron_energies_list : List
+        list of positron energies
+    gamma_ray_line_array_list : List
+        list of gamma ray lines
 
-    # Calculate number of packets per shell based on the mass of isotopes
-    number_of_isotopes = plasma.isotope_number_density * ejecta_volume
-    total_number_isotopes = number_of_isotopes.sum(axis=1)
+    """
 
-    inventories = raw_isotope_abundance.to_inventories()
     all_isotope_names = get_all_isotopes(raw_isotope_abundance)
     all_isotope_names.sort()
 
-    gamma_ray_lines = get_nuclear_lines_database(path_to_decay_data)
-
-    taus = {}
-    parents = {}
     gamma_ray_line_array_list = []
     average_energies_list = []
     average_positron_energies_list = []
 
-    for i, isotope in enumerate(all_isotope_names):
-        nuclide = rd.Nuclide(isotope)
-        taus[isotope] = nuclide.half_life() / np.log(2)
-        child = nuclide.progeny()
-        if child is not None:
-            for c in child:
-                if rd.Nuclide(c).half_life("readable") != "stable":
-                    parents[c] = isotope
+    gamma_ray_line_dict = {}
+    average_energies = {}
+    average_positron_energies = {}
 
+    for i, isotope in enumerate(all_isotope_names):
         energy, intensity = setup_input_energy(
             gamma_ray_lines[gamma_ray_lines.index == isotope.replace("-", "")],
             "g",
         )
+        average_energies_list.append(np.sum(energy * intensity))  # keV
         gamma_ray_line_array_list.append(np.stack([energy, intensity]))
-        average_energies_list.append(np.sum(energy * intensity))
+
         positron_energy, positron_intensity = setup_input_energy(
             gamma_ray_lines[gamma_ray_lines.index == isotope.replace("-", "")],
             "bp",
@@ -462,13 +436,8 @@ def main_gamma_ray_loop(
             np.sum(positron_energy * positron_intensity)
         )
 
-    # Construct Numba typed dicts
-    gamma_ray_line_arrays = {}
-    average_energies = {}
-    average_positron_energies = {}
-
     for iso, lines in zip(all_isotope_names, gamma_ray_line_array_list):
-        gamma_ray_line_arrays[iso] = lines
+        gamma_ray_line_dict[iso] = lines
 
     for iso, energy, positron_energy in zip(
         all_isotope_names, average_energies_list, average_positron_energies_list
@@ -476,222 +445,129 @@ def main_gamma_ray_loop(
         average_energies[iso] = energy
         average_positron_energies[iso] = positron_energy
 
-    # urilight chooses to have 0 as the baseline for this calculation
-    # but time_start may also be valid in which case decay time is time_end - time_start
-    total_energy_list = []
-
-    for shell, inv in enumerate(inventories):
-        decayed_energy = {}
-        total_decays = inv.cumulative_decays(time_end)
-        for nuclide in total_decays:
-            if nuclide in parents and nuclide != "Co-56" and nuclide != "Co-57":
-                parent = parents[nuclide]
-                if parent in parents:
-                    parent = parents[parent]
-                decayed_energy[parent] += (
-                    total_decays[nuclide]
-                    * average_energies[nuclide]
-                    * shell_masses[shell]
-                )
-            else:
-                decayed_energy[nuclide] = (
-                    total_decays[nuclide]
-                    * average_energies[nuclide]
-                    * shell_masses[shell]
-                )
-
-        total_energy_list.append(decayed_energy)
-
-    total_energy = pd.DataFrame(total_energy_list)
-
-    total_energy_columns = total_energy.columns.to_list()
-
-    total_energy = total_energy[
-        sorted(
-            total_energy_columns, key=get_nuclide_atomic_number, reverse=True
-        )
-    ]
-
-    energy_per_mass = total_energy.divide(
-        (raw_isotope_abundance * shell_masses).T.to_numpy(),
-        axis=0,
-    )
-
-    # Time averaged energy per mass for constant packet count
-    average_power_per_mass = energy_per_mass / (time_end - time_start)
-
-    energy_per_mass_norm = energy_per_mass.divide(
-        energy_per_mass.sum(axis=1), axis=0
-    )  # .cumsum(axis=1)
-
-    decayed_packet_count = num_decays * number_of_isotopes.divide(
-        total_number_isotopes, axis=0
-    )
-
-    packets_per_isotope = (
-        (energy_per_mass_norm * decayed_packet_count.T.values)
-        .round()
-        .fillna(0)
-        .astype(int)
+    return (
+        average_energies,
+        average_positron_energies,
+        gamma_ray_line_dict,
     )
 
-    print("Total gamma-ray energy")
-    print(total_energy.sum().sum() * u.keV.to("erg"))
 
-    print("Total positron energy")
-    print(total_energy["Co-56"].sum(axis=0) * 0.0337 * u.keV.to("erg"))
-
-    # Taking iron group to be elements 21-30
-    # Used as part of the approximations for photoabsorption and pair creation
-    # Dependent on atomic data
-    iron_group_fraction_per_shell = model.abundance.loc[(21):(30)].sum(axis=0)
+def get_taus(raw_isotope_abundance):
+    """
+    Function to calculate taus for each isotope
+    Parameters
+    ----------
+    raw_isotope_abundance : pd.DataFrame
+        isotope abundance in mass fractions
 
-    number_of_packets = packets_per_isotope.sum().sum()
-    print("Total packets:", number_of_packets)
+    Returns
+    -------
+    taus : Dict
+        dictionary of taus for each isotope
+    parents : Dict
+        dictionary of parents for each isotope
+    """
+    all_isotope_names = get_all_isotopes(raw_isotope_abundance)
+    all_isotope_names.sort()
 
-    packet_energy = total_energy.sum().sum() / number_of_packets
+    taus = {}
+    parents = {}
+    for isotope in all_isotope_names:
+        nuclide = rd.Nuclide(isotope)
+        taus[isotope] = nuclide.half_life() / np.log(2)
+        child = nuclide.progeny()
+        if child is not None:
+            for c in child:
+                if rd.Nuclide(c).half_life("readable") != "stable":
+                    parents[isotope] = c
 
-    print("Energy per packet", packet_energy)
+    return taus, parents
 
-    # Need to update volume for positron deposition to be time-dependent
-    print("Initializing packets")
-    (
-        packets,
-        energy_df_rows,
-        energy_plot_df_rows,
-        energy_plot_positron_rows,
-    ) = initialize_packets(
-        packets_per_isotope,
-        packet_energy,
-        gamma_ray_line_arrays,
-        positronium_fraction,
-        inner_velocities,
-        outer_velocities,
-        inv_volume_time,
-        times,
-        energy_df_rows,
-        effective_time_array,
-        taus,
-        parents,
-        average_positron_energies,
-        inventories,
-        average_power_per_mass,
-    )
 
-    print("Total positron energy from packets")
-    print((energy_df_rows).sum().sum() * u.eV.to("erg"))
+def decay_chain_energies(
+    average_energies,
+    total_decays,
+):
+    """
+    Function to calculate decay chain energies.
+    Parameters
+    ----------
+    raw_isotope_abundance : pd.DataFrame
+        isotope abundance in mass fractions
+    average_energies_list : List
+        list of gamma ray energies
+    average_positron_energies_list : List
+        list of positron energies
+    gamma_ray_line_array_list : List
+        list of gamma ray lines
+    total_decays : Dict
+        dictionary of total decays for each isotope in each shell
+    Returns
+    -------
+    decay_energy : Dict
+        dictionary of decay chain energies for each isotope in each shell
 
-    total_cmf_energy = 0
-    total_rf_energy = 0
+    """
+    decay_energy = {}
+    for shell, isotopes in total_decays.items():
+        decay_energy[shell] = {}
+        for name, isotope in isotopes.items():
+            decay_energy[shell][name] = {}
+            for iso, dps in isotope.items():
+                decay_energy[shell][name][iso] = dps * average_energies[iso]
 
-    for p in packets:
-        total_cmf_energy += p.energy_cmf
-        total_rf_energy += p.energy_rf
+    return decay_energy
 
-    print("Total CMF energy")
-    print(total_cmf_energy)
 
-    # Below is the Artis compensation for their method of packet rejection
+def calculate_energy_per_mass(
+    decay_energy, raw_isotope_abundance, shell_masses
+):
     """
-    energy_ratio = total_energy.sum().sum() / total_cmf_energy
-
-    print("Energy ratio")
-    print(energy_ratio)
-    
-    for p in packets:
-        p.energy_cmf *= energy_ratio
-        p.energy_rf *= energy_ratio
-
-    for e in energy_df_rows:
-        e *= energy_ratio
-    
-    for row in energy_plot_df_rows:
-        row[1] *= energy_ratio
+    Function to calculate decay energy per mass for each isotope chain.
+    Parameters
+    ----------
+    decay_energy : Dict
+        dictionary of decay chain energies for each isotope in each shell
+    raw_isotope_abundance : pd.DataFrame
+        isotope abundance in mass fractions.
+    shell_masses : numpy.ndarray
+        shell masses in units of g
+    Returns
+    -------
+    energy_per_mass : pd.DataFrame
+        decay energy per mass for each isotope chain
+        For e.g Ni56 has 2 decay chains:
+        Ni56 -> Co56 -> Fe56. It will calculate the decay energy per mass for each chain
+        and store as a dataframe.
     """
-    print("Total RF energy")
-    print(total_rf_energy)
-
-    energy_bins = np.logspace(2, 3.8, spectrum_bins)
-    energy_out = np.zeros((len(energy_bins - 1), time_steps))
-
-    # Process packets
-    (
-        energy_df_rows,
-        energy_plot_df_rows,
-        energy_out,
-        deposition_estimator,
-    ) = gamma_packet_loop(
-        packets,
-        grey_opacity,
-        photoabsorption_opacity,
-        pair_creation_opacity,
-        electron_number_density_time,
-        mass_density_time,
-        inv_volume_time,
-        iron_group_fraction_per_shell.to_numpy(),
-        inner_velocities,
-        outer_velocities,
-        times,
-        dt_array,
-        effective_time_array,
-        energy_bins,
-        energy_df_rows,
-        energy_plot_df_rows,
-        energy_out,
-    )
-
-    # DataFrame of energy information
-    energy_plot_df = pd.DataFrame(
-        data=energy_plot_df_rows,
-        columns=[
-            "packet_index",
-            "energy_input",
-            "energy_input_r",
-            "energy_input_time",
-            "energy_input_type",
-            "compton_opacity",
-            "photoabsorption_opacity",
-            "total_opacity",
-        ],
-    )
+    energy_dict = {}
+    for shell, isotopes in decay_energy.items():
+        energy_dict[shell] = {}
+        for name, isotope in isotopes.items():
+            energy_dict[shell][name] = sum(isotope.values())
+
+    energy_list = []
+    for shell, isotopes in energy_dict.items():
+        for isotope, energy in isotopes.items():
+            energy_list.append(
+                {
+                    "shell": shell,
+                    "atomic_number": isotope[0],
+                    "mass_number": isotope[1],
+                    "value": energy,
+                }
+            )
 
-    # DataFrame of positron energies
-    energy_plot_positrons = pd.DataFrame(
-        data=energy_plot_positron_rows,
-        columns=[
-            "packet_index",
-            "energy_input",
-            "energy_input_r",
-            "energy_input_time",
-        ],
+    df = pd.DataFrame(energy_list)
+    energy_df = pd.pivot_table(
+        df,
+        values="value",
+        index=["atomic_number", "mass_number"],
+        columns="shell",
     )
 
-    # DataFrame of estimated deposition
-    # Multiply dataframes by inv_volume_time array
-    # if per unit volume is needed
-    energy_estimated_deposition = (
-        pd.DataFrame(data=deposition_estimator, columns=times[:-1])
-    ) / dt_array
-
-    # Energy is eV/s
-    energy_df = pd.DataFrame(data=energy_df_rows, columns=times[:-1]) / dt_array
-
-    final_energy = 0
-    for p in packets:
-        final_energy += p.energy_rf
-
-    print("Final energy to test for conservation")
-    print(final_energy)
-
-    escape_energy = pd.DataFrame(
-        data=energy_out, columns=times[:-1], index=energy_bins
+    energy_per_mass = energy_df.divide(
+        (raw_isotope_abundance * shell_masses).to_numpy(), axis=0
     )
 
-    return (
-        energy_df,
-        energy_plot_df,
-        escape_energy,
-        decayed_packet_count,
-        energy_plot_positrons,
-        energy_estimated_deposition,
-    )
+    return energy_per_mass, energy_df
diff --git a/tardis/energy_input/tests/test_gamma_ray_transport.py b/tardis/energy_input/tests/test_gamma_ray_transport.py
new file mode 100644
index 00000000000..da17889fe13
--- /dev/null
+++ b/tardis/energy_input/tests/test_gamma_ray_transport.py
@@ -0,0 +1,306 @@
+import os
+import pytest
+import numpy as np
+import pandas as pd
+import numpy.testing as npt
+from pathlib import Path
+import radioactivedecay as rd
+from radioactivedecay import converters
+
+import tardis
+from tardis.model import SimulationState
+from tardis.io.configuration import config_reader
+from tardis.io.util import yaml_load_file, YAMLLoader
+from tardis.energy_input.energy_source import (
+    get_nuclear_lines_database,
+    get_all_isotopes,
+    setup_input_energy,
+)
+from tardis.energy_input.gamma_ray_transport import (
+    calculate_shell_masses,
+    create_isotope_dicts,
+    create_inventories_dict,
+    calculate_total_decays,
+    calculate_average_energies,
+    decay_chain_energies,
+    calculate_energy_per_mass,
+)
+import astropy.units as u
+import astropy.constants as c
+
+
+@pytest.fixture(scope="module")
+def gamma_ray_config(example_configuration_dir: Path):
+    """
+    Parameters
+    ----------
+    example_configuration_dir: Path to the configuration directory.
+
+    Returns
+    -------
+    Tardis configuration
+    """
+    yml_path = (
+        example_configuration_dir
+        / "tardis_configv1_density_exponential_nebular_multi_isotope.yml"
+    )
+
+    return config_reader.Configuration.from_yaml(yml_path)
+
+
+@pytest.fixture(scope="module")
+def simulation_setup(gamma_ray_config):
+    """
+    Parameters
+    ----------
+    gamma_ray_config:
+
+    Returns
+    -------
+    Tardis model
+    """
+    gamma_ray_config.model.structure.velocity.start = 1.0 * u.km / u.s
+    gamma_ray_config.model.structure.density.rho_0 = 5.0e2 * u.g / (u.cm**3)
+    gamma_ray_config.supernova.time_explosion = 1.0 * (u.d)
+    model = SimulationState.from_config(gamma_ray_config)
+    return model
+
+
+def test_calculate_shell_masses(simulation_setup):
+    """Function to test calculation of shell masses.
+    Parameters
+    ----------
+    simulation_setup: A simulation setup which returns a model.
+    """
+    model = simulation_setup
+    volume = 2.70936170e39  # cm^3
+    density = 5.24801665e-09  # g/cm^3
+    desired = volume * density
+
+    shell_masses = calculate_shell_masses(model)[0].value
+    npt.assert_allclose(shell_masses, desired)
+
+
+@pytest.mark.parametrize("nuclide_name", ["Ni-56", "Fe-52", "Cr-48"])
+def test_activity(simulation_setup, nuclide_name):
+    """
+    Function to test the decay of an atom in radioactivedecay with an analytical solution.
+    Parameters
+    ----------
+    simulation_setup: A simulation setup which returns a model.
+    nuclide_name: Name of the nuclide.
+    """
+    nuclide = rd.Nuclide(nuclide_name)
+    model = simulation_setup
+    t_half = nuclide.half_life() * u.s
+    decay_constant = np.log(2) / t_half
+    time_delta = 1.0 * u.s
+    shell_masses = calculate_shell_masses(model)
+    raw_isotope_abundance = model.raw_isotope_abundance
+    raw_isotope_abundance_mass = raw_isotope_abundance.apply(
+        lambda x: x * shell_masses, axis=1
+    )
+    mass = raw_isotope_abundance_mass.loc[nuclide.Z, nuclide.A][0]
+    iso_dict = create_isotope_dicts(raw_isotope_abundance, shell_masses)
+    inv_dict = create_inventories_dict(iso_dict)
+
+    total_decays = calculate_total_decays(inv_dict, time_delta)
+    actual = total_decays[0][nuclide.Z, nuclide.A][nuclide_name]
+
+    isotopic_mass = nuclide.atomic_mass * (u.g)
+    number_of_moles = mass * (u.g) / isotopic_mass
+    number_of_atoms = (number_of_moles * c.N_A).value
+    N1 = number_of_atoms * np.exp(-decay_constant * time_delta)
+    expected = number_of_atoms - N1.value
+
+    npt.assert_allclose(actual, expected)
+
+
+@pytest.mark.parametrize("nuclide_name", ["Ni-56", "Fe-52", "Cr-48"])
+def test_activity_chain(simulation_setup, nuclide_name):
+    """
+    Function to test two atom decay chain in radioactivedecay with an analytical solution.
+    Parameters
+    ----------
+    simulation_setup: A simulation setup which returns a model.
+    nuclide_name: Name of the nuclide.
+    """
+    nuclide = rd.Nuclide(nuclide_name)
+    model = simulation_setup
+    t_half = nuclide.half_life()
+    decay_constant = np.log(2) / t_half
+    time_delta = 1.0 * (u.d).to(u.s)
+    progeny = rd.Nuclide(nuclide.progeny()[0])
+    decay_constant_progeny = np.log(2) / progeny.half_life()
+    shell_masses = calculate_shell_masses(model)
+    raw_isotope_abundance = model.raw_isotope_abundance
+    raw_isotope_abundance_mass = raw_isotope_abundance.apply(
+        lambda x: x * shell_masses, axis=1
+    )
+    mass = raw_isotope_abundance_mass.loc[nuclide.Z, nuclide.A][0]
+    iso_dict = create_isotope_dicts(raw_isotope_abundance, shell_masses)
+    inv_dict = create_inventories_dict(iso_dict)
+    total_decays = calculate_total_decays(inv_dict, time_delta)
+    actual_parent = total_decays[0][nuclide.Z, nuclide.A][nuclide_name]
+    actual_progeny = total_decays[0][nuclide.Z, nuclide.A][nuclide.progeny()[0]]
+    isotopic_mass = nuclide.atomic_mass
+    number_of_moles = mass / isotopic_mass
+    number_of_atoms = number_of_moles * converters.AVOGADRO
+    decayed_parent = number_of_atoms * np.exp(-decay_constant * time_delta)
+    expected_parent = number_of_atoms * (
+        1 - np.exp(-decay_constant * time_delta)
+    )
+
+    npt.assert_almost_equal(expected_parent, actual_parent)
+    # npt.assert_allclose(actual_progeny, expected_progeny, rtol=1e-3)
+
+
+@pytest.mark.parametrize("nuclide_name", ["Ni-56", "Fe-52", "Cr-48"])
+def test_isotope_dicts(simulation_setup, nuclide_name):
+    """
+    Function to test if the right names for the isotopes are present as dictionary keys.
+    Parameters
+    ----------
+    simulation_setup: A simulation setup which returns a model.
+    nuclide_name: Name of the nuclide.
+    """
+    model = simulation_setup
+    nuclide = rd.Nuclide(nuclide_name)
+    raw_isotope_abundance = model.raw_isotope_abundance
+    shell_masses = calculate_shell_masses(model)
+    iso_dict = create_isotope_dicts(raw_isotope_abundance, shell_masses)
+
+    Z, A = nuclide.Z, nuclide.A
+
+    for shell_number, isotope_dict in iso_dict.items():
+        isotope_dict_key = isotope_dict[Z, A]
+        assert nuclide_name.replace("-", "") in isotope_dict_key.keys()
+
+
+@pytest.mark.parametrize("nuclide_name", ["Ni-56", "Fe-52", "Cr-48"])
+def test_inventories_dict(simulation_setup, nuclide_name):
+    """
+    Function to test if the inventories dictionary is created correctly.
+    Parameters
+    ----------
+    simulation_setup: A simulation setup which returns a model.
+    nuclide_name: Name of the nuclide.
+    """
+    model = simulation_setup
+    nuclide = rd.Nuclide(nuclide_name)
+    raw_isotope_abundance = model.raw_isotope_abundance
+    shell_masses = calculate_shell_masses(model)
+    iso_dict = create_isotope_dicts(raw_isotope_abundance, shell_masses)
+    inventories_dict = create_inventories_dict(iso_dict)
+
+    Z, A = nuclide.Z, nuclide.A
+    raw_isotope_abundance_mass = raw_isotope_abundance.apply(
+        lambda x: x * shell_masses, axis=1
+    )
+
+    mass = raw_isotope_abundance_mass.loc[Z, A][0]
+    inventory = rd.Inventory({nuclide.nuclide: mass}, "g")
+    assert inventories_dict[0][Z, A] == inventory
+
+
+def test_average_energies(simulation_setup, atomic_dataset):
+    """
+    Function to test if the energy from each isotope is there in the list.
+    Parameters
+    ----------
+    simulation_setup: A simulation setup which returns a model.
+    atomic_dataset: Tardis atomic and nuclear dataset.
+    """
+
+    model = simulation_setup
+    raw_isotope_abundance = model.raw_isotope_abundance
+    gamma_ray_lines = atomic_dataset.decay_radiation_data
+
+    all_isotope_names = get_all_isotopes(raw_isotope_abundance)
+
+    average_energies_list = []
+
+    for isotope_name in all_isotope_names:
+        energy, intensity = setup_input_energy(
+            gamma_ray_lines[
+                gamma_ray_lines.index == isotope_name.replace("-", "")
+            ],
+            "g",
+        )
+        average_energies_list.append(np.sum(energy * intensity))  # keV
+
+    assert len(average_energies_list) == len(all_isotope_names)
+
+
+@pytest.mark.parametrize("nuclide_name", ["Ni-56", "Fe-52", "Cr-48"])
+def test_decay_energy_chain(simulation_setup, atomic_dataset, nuclide_name):
+    """
+    This function tests if the decay energy is calculated correctly for a decay chain.
+    Parameters
+    ----------
+    simulation_setup: A simulation setup which returns a model.
+    atomic_dataset: Tardis atomic and nuclear dataset.
+    nuclide_name: Name of the nuclide.
+    """
+    model = simulation_setup
+    nuclide = rd.Nuclide(nuclide_name)
+    raw_isotope_abundance = model.raw_isotope_abundance
+    shell_masses = calculate_shell_masses(model)
+    iso_dict = create_isotope_dicts(raw_isotope_abundance, shell_masses)
+    inventories_dict = create_inventories_dict(iso_dict)
+    gamma_ray_lines = atomic_dataset.decay_radiation_data
+    all_isotope_names = get_all_isotopes(raw_isotope_abundance)
+    Z, A = nuclide.Z, nuclide.A
+
+    total_decays = calculate_total_decays(inventories_dict, 1.0 * u.s)
+
+    (
+        average_energies,
+        average_positron_energies,
+        gamma_ray_line_dict,
+    ) = calculate_average_energies(raw_isotope_abundance, gamma_ray_lines)
+
+    decay_chain_energy = decay_chain_energies(
+        average_energies,
+        total_decays,
+    )
+
+    expected = (
+        total_decays[0][Z, A][nuclide_name] * average_energies[nuclide_name]
+    )
+    actual = decay_chain_energy[0][Z, A][nuclide_name]
+
+    npt.assert_almost_equal(expected, actual)
+
+
+def test_energy_per_mass(simulation_setup, atomic_dataset):
+    """
+    This function tests if the energy per mass has the same dimensions as the raw_isotope_abundance.
+    This means for each decay chain we are calculating the energy per mass, by summing the energy from each isotope.
+    Parameters
+    ----------
+    simulation_setup: A simulation setup which returns a model.
+    atomic_dataset: Tardis atomic and nuclear dataset.
+    """
+    model = simulation_setup
+    raw_isotope_abundance = model.raw_isotope_abundance
+    shell_masses = calculate_shell_masses(model)
+    iso_dict = create_isotope_dicts(raw_isotope_abundance, shell_masses)
+    inventories_dict = create_inventories_dict(iso_dict)
+    total_decays = calculate_total_decays(inventories_dict, 1.0 * u.s)
+
+    gamma_ray_lines = atomic_dataset.decay_radiation_data
+    average_energies = calculate_average_energies(
+        raw_isotope_abundance, gamma_ray_lines
+    )
+    decay_energy = decay_chain_energies(
+        average_energies[0],
+        total_decays,
+    )
+    energy_per_mass = calculate_energy_per_mass(
+        decay_energy, raw_isotope_abundance, shell_masses
+    )
+    # If the shape is not same that means the code is not working with multiple isotopes
+    assert (
+        energy_per_mass[0].shape == (raw_isotope_abundance * shell_masses).shape
+    )
diff --git a/tardis/io/atom_data/base.py b/tardis/io/atom_data/base.py
index a1cf9e8d449..1c21373a184 100644
--- a/tardis/io/atom_data/base.py
+++ b/tardis/io/atom_data/base.py
@@ -139,6 +139,7 @@ class AtomData(object):
         "yg_data",
         "two_photon_data",
         "linelist",
+        "decay_radiation_data",
     ]
 
     # List of tuples of the related dataframes.
@@ -275,6 +276,7 @@ def __init__(
         yg_data=None,
         two_photon_data=None,
         linelist=None,
+        decay_radiation_data=None,
     ):
         self.prepared = False
 
@@ -337,6 +339,8 @@ def __init__(
         if linelist is not None:
             self.linelist = linelist
 
+        if decay_radiation_data is not None:
+            self.decay_radiation_data = decay_radiation_data
         self._check_related()
 
         self.symbol2atomic_number = OrderedDict(
diff --git a/tardis/io/configuration/tests/data/tardis_configv1_density_exponential_nebular_multi_isotope.yml b/tardis/io/configuration/tests/data/tardis_configv1_density_exponential_nebular_multi_isotope.yml
new file mode 100644
index 00000000000..4786784264a
--- /dev/null
+++ b/tardis/io/configuration/tests/data/tardis_configv1_density_exponential_nebular_multi_isotope.yml
@@ -0,0 +1,52 @@
+tardis_config_version: v1.0
+
+supernova:
+    luminosity_requested: 9.44 log_lsun
+    time_explosion: 150 day
+
+atom_data: kurucz_atom_pure.h5
+
+model:
+    structure:
+        type: specific
+
+        velocity:
+            start: 1.1e4 km/s
+            stop: 20000 km/s
+            num: 20
+
+        density:
+            type: exponential
+            time_0: 20. second
+            rho_0: 3.e2 g/cm^3
+            v_0: 3000.  km/s
+
+    abundances:
+        type: uniform
+        O: 0.19
+        Mg: 0.03
+        Si: 0.12
+        S: 0.09
+        Ar: 0.02
+        Ca: 0.03
+        Ni56: 0.5
+        Fe52: 0.01
+        Cr48: 0.01
+
+plasma:
+    ionization: nebular
+    excitation: dilute-lte
+    radiative_rates_type: dilute-blackbody
+    line_interaction_type: scatter
+
+montecarlo:
+    seed: 23111963
+    no_of_packets: 2.0e+5
+    iterations: 30
+    last_no_of_packets: 5.0e+5
+    no_of_virtual_packets: 5
+
+spectrum:
+    start: 500 angstrom
+    stop: 20000 angstrom
+    num: 10000
diff --git a/tardis/io/decay.py b/tardis/io/decay.py
index 7005fa50e6e..5b652f460a7 100644
--- a/tardis/io/decay.py
+++ b/tardis/io/decay.py
@@ -57,7 +57,7 @@ def id_to_tuple(atomic_id):
         nuclide = Nuclide(atomic_id)
         return nuclide.Z, nuclide.A
 
-    def to_inventories(self):
+    def to_inventories(self, shell_masses=None):
         """
         Convert DataFrame to a list of inventories interpreting the MultiIndex as
         atomic_number and mass_number
@@ -72,7 +72,12 @@ def to_inventories(self):
         for (atomic_number, mass_number), abundances in self.iterrows():
             nuclear_symbol = f"{Z_to_elem(atomic_number)}{mass_number}"
             for i in range(len(self.columns)):
-                comp_dicts[i][nuclear_symbol] = abundances[i]
+                if shell_masses is None:
+                    comp_dicts[i][nuclear_symbol] = abundances[i]
+                else:
+                    comp_dicts[i][nuclear_symbol] = (
+                        abundances[i] * shell_masses[i].to(u.g).value
+                    )
         return [Inventory(comp_dict, "g") for comp_dict in comp_dicts]
 
     def decay(self, t):