sampling positronium

tardis/energy_input/gamma_packet_loop.py

@@ -1,7 +1,16 @@
+import logging
 import numpy as np
 from numba import njit
 
-from tardis.energy_input.gamma_ray_estimators import deposition_estimator_kasen
+from tardis.transport.montecarlo import njit_dict_no_parallel
+from tardis.opacities.opacities import (
+    compton_opacity_calculation,
+    photoabsorption_opacity_calculation,
+    pair_creation_opacity_calculation,
+    kappa_calculation,
+    pair_creation_opacity_artis,
+    SIGMA_T,
+)
 from tardis.energy_input.gamma_ray_grid import (
     distance_trace,
     move_packet,
@@ -27,7 +36,7 @@
     pair_creation_opacity_calculation,
     photoabsorption_opacity_calculation,
 )
-from tardis.transport.montecarlo import njit_dict_no_parallel
+from tardis.energy_input.gamma_ray_estimators import deposition_estimator_kasen
 
 
 @njit(**njit_dict_no_parallel)
@@ -38,16 +47,13 @@
     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,
     effective_time_array,
     energy_bins,
-    energy_df_rows,
-    energy_plot_df_rows,
     energy_out,
     packets_info_array,
 ):
@@ -103,22 +109,18 @@
     escaped_packets = 0
     scattered_packets = 0
     packet_count = len(packets)
-    print("Entering gamma ray loop for " + str(packet_count) + " packets")
 
-    deposition_estimator = np.zeros_like(energy_df_rows)
+    print("Packet count:", packet_count)
 
     for i in range(packet_count):
         packet = packets[i]
-        time_index = get_index(packet.time_current, times)
-
+        time_index = get_index(packet.time_current, effective_time_array)
         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
-
         while packet.status == GXPacketStatus.IN_PROCESS:
             # Get delta-time value for this step
             dt = dt_array[time_index]
@@ -203,9 +205,8 @@
                 outer_velocities,
                 total_opacity,
                 effective_time_array[time_index],
-                times[time_index + 1],
+                effective_time_array[time_index + 1],
             )
-
             distance = min(
                 distance_interaction, distance_boundary, distance_time
             )
@@ -214,17 +215,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,27 +236,6 @@
 
                 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,
-                    ]
-                )
-
                 if packet.status == GXPacketStatus.PHOTOABSORPTION:
                     # Packet destroyed, go to the next packet
                     break
@@ -279,14 +248,16 @@
 
                 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
                     )
+                    luminosity = packet.energy_rf / dt
+
                     packet.status = GXPacketStatus.ESCAPED
                     escaped_packets += 1
                     if scattered:
@@ -303,7 +274,7 @@
                     packet.nu_cmf,
                     packet.nu_rf,
                     packet.energy_cmf,
-                    lum_rf,
+                    luminosity,
                     packet.energy_rf,
                     packet.shell,
                 ]
@@ -313,11 +284,7 @@
     print("Scattered packets:", scattered_packets)
 
     return (
-        energy_df_rows,
-        energy_plot_df_rows,
         energy_out,
-        deposition_estimator,
-        bin_width,
         packets_info_array,
     )
 

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)
@@ -75,18 +76,18 @@
         dictionary of inventories for each shell
 
     time_end : float
-        End time of simulation in days.
+        End time of simulation step in days.
 
 
     Returns
     -------
     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 = {}
 
@@ -109,7 +110,8 @@
 
 def create_isotope_decay_df(cumulative_decay_df, gamma_ray_lines):
     """
-    Function to create a dataframe of isotopes for each shell with their decay mode, number of decays, radiation type,
+    Function to create a dataframe of isotopes for each shell with their decay mode,
+    number of decays, radiation type,
     radiation energy and radiation intensity.
 
     Parameters
@@ -167,3 +169,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)
+
+    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

tardis/energy_input/gamma_ray_grid.py

@@ -50,8 +50,6 @@ def calculate_distance_radial(photon, r_inner, r_outer):
         outer_2 = -1
 
     distances = np.array([inner_1, inner_2, outer_1, outer_2])
-
-    # the correct distance is the shortest positive distance
     distance_list = [i for i in distances if i > 0]
 
     if not distance_list: