diff --git a/.mailmap b/.mailmap
index ed00506c00f..8607aa37dc7 100644
--- a/.mailmap
+++ b/.mailmap
@@ -47,6 +47,7 @@ Atharwa Kharkar <atharwa.kharkar19@vit.edu> atharwa_24 <atharwa.kharkar19@vit.ed
 Atul Kumar <kr.atul.atk@gmail.com>
 
 Anirban Dutta <anirbaniamdutta@gmail.com>
+Anirban Dutta <duttaan2@msu.edu>
 
 Barnabás Barna <bbarna@titan.physx.u-szeged.hu>
 
diff --git a/tardis/energy_input/GXPacket.py b/tardis/energy_input/GXPacket.py
index dd8792876bc..d24a803e05a 100644
--- a/tardis/energy_input/GXPacket.py
+++ b/tardis/energy_input/GXPacket.py
@@ -31,7 +31,8 @@ class GXPacketStatus(IntEnum):
     ("nu_cmf", float64),
     ("status", int64),
     ("shell", int64),
-    ("time_current", float64),
+    ("time_start", float64),
+    ("time_index", int64),
     ("tau", float64),
 ]
 
@@ -52,7 +53,8 @@ def __init__(
         nu_cmf,
         status,
         shell,
-        time_current,
+        time_start,
+        time_index,
     ):
         self.location = location
         self.direction = direction
@@ -62,7 +64,8 @@ def __init__(
         self.nu_cmf = nu_cmf
         self.status = status
         self.shell = shell
-        self.time_current = time_current
+        self.time_start = time_start
+        self.time_index = time_index
         # TODO: rename to tau_event
         self.tau = -np.log(np.random.random())
 
@@ -94,7 +97,8 @@ def __init__(
         nu_cmf,
         status,
         shell,
-        time_current,
+        time_start,
+        time_index,
     ):
         self.location = location
         self.direction = direction
@@ -104,9 +108,9 @@ def __init__(
         self.nu_cmf = nu_cmf
         self.status = status
         self.shell = shell
-        self.time_current = time_current
-        self.number_of_packets = len(self.energy_rf)
-        self.tau = -np.log(np.random.random(self.number_of_packets))
+        self.time_start = time_start
+        self.time_index = time_index
+        self.tau = -np.log(np.random.random())
 
 
 # @njit(**njit_dict_no_parallel)
diff --git a/tardis/energy_input/gamma_packet_loop.py b/tardis/energy_input/gamma_packet_loop.py
index 385febda679..35c84f384c4 100644
--- a/tardis/energy_input/gamma_packet_loop.py
+++ b/tardis/energy_input/gamma_packet_loop.py
@@ -1,7 +1,6 @@
 import numpy as np
 from numba import njit
 
-from tardis.energy_input.gamma_ray_estimators import deposition_estimator_kasen
 from tardis.energy_input.gamma_ray_grid import (
     distance_trace,
     move_packet,
@@ -38,17 +37,16 @@ def gamma_packet_loop(
     pair_creation_opacity_type,
     electron_number_density_time,
     mass_density_time,
-    inv_volume_time,
     iron_group_fraction_per_shell,
     inner_velocities,
     outer_velocities,
-    times,
     dt_array,
+    times,
     effective_time_array,
     energy_bins,
-    energy_df_rows,
-    energy_plot_df_rows,
     energy_out,
+    total_energy,
+    energy_deposited_gamma,
     packets_info_array,
 ):
     """Propagates packets through the simulation
@@ -103,21 +101,20 @@ def gamma_packet_loop(
     escaped_packets = 0
     scattered_packets = 0
     packet_count = len(packets)
+    # Logging does not work with numba. Using print instead.
     print("Entering gamma ray loop for " + str(packet_count) + " packets")
 
-    deposition_estimator = np.zeros_like(energy_df_rows)
-
     for i in range(packet_count):
         packet = packets[i]
-        time_index = get_index(packet.time_current, times)
+        time_index = packet.time_index
 
         if time_index < 0:
-            print(packet.time_current, time_index)
+            print(packet.time_start, time_index)
             raise ValueError("Packet time index less than 0!")
 
         scattered = False
-
-        initial_energy = packet.energy_cmf
+        # Not used now. Useful for the deposition estimator.
+        # initial_energy = packet.energy_cmf
 
         while packet.status == GXPacketStatus.IN_PROCESS:
             # Get delta-time value for this step
@@ -129,7 +126,7 @@ def gamma_packet_loop(
                 doppler_factor = doppler_factor_3d(
                     packet.direction,
                     packet.location,
-                    effective_time_array[time_index],
+                    times[time_index],
                 )
 
                 kappa = kappa_calculation(comoving_energy)
@@ -210,21 +207,10 @@ def gamma_packet_loop(
                 distance_interaction, distance_boundary, distance_time
             )
 
-            packet.time_current += distance / C_CGS
+            packet.time_start += distance / C_CGS
 
             packet = move_packet(packet, distance)
 
-            deposition_estimator[packet.shell, time_index] += (
-                (initial_energy * 1000)
-                * distance
-                * (packet.energy_cmf / initial_energy)
-                * deposition_estimator_kasen(
-                    comoving_energy,
-                    mass_density_time[packet.shell, time_index],
-                    iron_group_fraction_per_shell[packet.shell],
-                )
-            )
-
             if distance == distance_time:
                 time_index += 1
 
@@ -234,7 +220,7 @@ def gamma_packet_loop(
                 else:
                     packet.shell = get_index(
                         packet.get_location_r(),
-                        inner_velocities * effective_time_array[time_index],
+                        inner_velocities * times[time_index],
                     )
 
             elif distance == distance_interaction:
@@ -246,47 +232,36 @@ def gamma_packet_loop(
 
                 packet, ejecta_energy_gained = process_packet_path(packet)
 
-                # Save packets to dataframe rows
-                # convert KeV to eV / s / cm^3
-                energy_df_rows[packet.shell, time_index] += (
-                    ejecta_energy_gained * 1000
-                )
-
-                energy_plot_df_rows[i] = np.array(
-                    [
-                        i,
-                        ejecta_energy_gained * 1000
-                        # * inv_volume_time[packet.shell, time_index]
-                        / dt,
-                        packet.get_location_r(),
-                        packet.time_current,
-                        packet.shell,
-                        compton_opacity,
-                        photoabsorption_opacity,
-                        pair_creation_opacity,
-                    ]
-                )
+                # Ejecta gains energy from the packets (gamma-rays)
+                energy_deposited_gamma[
+                    packet.shell, time_index
+                ] += ejecta_energy_gained
+                # Ejecta gains energy from both gamma-rays and positrons
+                total_energy[packet.shell, time_index] += ejecta_energy_gained
 
                 if packet.status == GXPacketStatus.PHOTOABSORPTION:
                     # Packet destroyed, go to the next packet
                     break
-                else:
-                    packet.status = GXPacketStatus.IN_PROCESS
-                    scattered = True
+                packet.status = GXPacketStatus.IN_PROCESS
+                scattered = True
 
             else:
                 packet.shell += shell_change
 
                 if packet.shell > len(mass_density_time[:, 0]) - 1:
                     rest_energy = packet.nu_rf * H_CGS_KEV
-                    lum_rf = (packet.energy_rf * 1.6022e-9) / dt
                     bin_index = get_index(rest_energy, energy_bins)
                     bin_width = (
                         energy_bins[bin_index + 1] - energy_bins[bin_index]
                     )
-                    energy_out[bin_index, time_index] += rest_energy / (
-                        bin_width * dt
+                    freq_bin_width = bin_width / H_CGS_KEV
+                    energy_out[bin_index, time_index] += (
+                        packet.energy_rf
+                        / dt
+                        / freq_bin_width  # Take light crossing time into account
                     )
+
+                    luminosity = packet.energy_rf / dt
                     packet.status = GXPacketStatus.ESCAPED
                     escaped_packets += 1
                     if scattered:
@@ -303,22 +278,20 @@ def gamma_packet_loop(
                     packet.nu_cmf,
                     packet.nu_rf,
                     packet.energy_cmf,
-                    lum_rf,
+                    luminosity,
                     packet.energy_rf,
                     packet.shell,
                 ]
             )
 
-    print("Escaped packets:", escaped_packets)
-    print("Scattered packets:", scattered_packets)
+    print("Number of escaped packets:", escaped_packets)
+    print("Number of scattered packets:", scattered_packets)
 
     return (
-        energy_df_rows,
-        energy_plot_df_rows,
         energy_out,
-        deposition_estimator,
-        bin_width,
         packets_info_array,
+        energy_deposited_gamma,
+        total_energy,
     )
 
 
@@ -355,7 +328,7 @@ def process_packet_path(packet):
             doppler_factor = doppler_factor_3d(
                 packet.direction,
                 packet.location,
-                packet.time_current,
+                packet.time_start,
             )
 
             packet.nu_rf = packet.nu_cmf / doppler_factor
diff --git a/tardis/energy_input/gamma_ray_channel.py b/tardis/energy_input/gamma_ray_channel.py
index 6bbb7a9685a..eae5e5f9ec9 100644
--- a/tardis/energy_input/gamma_ray_channel.py
+++ b/tardis/energy_input/gamma_ray_channel.py
@@ -1,10 +1,12 @@
 import logging
+
+import astropy.units as u
 import numpy as np
 import pandas as pd
-import astropy.units as u
 import radioactivedecay as rd
 
 from tardis.energy_input.util import KEV2ERG
+from tardis.model.matter.decay import IsotopicMassFraction
 
 logger = logging.getLogger(__name__)
 logging.basicConfig(level=logging.INFO)
@@ -51,6 +53,7 @@ def create_inventories_dict(isotope_dict):
     ----------
     isotope_dict : Dict
         dictionary of isotopes for each shell with their ``masses``.
+
     Returns
     -------
         inv : Dict
@@ -83,10 +86,10 @@ def calculate_total_decays(inventories, time_delta):
     cumulative_decay_df : pd.DataFrame
         total decays for x g of isotope for time 't'
     """
-    time_delta = u.Quantity(time_delta, u.s)
+    time_delta = u.Quantity(time_delta, u.d)
     total_decays = {}
     for shell, inventory in inventories.items():
-        total_decays[shell] = inventory.cumulative_decays(time_delta.value)
+        total_decays[shell] = inventory.cumulative_decays(time_delta.value, "d")
 
     flattened_dict = {}
 
@@ -125,7 +128,6 @@ def create_isotope_decay_df(cumulative_decay_df, gamma_ray_lines):
         dataframe of isotopes for each shell with their decay mode, number of decays, radiation type,
         radiation energy and radiation intensity.
     """
-
     gamma_ray_lines = gamma_ray_lines.rename_axis(
         "isotope"
     )  # renaming "Isotope" in nndc to "isotope"
@@ -167,3 +169,98 @@ def create_isotope_decay_df(cumulative_decay_df, gamma_ray_lines):
     )
 
     return isotope_decay_df
+
+
+def time_evolve_mass_fraction(raw_isotope_mass_fraction, time_array):
+    """
+    Function to evolve the mass fraction of isotopes with time.
+
+    Parameters
+    ----------
+    raw_isotope_mass_fraction : pd.DataFrame
+        isotope mass fraction in mass fractions.
+    time_array : np.array
+        array of time in days.
+
+    Returns
+    -------
+    time_evolved_isotope_mass_fraction : pd.DataFrame
+        time evolved mass fraction of isotopes.
+    """
+    initial_isotope_mass_fraction = raw_isotope_mass_fraction
+    isotope_mass_fraction_list = []
+    dt = np.diff(time_array)
+    t_start = time_array[:-1]
+    t_end = time_array[1:]
+    t_start_index = np.indices(t_start.shape)[0]
+
+    for time in dt:
+        decayed_isotope_mass_fraction = IsotopicMassFraction(
+            initial_isotope_mass_fraction
+        ).decay(time)
+        isotope_mass_fraction_list.append(decayed_isotope_mass_fraction)
+        initial_isotope_mass_fraction = decayed_isotope_mass_fraction
+
+    time_keys = list(zip(t_start, t_end, t_start_index))
+
+    time_evolved_isotope_mass_fraction = pd.concat(
+        isotope_mass_fraction_list, keys=time_keys, names=["time_start", "time_end", "time_index"],
+    )
+
+    return time_evolved_isotope_mass_fraction
+
+
+def time_evolve_cumulative_decay(
+    raw_isotope_mass_fraction, shell_masses, gamma_ray_lines, time_array
+):
+    """
+    Function to calculate the total decays for each isotope for each shell at each time step.
+
+    Parameters
+    ----------
+    raw_isotope_mass_fraction : pd.DataFrame
+        isotope abundance in mass fractions.
+    shell_masses : numpy.ndarray
+        shell masses in units of g
+    gamma_ray_lines : pd.DataFrame
+        gamma ray lines from nndc stored as a pandas dataframe.
+    time_array : numpy.ndarray
+        array of time steps in days.
+
+    Returns
+    -------
+    time_evolve_decay_df : pd.DataFrame
+        dataframe of isotopes for each shell with their decay mode, number of decays, radiation type,
+        radiation energy and radiation intensity at each time step.
+
+    """
+    isotope_decay_df_list = []
+    initial_isotope_mass_fraction = raw_isotope_mass_fraction
+
+    dt = np.diff(time_array)
+    t_start = time_array[:-1]
+    t_end = time_array[1:]
+    t_start_index = np.indices(time_array.shape)[0]
+    # Do not append dataframes to empty lists. This will be modified in future.
+    for time in dt:
+        isotope_dict = create_isotope_dicts(
+            initial_isotope_mass_fraction, shell_masses
+        )
+        inventories = create_inventories_dict(isotope_dict)
+        total_decays = calculate_total_decays(inventories, time)
+        isotope_df_time = create_isotope_decay_df(total_decays, gamma_ray_lines)
+        isotope_decay_df_list.append(isotope_df_time)
+
+        decayed_isotope_mass_fraction = IsotopicMassFraction(
+            initial_isotope_mass_fraction
+        ).decay(time)
+
+        initial_isotope_mass_fraction = decayed_isotope_mass_fraction
+
+    time_keys = list(zip(t_start, t_end, t_start_index))
+
+    time_evolved_decay_df = pd.concat(
+        isotope_decay_df_list, keys=time_keys, names=["time_start", "time_end", "time_index"],
+    )
+
+    return time_evolved_decay_df
diff --git a/tardis/energy_input/gamma_ray_grid.py b/tardis/energy_input/gamma_ray_grid.py
index 1f86be634af..4321161a042 100644
--- a/tardis/energy_input/gamma_ray_grid.py
+++ b/tardis/energy_input/gamma_ray_grid.py
@@ -1,12 +1,12 @@
 import numpy as np
 from numba import njit
 
-from tardis.transport.montecarlo import njit_dict_no_parallel
 from tardis.energy_input.util import (
+    C_CGS,
     doppler_factor_3d,
     solve_quadratic_equation,
-    C_CGS,
 )
+from tardis.transport.montecarlo import njit_dict_no_parallel
 
 
 @njit(**njit_dict_no_parallel)
@@ -25,7 +25,6 @@ def calculate_distance_radial(photon, r_inner, r_outer):
     distance : float
 
     """
-
     # solve the quadratic distance equation for the inner and
     # outer shell boundaries
     inner_1, inner_2 = solve_quadratic_equation(
@@ -108,7 +107,7 @@ def distance_trace(
     )
 
     distance_interaction = photon.tau / total_opacity
-    distance_time = (next_time - photon.time_current) * C_CGS
+    distance_time = (next_time - photon.time_start) * C_CGS
     return distance_interaction, distance_boundary, distance_time, shell_change
 
 
@@ -135,7 +134,7 @@ def move_packet(packet, distance):
     packet.location = location_new
 
     doppler_factor = doppler_factor_3d(
-        packet.direction, packet.location, packet.time_current
+        packet.direction, packet.location, packet.time_start
     )
 
     packet.nu_cmf = packet.nu_rf * doppler_factor
diff --git a/tardis/energy_input/gamma_ray_interactions.py b/tardis/energy_input/gamma_ray_interactions.py
index 649a526af72..46a9d2ad3f2 100644
--- a/tardis/energy_input/gamma_ray_interactions.py
+++ b/tardis/energy_input/gamma_ray_interactions.py
@@ -197,7 +197,7 @@ def compton_scatter(photon, compton_angle):
     """
     # get comoving frame direction
     comov_direction = angle_aberration_gamma(
-        photon.direction, photon.location, photon.time_current
+        photon.direction, photon.location, photon.time_start
     )
 
     # compute an arbitrary perpendicular vector to the comoving direction
@@ -231,7 +231,7 @@ def compton_scatter(photon, compton_angle):
 
     # Calculate the angle aberration of the final direction
     final_direction = angle_aberration_gamma(
-        final_compton_scattered_vector, photon.location, photon.time_current
+        final_compton_scattered_vector, photon.location, photon.time_start
     )
 
     return final_direction
@@ -265,13 +265,13 @@ def pair_creation_packet(packet):
 
     # Calculate aberration of the random angle for the rest frame
     final_direction = angle_aberration_gamma(
-        new_direction, packet.location, -1 * packet.time_current
+        new_direction, packet.location, -1 * packet.time_start
     )
 
     packet.direction = final_direction
 
     doppler_factor = doppler_factor_3d(
-        packet.direction, packet.location, packet.time_current
+        packet.direction, packet.location, packet.time_start
     )
 
     packet.nu_cmf = ELECTRON_MASS_ENERGY_KEV / H_CGS_KEV
diff --git a/tardis/energy_input/gamma_ray_packet_source.py b/tardis/energy_input/gamma_ray_packet_source.py
index c880bdf3632..4285ed2cbe5 100644
--- a/tardis/energy_input/gamma_ray_packet_source.py
+++ b/tardis/energy_input/gamma_ray_packet_source.py
@@ -1,21 +1,25 @@
+import astropy.units as u
 import numpy as np
-import pandas as pd
 
-from tardis.energy_input.energy_source import (
-    positronium_continuum,
-)
 from tardis.energy_input.GXPacket import (
     GXPacketCollection,
 )
-from tardis.energy_input.samplers import sample_energy
+from tardis.energy_input.samplers import (
+    PositroniumSampler,
+    sample_energy,
+)
 from tardis.energy_input.util import (
     H_CGS_KEV,
     doppler_factor_3d,
+    doppler_factor_3D_all_packets,
     get_index,
     get_random_unit_vector,
 )
 from tardis.transport.montecarlo.packet_source import BasePacketSource
 
+POSITRON_ANNIHILATION_LINE = 511.0
+PARA_TO_ORTHO_RATIO = 0.25
+
 
 class RadioactivePacketSource(BasePacketSource):
     def __init__(
@@ -487,34 +491,56 @@ class GammaRayPacketSource(BasePacketSource):
     def __init__(
         self,
         packet_energy,
-        gamma_ray_lines,
+        isotope_decay_df,
         positronium_fraction,
         inner_velocities,
         outer_velocities,
         inv_volume_time,
         times,
-        energy_df_rows,
         effective_times,
         taus,
         parents,
-        average_positron_energies,
-        average_power_per_mass,
         **kwargs,
     ):
+        """
+        New Gamma ray packet source class
+
+        Parameters
+        ----------
+        packet_energy : float
+            Energy of the gamma ray packet
+        isotope_decay_df : pd.DataFrame
+            DataFrame of isotope decay data
+        positronium_fraction : float
+            Fraction of positrons that form positronium
+        inner_velocities : array
+            Array of inner shell velocities
+        outer_velocities : array
+            Array of outer shell velocities
+        inv_volume_time : array
+            Array of inverse volume times
+            1 / ((4 * np.pi)/3 * (vt) ** 3)
+            Indicates how the ejecta volume changes with time
+        times : array
+            Array of time steps
+        effective_times : array
+            Array of effective time steps
+        taus : dict
+            Dictionary of isotope mean lifetimes in seconds
+        parents : dict
+            Dictionary of isotope parents
+
+        """
         self.packet_energy = packet_energy
-        self.gamma_ray_lines = gamma_ray_lines
+        self.isotope_decay_df = isotope_decay_df
         self.positronium_fraction = positronium_fraction
         self.inner_velocities = inner_velocities
         self.outer_velocities = outer_velocities
         self.inv_volume_time = inv_volume_time
         self.times = times
-        self.energy_df_rows = energy_df_rows
         self.effective_times = effective_times
         self.taus = taus
         self.parents = parents
-        self.average_positron_energies = average_positron_energies
-        self.average_power_per_mass = average_power_per_mass
-        self.energy_plot_positron_rows = np.empty(0)
         super().__init__(**kwargs)
 
     def create_packet_mus(self, no_of_packets, *args, **kwargs):
@@ -545,11 +571,9 @@ def create_packet_radii(self, sampled_packets_df):
 
     def create_packet_nus(
         self,
-        no_of_packets,
         packets,
         positronium_fraction,
-        positronium_energy,
-        positronium_intensity,
+        number_of_packets,
     ):
         """Create an array of packet frequency-energies (i.e. E = h * nu)
 
@@ -561,30 +585,40 @@ def create_packet_nus(
             DataFrame of packets
         positronium_fraction : float
             The fraction of positrons that form positronium
-        positronium_energy : array
-            Array of positronium frequency-energies to sample
-        positronium_intensity : array
-            Array of positronium intensities to sample
+            default is 0.0
 
         Returns
         -------
         array
             Array of sampled frequency-energies
         """
-        energy_array = np.zeros(no_of_packets)
-        zs = np.random.random(no_of_packets)
-        for i in range(no_of_packets):
-            # positron
-            if packets.iloc[i]["decay_type"] == "bp":
-                # positronium formation 75% of the time if fraction is 1
-                if zs[i] < positronium_fraction and np.random.random() < 0.75:
-                    energy_array[i] = sample_energy(
-                        positronium_energy, positronium_intensity
-                    )
-                else:
-                    energy_array[i] = 511
-            else:
-                energy_array[i] = packets.iloc[i]["radiation_energy_kev"]
+        energy_array = np.zeros(number_of_packets)
+
+        all_packets = np.array([True] * number_of_packets)
+
+        # positronium formation if fraction is greater than zero
+        positronium_formation = (
+            np.random.uniform(0, 1, number_of_packets) < positronium_fraction
+        )
+        # annihilation line of positrons
+        annihilation_line = packets["radiation_energy_keV"] == POSITRON_ANNIHILATION_LINE
+        # three photon decay of positronium
+        three_photon_decay = np.random.random(number_of_packets) > PARA_TO_ORTHO_RATIO
+
+        energy_array[all_packets] = packets.loc[
+            all_packets, "radiation_energy_keV"
+        ]
+
+        energy_array[
+            positronium_formation & annihilation_line & three_photon_decay
+        ] = PositroniumSampler().sample_energy(
+            n_samples=np.sum(
+                positronium_formation & annihilation_line & three_photon_decay
+            )
+        )
+        energy_array[
+            positronium_formation & annihilation_line & ~three_photon_decay
+        ] = POSITRON_ANNIHILATION_LINE
 
         return energy_array
 
@@ -704,61 +738,46 @@ def create_packets(
         packet_energies_cmf = np.zeros(number_of_packets)
         nus_rf = np.zeros(number_of_packets)
         nus_cmf = np.zeros(number_of_packets)
-        times = np.zeros(number_of_packets)
-        # set packets to IN_PROCESS status
         statuses = np.ones(number_of_packets, dtype=np.int64) * 3
 
-        self.energy_plot_positron_rows = np.zeros((number_of_packets, 4))
-
-        # compute positronium continuum
-        positronium_energy, positronium_intensity = positronium_continuum()
+        # sample packets from the gamma-ray lines only (include X-rays!)
+        sampled_packets_df_gamma = decays_per_isotope[
+            decays_per_isotope["radiation"] == "g"
+        ]
 
-        # sample packets from dataframe, returning a dataframe where each row is
-        # a sampled packet
-        sampled_packets_df = decays_per_isotope.sample(
+        # sample packets from the time evolving dataframe
+        sampled_packets_df = sampled_packets_df_gamma.sample(
             n=number_of_packets,
             weights="decay_energy_erg",
             replace=True,
             random_state=np.random.RandomState(self.base_seed),
         )
-        # get unique isotopes that have produced packets
-        isotopes = pd.unique(sampled_packets_df.index.get_level_values(2))
 
-        # compute the positron fraction for unique isotopes
-        isotope_positron_fraction = self.calculate_positron_fraction(isotopes)
-
-        # get the packet shell index
-        shells = sampled_packets_df.index.get_level_values(1)
+        # get the isotopes and shells of the sampled packets
+        isotopes = sampled_packets_df.index.get_level_values("isotope")
+        isotope_positron_fraction = self.calculate_positron_fraction(
+            isotopes, number_of_packets
+        )
+        shells = sampled_packets_df.index.get_level_values("shell_number")
 
         # get the inner and outer velocity boundaries for each packet to compute
-        # the initial radii
         sampled_packets_df["inner_velocity"] = self.inner_velocities[shells]
         sampled_packets_df["outer_velocity"] = self.outer_velocities[shells]
 
-        # sample radii at time = 0
+        # The radii of the packets at what ever time they are emitted
         initial_radii = self.create_packet_radii(sampled_packets_df)
 
         # get the time step index of the packets
-        initial_time_indexes = sampled_packets_df.index.get_level_values(0)
+        decay_time_indices = sampled_packets_df.index.get_level_values("time_index")
 
-        # get the time of the middle of the step for each packet
-        packet_effective_times = np.array(
-            [self.effective_times[i] for i in initial_time_indexes]
-        )
-
-        # packet decay time
-        times = self.create_packet_times_uniform_energy(
-            number_of_packets,
-            sampled_packets_df.index.get_level_values(2),
-            packet_effective_times,
-        )
+        effective_decay_times = self.times[decay_time_indices]
 
         # scale radius by packet decay time. This could be replaced with
         # Geometry object calculations. Note that this also adds a random
         # unit vector multiplication for 3D. May not be needed.
         locations = (
-            initial_radii
-            * packet_effective_times
+            initial_radii.values
+            * effective_decay_times
             * self.create_packet_directions(number_of_packets)
         )
 
@@ -768,55 +787,25 @@ def create_packets(
         # the individual gamma-ray energy that makes up a packet
         # co-moving frame, including positronium formation
         nu_energies_cmf = self.create_packet_nus(
-            number_of_packets,
             sampled_packets_df,
             self.positronium_fraction,
-            positronium_energy,
-            positronium_intensity,
+            number_of_packets,
         )
 
-        # equivalent frequencies
         nus_cmf = nu_energies_cmf / H_CGS_KEV
 
-        # per packet co-moving frame total energy
         packet_energies_cmf = self.create_packet_energies(
             number_of_packets, self.packet_energy
         )
-
-        # rest frame gamma-ray energy and frequency
-        # this probably works fine without the loop
-        # non-relativistic
         packet_energies_rf = np.zeros(number_of_packets)
         nus_rf = np.zeros(number_of_packets)
-        for i in range(number_of_packets):
-            doppler_factor = doppler_factor_3d(
-                directions[:, i],
-                locations[:, i],
-                times[i],
-            )
-            packet_energies_rf[i] = packet_energies_cmf[i] / doppler_factor
-            nus_rf[i] = nus_cmf[i] / doppler_factor
-
-            # deposit positron energy in both output arrays
-            # this is an average across all packets that are created
-            # it could be changed to be only for packets that are from positrons
-            self.energy_plot_positron_rows[i] = np.array(
-                [
-                    i,
-                    isotope_positron_fraction[sampled_packets_df["isotopes"][i]]
-                    * packet_energies_cmf[i],
-                    # this needs to be sqrt(sum of squares) to get radius
-                    np.linalg.norm(locations[i]),
-                    times[i],
-                ]
-            )
 
-            # this is an average across all packets that are created
-            # it could be changed to be only for packets that are from positrons
-            self.energy_df_rows[shells[i], times[i]] += (
-                isotope_positron_fraction[sampled_packets_df["isotopes"][i]]
-                * packet_energies_cmf[i]
-            )
+        doppler_factors = doppler_factor_3D_all_packets(
+            directions, locations, effective_decay_times
+        )
+
+        packet_energies_rf = packet_energies_cmf / doppler_factors
+        nus_rf = nus_cmf / doppler_factors
 
         return GXPacketCollection(
             locations,
@@ -827,29 +816,51 @@ def create_packets(
             nus_cmf,
             statuses,
             shells,
-            times,
-        )
+            effective_decay_times,
+            decay_time_indices,
+        ), isotope_positron_fraction
 
-    def calculate_positron_fraction(self, isotopes):
+    def calculate_positron_fraction(self, isotopes, number_of_packets):
         """Calculate the fraction of energy that an isotope
-        releases as positron kinetic energy
+        releases as positron kinetic energy compared to gamma-ray energy
 
         Parameters
         ----------
         isotopes : array
-            Array of isotope names as strings
+            Array of isotope names as strings. Here each isotope is associated with a packet.
+        number_of_packets : int
+            Number of gamma-ray packets
 
         Returns
         -------
         dict
             Fraction of energy released as positron kinetic energy per isotope
         """
-        positron_fraction = {}
-
-        for isotope in isotopes:
-            isotope_energy = self.gamma_ray_lines[isotope][0, :]
-            isotope_intensity = self.gamma_ray_lines[isotope][1, :]
-            positron_fraction[isotope] = self.average_positron_energies[
-                isotope
-            ] / np.sum(isotope_energy * isotope_intensity)
-        return positron_fraction
+        isotope_positron_fraction = np.zeros(number_of_packets)
+
+        # Find the positron fraction from the zeroth shell of the dataframe
+        # this is because the total positron kinetic energy is the same for all shells
+        shell_number_0 = self.isotope_decay_df[
+            self.isotope_decay_df.index.get_level_values("shell_number") == 0
+        ]
+
+        gamma_decay_df = shell_number_0[shell_number_0["radiation"] == "g"]
+
+        positrons_decay_df = shell_number_0[shell_number_0["radiation"] == "bp"]
+        # Find the total energy released from positrons per isotope from the dataframe
+        positron_energy_per_isotope = positrons_decay_df.groupby("isotope")[
+            "energy_per_channel_keV"
+        ].sum()
+        # Find the total energy released from gamma-ray per isotope from the dataframe
+        # TODO: Can be tested with total energy released from all radiation types
+        gamma_energy_per_isotope = gamma_decay_df.groupby("isotope")[
+            "energy_per_channel_keV"
+        ].sum()
+        # TODO: Possibly move this for loop
+        for i, isotope in enumerate(isotopes):
+            if isotope in positron_energy_per_isotope: # check if isotope is in the dataframe
+                isotope_positron_fraction[i] = (
+                    positron_energy_per_isotope[isotope]
+                    / gamma_energy_per_isotope[isotope]
+                )
+        return isotope_positron_fraction
diff --git a/tardis/energy_input/main_gamma_ray_loop.py b/tardis/energy_input/main_gamma_ray_loop.py
index 21934e948c0..3746f9cf09e 100644
--- a/tardis/energy_input/main_gamma_ray_loop.py
+++ b/tardis/energy_input/main_gamma_ray_loop.py
@@ -4,49 +4,62 @@
 import numpy as np
 import pandas as pd
 
-from tardis.energy_input.energy_source import (
-    get_nuclear_lines_database,
-)
 from tardis.energy_input.gamma_packet_loop import gamma_packet_loop
-from tardis.energy_input.gamma_ray_channel import (
-    calculate_total_decays,
-    create_inventories_dict,
-    create_isotope_dicts,
-)
-
-from tardis.energy_input.gamma_ray_transport import (
-    calculate_total_decays_old,
-    create_isotope_dicts_old,
-    create_inventories_dict_old,
-)
-from tardis.energy_input.gamma_ray_packet_source import RadioactivePacketSource
+from tardis.energy_input.gamma_ray_packet_source import GammaRayPacketSource
 from tardis.energy_input.gamma_ray_transport import (
-    calculate_average_energies,
-    calculate_average_power_per_mass,
     calculate_ejecta_velocity_volume,
-    calculate_energy_per_mass,
-    decay_chain_energies,
-    distribute_packets,
     get_taus,
     iron_group_fraction_per_shell,
 )
 from tardis.energy_input.GXPacket import GXPacket
+from tardis.energy_input.util import get_index, make_isotope_string_tardis_like
 
 logger = logging.getLogger(__name__)
 logging.basicConfig(level=logging.INFO)
 
 
+def get_effective_time_array(time_start, time_end, time_space, time_steps):
+    """
+    Function to get the effective time array for the gamma-ray loop.
+
+    Parameters
+    ----------
+    time_start : float
+        start time in days.
+    time_end : float
+        end time in days.
+    time_space : str
+        linear or log.
+    time_steps : int
+        number of time steps.
+
+    Returns
+    -------
+    times : np.ndarray
+        array of times in secs.
+    effective_time_array : np.ndarray
+        effective time array in secs.
+    """
+    assert time_start < time_end, "time_start must be smaller than time_end!"
+    if time_space == "log":
+        times = np.geomspace(time_start, time_end, time_steps + 1)
+        effective_time_array = np.sqrt(times[:-1] * times[1:])
+    else:
+        times = np.linspace(time_start, time_end, time_steps + 1)
+        effective_time_array = 0.5 * (times[:-1] + times[1:])
+
+    return times, effective_time_array
+
+
 def run_gamma_ray_loop(
     model,
-    plasma,
+    isotope_decay_df,
+    cumulative_decays_df,
     num_decays,
-    time_start,
-    time_end,
-    time_space,
-    time_steps,
+    times,
+    effective_time_array,
     seed,
     positronium_fraction,
-    atom_data,
     spectrum_bins,
     grey_opacity,
     photoabsorption_opacity="tardis",
@@ -54,27 +67,58 @@ def run_gamma_ray_loop(
 ):
     """
     Main loop to determine the gamma-ray propagation through the ejecta.
+
+    Parameters
+    ----------
+    model : tardis.model.Radial1DModel
+            Tardis model object.
+    plasma : tardis.plasma.standard_plasmas.BasePlasma
+            Tardis plasma object.
+    isotope_decay_df : pd.DataFrame
+            DataFrame containing the cumulative decay data.
+    cumulative_decays_df : pd.DataFrame
+            DataFrame containing the time evolving mass fractions.
+    num_decays : int
+            Number of packets to decay.
+    times : np.ndarray
+            Array of times in days.
+    effective_time_array : np.ndarray
+            Effective time array in days.
+    seed : int
+            Seed for the random number generator.
+    positronium_fraction : float
+            Fraction of positronium.
+    spectrum_bins : int
+            Number of spectrum bins.
+    grey_opacity : float
+            Grey opacity.
+    photoabsorption_opacity : str
+            Photoabsorption opacity.
+    pair_creation_opacity : str
+            Pair creation opacity.
+
+
+    Returns
+    -------
+    escape_energy : pd.DataFrame
+            DataFrame containing the energy escaping the ejecta.
+    packets_df_escaped : pd.DataFrame
+            DataFrame containing the packets info that escaped the ejecta.
+
+
     """
     np.random.seed(seed)
-    time_explosion = model.time_explosion.to(u.s).value
+    times = times * u.d.to(u.s)
+    effective_time_array = effective_time_array * u.d.to(u.s)
     inner_velocities = model.v_inner.to("cm/s").value
     outer_velocities = model.v_outer.to("cm/s").value
     ejecta_volume = model.volume.to("cm^3").value
-    number_of_shells = model.no_of_shells
     shell_masses = model.volume * model.density
+    number_of_shells = len(shell_masses)
+    # TODO: decaying upto times[0]. raw_isotope_abundance is possibly not the best name
     raw_isotope_abundance = model.composition.raw_isotope_abundance.sort_values(
         by=["atomic_number", "mass_number"], ascending=False
     )
-    time_start *= u.d.to(u.s)
-    time_end *= u.d.to(u.s)
-
-    assert time_start < time_end, "time_start must be smaller than time_end!"
-    if time_space == "log":
-        times = np.geomspace(time_start, time_end, time_steps + 1)
-        effective_time_array = np.sqrt(times[:-1] * times[1:])
-    else:
-        times = np.linspace(time_start, time_end, time_steps + 1)
-        effective_time_array = 0.5 * (times[:-1] + times[1:])
 
     dt_array = np.diff(times)
 
@@ -84,108 +128,70 @@ def run_gamma_ray_loop(
         1.0 / ejecta_velocity_volume[:, np.newaxis]
     ) / effective_time_array**3.0
 
-    energy_df_rows = np.zeros((number_of_shells, 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[0] > 1:
-            plasma.isotope_number_density.loc[atom_number].update = np.sum(
+    for (
+        atom_number
+    ) in model.composition.isotopic_number_density.index.get_level_values(0):
+        values = model.composition.isotopic_number_density.loc[
+            atom_number
+        ].values
+        if values.shape[1] > 1:
+            model.elemental_number_density.loc[atom_number] = np.sum(
                 values, axis=0
             )
         else:
-            plasma.isotope_number_density.loc[atom_number].update = values
+            model.elemental_number_density.loc[atom_number]= values
 
     # Electron number density
-    electron_number_density = plasma.number_density.mul(
-        plasma.number_density.index, axis=0
+    electron_number_density = model.elemental_number_density.mul(
+        model.elemental_number_density.index,
+        axis=0,
     ).sum()
-    taus, parents = get_taus(raw_isotope_abundance)
-    # inventories = raw_isotope_abundance.to_inventories()
     electron_number = np.array(electron_number_density * ejecta_volume)
+
+    # Evolve electron number and mass density with time
     electron_number_density_time = (
         electron_number[:, np.newaxis] * inv_volume_time
     )
-
-    # Calculate decay chain energies
-
     mass_density_time = shell_masses[:, np.newaxis] * inv_volume_time
-    gamma_ray_lines = atom_data.decay_radiation_data
-    isotope_dict = create_isotope_dicts_old(raw_isotope_abundance, shell_masses)
-    inventories_dict = create_inventories_dict_old(isotope_dict)
-    total_decays = calculate_total_decays_old(
-        inventories_dict, time_end - time_start
-    )
-
-    (
-        average_energies,
-        average_positron_energies,
-        gamma_ray_line_dict,
-    ) = calculate_average_energies(raw_isotope_abundance, gamma_ray_lines)
-
-    decayed_energy = decay_chain_energies(
-        average_energies,
-        total_decays,
-    )
-    energy_per_mass, energy_df = calculate_energy_per_mass(
-        decayed_energy, raw_isotope_abundance, shell_masses
-    )
-    average_power_per_mass = calculate_average_power_per_mass(
-        energy_per_mass, time_end - time_start
-    )
-    number_of_isotopes = plasma.isotope_number_density * ejecta_volume
-    total_isotope_number = number_of_isotopes.sum().sum()
-    decayed_packet_count = num_decays * number_of_isotopes.divide(
-        total_isotope_number, axis=0
-    )
-
-    total_energy = energy_df.sum().sum()
-    energy_per_packet = total_energy / num_decays
-    packets_per_isotope_df = (
-        distribute_packets(decayed_energy, total_energy, num_decays)
-        .round()
-        .fillna(0)
-        .astype(int)
-    )
 
-    total_energy = total_energy * u.eV.to("erg")
+    taus, parents = get_taus(raw_isotope_abundance)
+    # Need to get the strings for the isotopes without the dashes
+    taus = make_isotope_string_tardis_like(taus)
 
-    logger.info(f"Total gamma-ray energy is {total_energy}")
+    gamma_df = isotope_decay_df[isotope_decay_df["radiation"] == "g"]
+    gamma_df_decay_energy_keV = gamma_df["decay_energy_keV"]
+    total_energy_gamma = gamma_df["decay_energy_erg"].sum()
 
-    iron_group_fraction = iron_group_fraction_per_shell(model)
-    number_of_packets = packets_per_isotope_df.sum().sum()
-    logger.info(f"Total number of packets is {number_of_packets}")
-    individual_packet_energy = total_energy / number_of_packets
-    logger.info(f"Energy per packet is {individual_packet_energy}")
+    energy_per_packet = total_energy_gamma / num_decays
 
-    logger.info("Initializing packets")
+    logger.info(f"Total energy in gamma-rays is {total_energy_gamma}")
+    logger.info(f"Energy per packet is {energy_per_packet}")
 
-    packet_source = RadioactivePacketSource(
-        individual_packet_energy,
-        gamma_ray_line_dict,
+    packet_source = GammaRayPacketSource(
+        energy_per_packet,
+        isotope_decay_df,
         positronium_fraction,
         inner_velocities,
         outer_velocities,
         inv_volume_time,
         times,
-        energy_df_rows,
         effective_time_array,
         taus,
         parents,
-        average_positron_energies,
-        average_power_per_mass,
     )
 
-    packet_collection = packet_source.create_packets(packets_per_isotope_df)
+    logger.info("Creating packets")
+    packet_collection, isotope_positron_fraction = packet_source.create_packets(
+        cumulative_decays_df, num_decays, seed
+    )
 
-    energy_df_rows = packet_source.energy_df_rows
-    energy_plot_df_rows = np.zeros((number_of_packets, 8))
+    total_energy = np.zeros((number_of_shells, len(times) - 1))
 
     logger.info("Creating packet list")
     packets = []
-    total_cmf_energy = packet_collection.energy_cmf.sum()
-    total_rf_energy = packet_collection.energy_rf.sum()
-    for i in range(number_of_packets):
-        packet = GXPacket(
+    # This for loop is expensive. Need to rewrite GX packet to handle arrays
+    packets = [
+        GXPacket(
             packet_collection.location[:, i],
             packet_collection.direction[:, i],
             packet_collection.energy_rf[i],
@@ -194,36 +200,40 @@ def run_gamma_ray_loop(
             packet_collection.nu_cmf[i],
             packet_collection.status[i],
             packet_collection.shell[i],
-            packet_collection.time_current[i],
-        )
-        packets.append(packet)
-        energy_plot_df_rows[i] = np.array(
-            [
-                i,
-                packet.energy_rf,
-                packet.get_location_r(),
-                packet.time_current,
-                int(packet.status),
-                0,
-                0,
-                0,
-            ]
+            packet_collection.time_start[i],
+            packet_collection.time_index[i],
         )
+        for i in range(num_decays)
+    ]
 
-    logger.info(f"Total cmf energy is {total_cmf_energy}")
-    logger.info(f"Total rf energy is {total_rf_energy}")
+    for i, p in enumerate(packets):
+        total_energy[p.shell, p.time_index] += isotope_positron_fraction[i] * energy_per_packet
+
+    logger.info(f"Total energy deposited by the positrons is {total_energy.sum().sum()}")
 
     energy_bins = np.logspace(2, 3.8, spectrum_bins)
-    energy_out = np.zeros((len(energy_bins - 1), time_steps))
+    energy_out = np.zeros((len(energy_bins - 1), len(times) - 1))
+    energy_deposited = np.zeros((number_of_shells, len(times) - 1))
     packets_info_array = np.zeros((int(num_decays), 8))
+    iron_group_fraction = iron_group_fraction_per_shell(model)
+
+    logger.info("Entering the main gamma-ray loop")
+
+    total_cmf_energy = 0
+    total_rf_energy = 0
+
+    for p in packets:
+        total_cmf_energy += p.energy_cmf
+        total_rf_energy += p.energy_rf
+
+    logger.info(f"Total CMF energy is {total_cmf_energy}")
+    logger.info(f"Total RF energy is {total_rf_energy}")
 
     (
-        energy_df_rows,
-        energy_plot_df_rows,
         energy_out,
-        deposition_estimator,
-        bin_width,
         packets_array,
+        energy_deposited_gamma,
+        total_energy,
     ) = gamma_packet_loop(
         packets,
         grey_opacity,
@@ -231,58 +241,105 @@ def run_gamma_ray_loop(
         pair_creation_opacity,
         electron_number_density_time,
         mass_density_time,
-        inv_volume_time,
         iron_group_fraction.to_numpy(),
         inner_velocities,
         outer_velocities,
-        times,
         dt_array,
+        times,
         effective_time_array,
         energy_bins,
-        energy_df_rows,
-        energy_plot_df_rows,
         energy_out,
+        total_energy,
+        energy_deposited,
         packets_info_array,
     )
 
-    energy_plot_df = pd.DataFrame(
-        data=energy_plot_df_rows,
+    packets_df_escaped = pd.DataFrame(
+        data=packets_array,
         columns=[
             "packet_index",
-            "energy_input",
-            "energy_input_r",
-            "energy_input_time",
-            "energy_input_type",
-            "compton_opacity",
-            "photoabsorption_opacity",
-            "total_opacity",
+            "status",
+            "nu_cmf",
+            "nu_rf",
+            "energy_cmf",
+            "luminosity",
+            "energy_rf",
+            "shell_number",
         ],
     )
 
-    energy_plot_positrons = pd.DataFrame(
-        data=packet_source.energy_plot_positron_rows,
-        columns=[
-            "packet_index",
-            "energy_input",
-            "energy_input_r",
-            "energy_input_time",
-        ],
-    )
-
-    energy_estimated_deposition = (
-        pd.DataFrame(data=deposition_estimator, columns=times[:-1])
-    ) / dt_array
-
-    energy_df = pd.DataFrame(data=energy_df_rows, columns=times[:-1]) / dt_array
     escape_energy = pd.DataFrame(
         data=energy_out, columns=times[:-1], index=energy_bins
     )
 
+    # deposited energy by gamma-rays in ergs
+    gamma_ray_deposited_energy = pd.DataFrame(
+        data=energy_deposited_gamma, columns=times[:-1]
+    )
+    # deposited energy by positrons and gamma-rays in ergs
+    total_energy = pd.DataFrame(
+        data=total_energy, columns=times[:-1]
+    )
+
+    logger.info(f"Total energy deposited by gamma rays and positrons is {total_energy.sum().sum()}")
+
+    total_deposited_energy = total_energy / dt_array
+
     return (
-        energy_df,
-        energy_plot_df,
         escape_energy,
-        decayed_packet_count,
-        energy_plot_positrons,
-        energy_estimated_deposition,
+        packets_df_escaped,
+        gamma_ray_deposited_energy,
+        total_energy,
+    )
+
+
+def get_packet_properties(number_of_shells, times, time_steps, packets):
+    """
+    Function to get the properties of the packets.
+
+    Parameters
+    ----------
+    packets : list
+        List of packets.
+
+    Returns
+    -------
+    packets_nu_cmf_array : np.ndarray
+        Array of packets in cmf.
+    packets_nu_rf_array : np.ndarray
+        Array of packets in rf.
+    packets_energy_cmf_array : np.ndarray
+        Array of packets energy in cmf.
+    packets_energy_rf_array : np.ndarray
+        Array of packets energy in rf.
+    packets_positron_energy_array : np.ndarray
+        Array of packets positron energy.
+    """
+    shell_number = []
+    time_current = []
+
+    # Bin the frequency of the packets in shell and time
+
+    packets_nu_cmf_array = np.zeros((number_of_shells, time_steps))
+    packets_nu_rf_array = np.zeros((number_of_shells, time_steps))
+    packets_energy_cmf_array = np.zeros((number_of_shells, time_steps))
+    packets_energy_rf_array = np.zeros((number_of_shells, time_steps))
+    packets_positron_energy_array = np.zeros((number_of_shells, time_steps))
+
+    for p in packets:
+        time_index = get_index(p.time_current, times)
+        shell_number.append(p.shell)
+        time_current.append(p.time_current)
+        packets_nu_cmf_array[p.shell, time_index] += p.nu_cmf
+        packets_nu_rf_array[p.shell, time_index] += p.nu_rf
+        packets_energy_cmf_array[p.shell, time_index] += p.energy_cmf
+        packets_energy_rf_array[p.shell, time_index] += p.energy_rf
+        packets_positron_energy_array[p.shell, time_index] += p.positron_energy
+
+    return (
+        packets_nu_cmf_array,
+        packets_nu_rf_array,
+        packets_energy_cmf_array,
+        packets_energy_rf_array,
+        packets_positron_energy_array,
     )
diff --git a/tardis/energy_input/samplers.py b/tardis/energy_input/samplers.py
index 674345c649f..aa2b5baae6f 100644
--- a/tardis/energy_input/samplers.py
+++ b/tardis/energy_input/samplers.py
@@ -1,3 +1,5 @@
+import astropy.constants as const
+import astropy.units as u
 import numpy as np
 from numba import njit
 
@@ -138,3 +140,53 @@ def sample_decay_time(
             np.random.random()
         )
     return decay_time
+
+
+class PositroniumSampler:
+    def __init__(self, n_grid=1000):
+        """
+        Parameters
+        ----------
+        n_grid : int, optional
+            Number of grid points for the CDF, by default 1000
+        """
+        self.x_grid = np.linspace(0.01, 0.99, n_grid)
+        self.cdf_grid = np.array(
+            [
+                np.trapz(self.pdf(self.x_grid[:i]), self.x_grid[:i])
+                for i in range(len(self.x_grid))
+            ]
+        )
+        self.cdf_grid /= self.cdf_grid[-1]
+
+    @staticmethod
+    def pdf(x):
+        """
+        Parameters
+        ----------
+        x : float
+            E / m_e c^2
+
+        Returns
+        -------
+        float
+            PDF of positronium energy from Ore amnd Powell (1949)
+        """
+        first_term = x * (1 - x) / (2 - x) ** 2
+        second_term = 2 * (1 - x) ** 2 * np.log(1 - x) / (2 - x) ** 2
+        third_term = (2 - x) / x
+        fourth_term = 2 * (1 - x) * np.log(1 - x) / x**2
+
+        return 2 * (first_term - second_term + third_term + fourth_term)
+
+    def quantile_function(self, p):
+
+        return np.interp(p, self.cdf_grid, self.x_grid)
+
+    def sample_energy(self, n_samples=1):
+
+        return (
+            self.quantile_function(np.random.random(n_samples))
+            * const.m_e.cgs.value
+            * const.c.cgs.value**2
+        ) * u.erg.to(u.keV)
diff --git a/tardis/energy_input/tests/conftest.py b/tardis/energy_input/tests/conftest.py
index b36623ef660..c9802bf5c9a 100644
--- a/tardis/energy_input/tests/conftest.py
+++ b/tardis/energy_input/tests/conftest.py
@@ -22,5 +22,6 @@ def basic_gamma_ray():
         nu_cmf=1000.0e3 / H_CGS_KEV,
         status=GXPacketStatus.IN_PROCESS,
         shell=1,
-        time_current=1000,
+        time_start=2.0,
+        time_index=0,
     )
diff --git a/tardis/energy_input/tests/test_gamma_ray_channel.py b/tardis/energy_input/tests/test_gamma_ray_channel.py
index 5842ba67783..0ef52df32a5 100644
--- a/tardis/energy_input/tests/test_gamma_ray_channel.py
+++ b/tardis/energy_input/tests/test_gamma_ray_channel.py
@@ -1,23 +1,24 @@
-import pytest
-import numpy as np
 from pathlib import Path
+
+import astropy.constants as const
 import astropy.units as u
+import numpy as np
 import numpy.testing as npt
+import pytest
 import radioactivedecay as rd
-import astropy.constants as const
-from radioactivedecay import converters
 
-from tardis.model import SimulationState
-from tardis.io.configuration import config_reader
-from tardis.energy_input.energy_source import (
-    get_nuclear_lines_database,
-)
 from tardis.energy_input.gamma_ray_channel import (
-    create_isotope_dicts,
-    create_inventories_dict,
     calculate_total_decays,
+    create_inventories_dict,
     create_isotope_decay_df,
+    create_isotope_dicts,
+    time_evolve_cumulative_decay,
 )
+from tardis.energy_input.gamma_ray_transport import get_taus
+from tardis.energy_input.main_gamma_ray_loop import get_effective_time_array
+from tardis.energy_input.util import KEV2ERG
+from tardis.io.configuration import config_reader
+from tardis.model import SimulationState
 
 
 @pytest.fixture(scope="module")
@@ -234,3 +235,124 @@ def test_activity(gamma_ray_test_composition, nuclide_name):
     expected = number_of_atoms * (1 - np.exp(-decay_constant * time_delta))
 
     npt.assert_allclose(actual, expected)
+
+
+def test_total_energy_production(gamma_ray_test_composition, atomic_dataset):
+    """
+    Function to test the total energy production with equation 18 of Nadyozhin 1994. This is only for Ni56 now.
+    Parameters
+    ----------
+    gamma_ray_test_composition: Function holding the composition.
+    """
+    time_start = 0.0 * u.d
+    time_end = np.inf * u.d
+    time_delta = (time_end - time_start).value
+
+    gamma_ray_lines = atomic_dataset.decay_radiation_data
+    raw_isotopic_mass_fraction, cell_masses = gamma_ray_test_composition
+    isotope_dict = create_isotope_dicts(raw_isotopic_mass_fraction, cell_masses)
+    inventories_dict = create_inventories_dict(isotope_dict)
+    total_decays = calculate_total_decays(inventories_dict, time_delta)
+    isotope_decay_df = create_isotope_decay_df(total_decays, gamma_ray_lines)
+    taus, parents = get_taus(raw_isotopic_mass_fraction)
+
+    ni56_nuclide = rd.Nuclide("Ni56")
+    atomic_mass_unit = const.u.cgs.value
+    tau_ni56 = taus["Ni-56"]
+    tau_co56 = taus["Co-56"]
+
+    ni56_mass_fraction = raw_isotopic_mass_fraction.loc[(28, 56)]
+
+    ni_56_mass = sum(ni56_mass_fraction * cell_masses)
+    ni_56_mass_solar = ni_56_mass / const.M_sun.cgs.value
+
+    shell_number_0 = isotope_decay_df[
+        isotope_decay_df.index.get_level_values("shell_number") == 0
+    ]
+
+    ni56 = shell_number_0[shell_number_0.index.get_level_values(1) == "Ni56"]
+    ni_56_energy = ni56["energy_per_channel_keV"].sum()
+
+    co_56 = shell_number_0[shell_number_0.index.get_level_values(1) == "Co56"]
+    co_56_energy = co_56["energy_per_channel_keV"].sum()
+
+    first_term = const.M_sun.cgs.value / (
+        (tau_co56 - tau_ni56) * ni56_nuclide.atomic_mass * atomic_mass_unit
+    )
+    ni_term = (
+        (ni_56_energy * (tau_co56 / tau_ni56 - 1) - co_56_energy)
+        * first_term
+        * KEV2ERG
+    )
+
+    co_term = co_56_energy * first_term * KEV2ERG
+
+    expected = (
+        ni_term
+        * tau_ni56
+        * (
+            np.exp(-time_start.value / tau_ni56)
+            - np.exp(-time_end.value / tau_ni56)
+        )
+        + co_term
+        * tau_co56
+        * (
+            np.exp(-time_start.value / tau_co56)
+            - np.exp(-time_end.value / tau_co56)
+        )
+    ) * ni_56_mass_solar
+
+    ni56_df = isotope_decay_df[
+        isotope_decay_df.index.get_level_values(1) == "Ni56"
+    ]
+    ni56_energy = ni56_df["decay_energy_erg"].sum()
+    co_56_df = isotope_decay_df[
+        isotope_decay_df.index.get_level_values(1) == "Co56"
+    ]
+    co56_energy = co_56_df["decay_energy_erg"].sum()
+    actual = ni56_energy + co56_energy
+
+    npt.assert_allclose(actual, expected)
+
+
+def test_cumulative_decays(gamma_ray_test_composition, atomic_dataset):
+    """
+    Function to test that the total energy calculated from summing all the decays
+    from the entire time range of simulation is the same as decay energy from individual
+    time steps considering that after each time step the composition (mass fractions) changes.
+    Tested for Ni56, Cr48, Fe52.
+    Parameters
+    ----------
+    gamma_ray_simulation_state: Tardis simulation state
+    atomic_dataset: Tardis atomic-nuclear dataset
+    """
+
+    time_start = 0.1 * u.d
+    time_end = 100 * u.d
+    time_steps = 3
+    time_space = "linear"
+    time_delta = (time_end - time_start).value
+
+    gamma_ray_lines = atomic_dataset.decay_radiation_data
+    raw_isotopic_mass_fraction, cell_masses = gamma_ray_test_composition
+    isotope_dict = create_isotope_dicts(raw_isotopic_mass_fraction, cell_masses)
+    inventories_dict = create_inventories_dict(isotope_dict)
+    total_decays = calculate_total_decays(inventories_dict, time_delta)
+    isotope_decay_df = create_isotope_decay_df(total_decays, gamma_ray_lines)
+
+    times, effective_times = get_effective_time_array(
+        time_start.value, time_end.value, time_space, time_steps
+    )
+    # total decay energy in the entire time range
+    actual = isotope_decay_df["decay_energy_erg"].sum()
+
+    # time evolve the decay energy
+    evolve_decays_with_time = time_evolve_cumulative_decay(
+        raw_isotopic_mass_fraction, cell_masses, gamma_ray_lines, times
+    )
+    expected = evolve_decays_with_time["decay_energy_erg"].sum()
+
+    # This rtol is set since the decay energy is calculated with Fe52 (which has Mn-52m as a daughter)
+    # The data is not available for Mn-52m in the decay_radiation_data
+    # If we use any other isotope without a metastable state, the total decay energy matches exactly.
+    npt.assert_allclose(actual, expected, rtol=1e-4)
diff --git a/tardis/energy_input/tests/test_gamma_ray_grid.py b/tardis/energy_input/tests/test_gamma_ray_grid.py
index 11a6f394ae5..269594d219c 100644
--- a/tardis/energy_input/tests/test_gamma_ray_grid.py
+++ b/tardis/energy_input/tests/test_gamma_ray_grid.py
@@ -1,10 +1,7 @@
-import pytest
 import numpy.testing as npt
-import numpy as np
+import pytest
 
 from tardis.energy_input.gamma_ray_grid import (
-    calculate_distance_radial,
-    distance_trace,
     move_packet,
 )
 
diff --git a/tardis/energy_input/tests/test_gamma_ray_interactions.py b/tardis/energy_input/tests/test_gamma_ray_interactions.py
index fedc16d9619..dc96eabb3de 100644
--- a/tardis/energy_input/tests/test_gamma_ray_interactions.py
+++ b/tardis/energy_input/tests/test_gamma_ray_interactions.py
@@ -1,15 +1,13 @@
-import pytest
 import numpy as np
 import numpy.testing as npt
+import pytest
 
 from tardis.energy_input.gamma_ray_interactions import (
-    compton_scatter,
     pair_creation_packet,
     scatter_type,
 )
-from tardis.energy_input.util import ELECTRON_MASS_ENERGY_KEV, H_CGS_KEV
-
 from tardis.energy_input.GXPacket import GXPacketStatus
+from tardis.energy_input.util import ELECTRON_MASS_ENERGY_KEV, H_CGS_KEV
 
 
 @pytest.mark.xfail(reason="To be implemented")
diff --git a/tardis/energy_input/util.py b/tardis/energy_input/util.py
index a75beae4837..61ea173f409 100644
--- a/tardis/energy_input/util.py
+++ b/tardis/energy_input/util.py
@@ -77,6 +77,28 @@ def doppler_factor_3d(direction_vector, position_vector, time):
     )
 
 
+@njit(**njit_dict_no_parallel)
+def doppler_factor_3D_all_packets(direction_vectors, position_vectors, times):
+    """Doppler shift for photons in 3D
+
+    Parameters
+    ----------
+    direction_vectors : array
+    position_vectors : array
+    times : array
+
+    Returns
+    -------
+    array
+        Doppler factors
+    """
+    velocity_vector = position_vectors / times
+    vel_mul_dir = np.multiply(velocity_vector, direction_vectors)
+    doppler_factors = 1 - (np.sum(vel_mul_dir, axis=0) / C_CGS)
+
+    return doppler_factors
+
+
 @njit(**njit_dict_no_parallel)
 def angle_aberration_gamma(direction_vector, position_vector, time):
     """Angle aberration formula for photons in 3D
@@ -169,13 +191,13 @@ def solve_quadratic_equation(position, direction, radius):
     a = np.sum(direction**2)
     b = 2.0 * np.sum(position * direction)
     c = -(radius**2) + np.sum(position**2)
-    root = b**2 - 4 * a * c
+    discriminant = b**2 - 4 * a * c
     solution_1 = -np.inf
     solution_2 = -np.inf
-    if root > 0.0:
-        solution_1 = (-b + np.sqrt(root)) / (2 * a)
-        solution_2 = (-b - np.sqrt(root)) / (2 * a)
-    elif root == 0:
+    if discriminant > 0.0:
+        solution_1 = (-b + np.sqrt(discriminant)) / (2 * a)
+        solution_2 = (-b - np.sqrt(discriminant)) / (2 * a)
+    elif discriminant == 0:
         solution_1 = -b / (2 * a)
 
     return solution_1, solution_2
@@ -203,13 +225,13 @@ def solve_quadratic_equation_expanding(position, direction, time, radius):
     a = np.dot(direction, direction) - (radius / light_distance) ** 2.0
     b = 2.0 * (np.dot(position, direction) - radius**2.0 / light_distance)
     c = np.dot(position, position) - radius**2.0
-    root = b**2.0 - 4.0 * a * c
+    discriminant = b**2.0 - 4.0 * a * c
     solution_1 = -np.inf
     solution_2 = -np.inf
-    if root > 0.0:
-        solution_1 = (-b + np.sqrt(root)) / (2.0 * a)
-        solution_2 = (-b - np.sqrt(root)) / (2.0 * a)
-    elif root == 0:
+    if discriminant > 0.0:
+        solution_1 = (-b + np.sqrt(discriminant)) / (2.0 * a)
+        solution_2 = (-b - np.sqrt(discriminant)) / (2.0 * a)
+    elif discriminant == 0:
         solution_1 = -b / (2.0 * a)
 
     return solution_1, solution_2
@@ -392,3 +414,26 @@ def get_index(value, array):
         i += 1
 
     return i
+
+
+def make_isotope_string_tardis_like(isotope_dict):
+    """Converts isotope string to TARDIS format
+        Ni-56 -> Ni56, Co-56 -> Co56
+    Parameters
+    ----------
+    isotope : str
+        Isotope string
+
+    Returns
+    -------
+    str
+        TARDIS-like isotope string
+    """
+
+    new_isotope_dict = {}
+
+    for key in isotope_dict.keys():
+        new_key = key.replace("-", "")
+        new_isotope_dict[new_key] = isotope_dict[key]
+
+    return new_isotope_dict