Gamma ray spectra dep

.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>
 

tardis/energy_input/GXPacket.py

@@ -32,6 +32,7 @@
     ("status", int64),
     ("shell", int64),
     ("time_current", float64),
+    ("positron_energy", float64),
     ("tau", float64),
 ]
 
@@ -53,6 +54,7 @@
         status,
         shell,
         time_current,
+        positron_energy,
     ):
         self.location = location
         self.direction = direction
@@ -64,6 +66,7 @@
         self.shell = shell
         self.time_current = time_current
         # TODO: rename to tau_event
+        self.positron_energy = positron_energy
         self.tau = -np.log(np.random.random())
 
     def get_location_r(self):
@@ -95,6 +98,7 @@
         status,
         shell,
         time_current,
+        positron_energy,
     ):
         self.location = location
         self.direction = direction
@@ -106,6 +110,7 @@
         self.shell = shell
         self.time_current = time_current
         self.number_of_packets = len(self.energy_rf)
+        self.positron_energy = positron_energy
         self.tau = -np.log(np.random.random(self.number_of_packets))
 
 

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 @@
     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,
+    energy_deposited_gamma,
+    energy_deposited_positron,
     packets_info_array,
 ):
     """Propagates packets through the simulation
@@ -103,21 +101,23 @@
     escaped_packets = 0
     scattered_packets = 0
     packet_count = len(packets)
+    # Logging does not work with numba
     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)
+        energy_deposited_positron[
+            packet.shell, time_index
+        ] += packet.positron_energy
 
         if time_index < 0:
             print(packet.time_current, 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
@@ -214,17 +214,6 @@
 
             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
 
@@ -246,26 +235,9 @@
 
                 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,
-                    ]
-                )
+                energy_deposited_gamma[
+                    packet.shell, time_index
+                ] += ejecta_energy_gained
 
                 if packet.status == GXPacketStatus.PHOTOABSORPTION:
                     # Packet destroyed, go to the next packet
@@ -279,14 +251,18 @@
 
                 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,7 +279,7 @@
                     packet.nu_cmf,
                     packet.nu_rf,
                     packet.energy_cmf,
-                    lum_rf,
+                    luminosity,
                     packet.energy_rf,
                     packet.shell,
                 ]
@@ -313,12 +289,10 @@
     print("Scattered packets:", scattered_packets)
 
     return (
-        energy_df_rows,
-        energy_plot_df_rows,
         energy_out,
-        deposition_estimator,
-        bin_width,
         packets_info_array,
+        energy_deposited_gamma,
+        energy_deposited_positron,
     )
 
 

tardis/energy_input/gamma_ray_channel.py

@@ -5,6 +5,7 @@
 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)
@@ -83,10 +84,10 @@
     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 = {}
 
@@ -167,3 +168,91 @@
     )
 
     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 = []
+
+    for time in time_array:
+
+        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_evolved_isotope_mass_fraction = pd.concat(
+        isotope_mass_fraction_list, keys=time_array, names=["time"]
+    )
+
+    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)
+    # Do not append dataframes to empty lists
+    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_evolved_decay_df = pd.concat(
+        isotope_decay_df_list, keys=time_array, names=["time"]
+    )
+
+    return time_evolved_decay_df