Gamma ray spectra dep

.mailmap

@@ -47,6 +47,7 @@ Atharwa Kharkar <[email protected]> atharwa_24 <[email protected]
 Atul Kumar <[email protected]>
 
 Anirban Dutta <[email protected]>
+Anirban Dutta <[email protected]>
 
 Barnabás Barna <[email protected]>
 

tardis/energy_input/GXPacket.py

@@ -31,7 +31,8 @@
     ("nu_cmf", float64),
     ("status", int64),
     ("shell", int64),
-    ("time_current", float64),
+    ("time_start", float64),
+    ("time_index", float64),
     ("tau", float64),
 ]
 
@@ -52,7 +53,8 @@
         nu_cmf,
         status,
         shell,
-        time_current,
+        time_start,
+        time_index,
     ):
         self.location = location
         self.direction = direction
@@ -62,7 +64,8 @@
         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 @@
         nu_cmf,
         status,
         shell,
-        time_current,
+        time_start,
+        time_index,
     ):
         self.location = location
         self.direction = direction
@@ -104,8 +108,8 @@
         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.time_current = time_start
+        self.time_index = time_index
         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,
+    total_energy,
+    energy_deposited_gamma,
     packets_info_array,
 ):
     """Propagates packets through the simulation
@@ -103,21 +101,20 @@
     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 @@
                 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 @@
                 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 @@
                 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 @@
 
                 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 @@
                     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 @@
             doppler_factor = doppler_factor_3d(
                 packet.direction,
                 packet.location,
-                packet.time_current,
+                packet.time_start,
             )
 
             packet.nu_rf = packet.nu_cmf / doppler_factor

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 @@
     ----------
     isotope_dict : Dict
         dictionary of isotopes for each shell with their ``masses``.
+
     Returns
     -------
         inv : Dict
@@ -83,10 +86,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 = {}
 
@@ -125,7 +128,6 @@
         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 @@
     )
 
     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