Last interaction radius

tardis/montecarlo/base.py

@@ -221,6 +221,7 @@ def _initialize_packets(self, T, no_of_packets, iteration, radius):
         self.input_energy = energies
 
         self._output_nu = np.ones(no_of_packets, dtype=np.float64) * -99.0
+        self._output_r = np.ones(no_of_packets, dtype=np.float64) * -99.0
         self._output_energy = np.ones(no_of_packets, dtype=np.float64) * -99.0
 
         self.last_line_interaction_in_id = -1 * np.ones(
@@ -385,6 +386,10 @@ def get_line_interaction_id(self, line_interaction_type):
     def output_nu(self):
         return u.Quantity(self._output_nu, u.Hz)
 
+    @property
+    def output_r(self):
+        return u.Quantity(self._output_r, u.cm)
+
     @property
     def output_energy(self):
         return u.Quantity(self._output_energy, u.erg)

tardis/montecarlo/montecarlo_numba/base.py

@@ -47,6 +47,7 @@ def montecarlo_radial1d(
         runner.input_mu,
         runner.input_energy,
         runner._output_nu,
+        runner._output_r,
         runner._output_energy,
     )
 
@@ -76,13 +77,15 @@ def montecarlo_radial1d(
         v_packets_energy_hist,
         last_interaction_type,
         last_interaction_in_nu,
+        last_interaction_in_r,
         last_line_interaction_in_id,
         last_line_interaction_out_id,
         virt_packet_nus,
         virt_packet_energies,
         virt_packet_initial_mus,
         virt_packet_initial_rs,
         virt_packet_last_interaction_in_nu,
+        virt_packet_last_interaction_in_r,
         virt_packet_last_interaction_type,
         virt_packet_last_line_interaction_in_id,
         virt_packet_last_line_interaction_out_id,
@@ -104,6 +107,7 @@ def montecarlo_radial1d(
     runner._montecarlo_virtual_luminosity.value[:] = v_packets_energy_hist
     runner.last_interaction_type = last_interaction_type
     runner.last_interaction_in_nu = last_interaction_in_nu
+    runner.last_interaction_in_r = last_interaction_in_r
     runner.last_line_interaction_in_id = last_line_interaction_in_id
     runner.last_line_interaction_out_id = last_line_interaction_out_id
 
@@ -121,6 +125,9 @@ def montecarlo_radial1d(
         runner.virt_packet_last_interaction_in_nu = np.concatenate(
             virt_packet_last_interaction_in_nu
         ).ravel()
+        runner.virt_packet_last_interaction_in_r = np.concatenate(
+            np.array(virt_packet_last_interaction_in_r)
+        ).ravel()
         runner.virt_packet_last_interaction_type = np.concatenate(
             virt_packet_last_interaction_type
         ).ravel()
@@ -170,12 +177,14 @@ def montecarlo_main_loop(
         Option to enable virtual packet logging.
     """
     output_nus = np.empty_like(packet_collection.packets_input_nu)
+    output_rs = np.empty_like(packet_collection.packets_output_r)
     last_interaction_types = (
         np.ones_like(packet_collection.packets_output_nu, dtype=np.int64) * -1
     )
     output_energies = np.empty_like(packet_collection.packets_output_nu)
 
     last_interaction_in_nus = np.empty_like(packet_collection.packets_output_nu)
+    last_interaction_in_rs = np.empty_like(packet_collection.packets_output_r)
     last_line_interaction_in_ids = (
         np.ones_like(packet_collection.packets_output_nu, dtype=np.int64) * -1
     )
@@ -227,6 +236,7 @@ def montecarlo_main_loop(
     virt_packet_initial_mus = []
     virt_packet_initial_rs = []
     virt_packet_last_interaction_in_nu = []
+    virt_packet_last_interaction_in_r = []
     virt_packet_last_interaction_type = []
     virt_packet_last_line_interaction_in_id = []
     virt_packet_last_line_interaction_out_id = []
@@ -268,7 +278,9 @@ def montecarlo_main_loop(
         )
 
         output_nus[i] = r_packet.nu
+        output_rs[i] = r_packet.r
         last_interaction_in_nus[i] = r_packet.last_interaction_in_nu
+        last_interaction_in_rs[i] = r_packet.last_interaction_in_r
         last_line_interaction_in_ids[i] = r_packet.last_line_interaction_in_id
         last_line_interaction_out_ids[i] = r_packet.last_line_interaction_out_id
 
@@ -329,6 +341,11 @@ def montecarlo_main_loop(
                     ]
                 )
             )
+            virt_packet_last_interaction_in_r.append(
+                vpacket_collection.last_interaction_in_r[
+                    : vpacket_collection.idx
+                ]
+            )
             virt_packet_last_interaction_type.append(
                 np.ascontiguousarray(
                     vpacket_collection.last_interaction_type[
@@ -357,17 +374,20 @@ def montecarlo_main_loop(
 
     packet_collection.packets_output_energy[:] = output_energies[:]
     packet_collection.packets_output_nu[:] = output_nus[:]
+    packet_collection.packets_output_r[:] = output_rs[:]
     return (
         v_packets_energy_hist,
         last_interaction_types,
         last_interaction_in_nus,
+        last_interaction_in_rs,
         last_line_interaction_in_ids,
         last_line_interaction_out_ids,
         virt_packet_nus,
         virt_packet_energies,
         virt_packet_initial_mus,
         virt_packet_initial_rs,
         virt_packet_last_interaction_in_nu,
+        virt_packet_last_interaction_in_r,
         virt_packet_last_interaction_type,
         virt_packet_last_line_interaction_in_id,
         virt_packet_last_line_interaction_out_id,

tardis/montecarlo/montecarlo_numba/numba_interface.py

@@ -282,6 +282,7 @@ def numba_plasma_initialize(plasma, line_interaction_type):
     ("packets_input_mu", float64[:]),
     ("packets_input_energy", float64[:]),
     ("packets_output_nu", float64[:]),
+    ("packets_output_r", float64[:]),
     ("packets_output_energy", float64[:]),
 ]
 
@@ -295,13 +296,15 @@ def __init__(
         packets_input_mu,
         packets_input_energy,
         packets_output_nu,
+        packets_output_r,
         packets_output_energy,
     ):
         self.packets_input_radius = packets_input_radius
         self.packets_input_nu = packets_input_nu
         self.packets_input_mu = packets_input_mu
         self.packets_input_energy = packets_input_energy
         self.packets_output_nu = packets_output_nu
+        self.packets_output_r = packets_output_r
         self.packets_output_energy = packets_output_energy
 
 
@@ -318,6 +321,7 @@ def __init__(
     ("number_of_vpackets", int64),
     ("length", int64),
     ("last_interaction_in_nu", float64[:]),
+    ("last_interaction_in_r", float64[:]),
     ("last_interaction_type", int64[:]),
     ("last_interaction_in_id", int64[:]),
     ("last_interaction_out_id", int64[:]),
@@ -346,6 +350,9 @@ def __init__(
         self.last_interaction_in_nu = np.zeros(
             temporary_v_packet_bins, dtype=np.float64
         )
+        self.last_interaction_in_r = np.zeros(
+            temporary_v_packets_bins, dtype=np.float64
+        )
         self.last_interaction_type = -1 * np.ones(
             temporary_v_packet_bins, dtype=np.int64
         )
@@ -366,6 +373,7 @@ def set_properties(
         initial_mu,
         initial_r,
         last_interaction_in_nu,
+        last_interaction_in_r,
         last_interaction_type,
         last_interaction_in_id,
         last_interaction_out_id,
@@ -379,6 +387,7 @@ def set_properties(
             temp_last_interaction_in_nu = np.empty(
                 temp_length, dtype=np.float64
             )
+            temp_last_interaction_in_r = np.empty(temp_length, dtype=np.float64)
             temp_last_interaction_type = np.empty(temp_length, dtype=np.int64)
             temp_last_interaction_in_id = np.empty(temp_length, dtype=np.int64)
             temp_last_interaction_out_id = np.empty(temp_length, dtype=np.int64)
@@ -390,6 +399,9 @@ def set_properties(
             temp_last_interaction_in_nu[
                 : self.length
             ] = self.last_interaction_in_nu
+            temp_last_interaction_in_r[
+                : self.length
+            ] = self.last_interaction_in_r
             temp_last_interaction_type[
                 : self.length
             ] = self.last_interaction_type
@@ -405,7 +417,8 @@ def set_properties(
             self.initial_mus = temp_initial_mus
             self.initial_rs = temp_initial_rs
             self.last_interaction_in_nu = temp_last_interaction_in_nu
-            self.last_interaction_type = temp_last_interaction_type
+            self.last_interaction_in_r = temp_last_interaction_in_r
+            self.last_interaction_type = temp_last_interaction_type 
             self.last_interaction_in_id = temp_last_interaction_in_id
             self.last_interaction_out_id = temp_last_interaction_out_id
             self.length = temp_length
@@ -415,6 +428,7 @@ def set_properties(
         self.initial_mus[self.idx] = initial_mu
         self.initial_rs[self.idx] = initial_r
         self.last_interaction_in_nu[self.idx] = last_interaction_in_nu
+        self.last_interaction_in_r[self.idx] = last_interaction_in_r
         self.last_interaction_type[self.idx] = last_interaction_type
         self.last_interaction_in_id[self.idx] = last_interaction_in_id
         self.last_interaction_out_id[self.idx] = last_interaction_out_id

tardis/montecarlo/montecarlo_numba/r_packet.py

@@ -40,6 +40,7 @@ class PacketStatus(IntEnum):
     ("index", int64),
     ("last_interaction_type", int64),
     ("last_interaction_in_nu", float64),
+    ("last_interaction_in_r", float64),
     ("last_line_interaction_in_id", int64),
     ("last_line_interaction_out_id", int64),
 ]
@@ -58,6 +59,7 @@ def __init__(self, r, mu, nu, energy, seed, index=0):
         self.index = index
         self.last_interaction_type = -1
         self.last_interaction_in_nu = 0.0
+        self.last_interaction_in_r = 0.0
         self.last_line_interaction_in_id = -1
         self.last_line_interaction_out_id = -1
 

tardis/montecarlo/montecarlo_numba/vpacket.py

@@ -282,6 +282,7 @@ def trace_vpacket_volley(
             v_packet_mu,
             r_packet.r,
             r_packet.last_interaction_in_nu,
+            r_packet.last_interaction_in_r,
             r_packet.last_interaction_type,
             r_packet.last_line_interaction_in_id,
             r_packet.last_line_interaction_out_id,

tardis/visualization/__init__.py

@@ -9,3 +9,4 @@
 from tardis.visualization.widgets.line_info import LineInfoWidget
 from tardis.visualization.widgets.custom_abundance import CustomAbundanceWidget
 from tardis.visualization.tools.sdec_plot import SDECPlotter
+from tardis.visualization.tools.interaction_radius_plot import InteractionRadiusPlotter