From 21602dceeb0db9e38b29062171d28a2aae7d0894 Mon Sep 17 00:00:00 2001
From: rodot- <jto33@drexel.edu>
Date: Wed, 4 Aug 2021 09:28:36 -0400
Subject: [PATCH 01/17] Reworked the continuum process selection function a
 little bit to work with the new continuum jitclass

---
 .../montecarlo_numba/interaction.py           | 40 +++++++++----------
 1 file changed, 20 insertions(+), 20 deletions(-)

diff --git a/tardis/montecarlo/montecarlo_numba/interaction.py b/tardis/montecarlo/montecarlo_numba/interaction.py
index 40288fb1751..b11eb5169d3 100644
--- a/tardis/montecarlo/montecarlo_numba/interaction.py
+++ b/tardis/montecarlo/montecarlo_numba/interaction.py
@@ -34,7 +34,7 @@ def scatter(r_packet, time_explosion):
  
     return comov_nu, inverse_new_doppler_factor
 
-def continuum_event(r_packet, time_explosion, plasma, chi_continuum_calculator):
+def continuum_event(r_packet, time_explosion, continuum):
 
     comov_nu, inverse_new_doppler_factor = scatter(r_packet, time_explosion)
 
@@ -43,25 +43,25 @@ def continuum_event(r_packet, time_explosion, plasma, chi_continuum_calculator):
 
     zrand = np.random.random()
 
-    (
-        chi_bf,
-        chi_bf_contributions,
-        current_continua,
-        chi_ff,
-    ) = chi_continuum_calculator(comov_nu, r_packet.current_shell_id)
-
-    cur_electron_density = numba_plasma.electron_density[
-        r_packet.current_shell_id
-    ]
-
-    chi_e = cur_electron_density * SIGMA_THOMSON
-
-    chi_nu = chi_bf + chi_ff + chi_e
-
-
-    if zrand < SIGMA_THOMSON / chi_nu:
-        thomson_scatter(r_packet, time_explosion)
-    elif zrand < (SIGMA_THOMSON + chi_bf) / chi_nu:
+    # Does this need to be re-calculated?
+    continuum.calculate(comov_nu, r_packet.current_shell_id)
+    chi_continuum = continuum.chi_bf + continuum.chi_ff
+
+    # Since trace_packet differentiates between thomson scattering
+    # and other continuum processes, we need to renormalize
+    # our odds of selecting each continuum process given that
+    # we are not thomson scattering
+    # P(bf|~e_scat) = P(~e_scat|bf) P(bf) / P(~e_scat)
+    # P(~e_scat) = 1 - P(e_scat) = 1 - chi_e / chi_nu
+    # P(bf) = chi_bf / chi_nu
+    # P(bf|~e_scat) = chi_bf / chi_nu / (1 - chi_e / chi_nu)
+    # P(bf|~e_scat) = chi_bf / (chi_nu - chi_e)
+    # Since chi_nu is the sum of ff, bf, and e_scat opacities
+    # we can just set chi_continuum = chi_nu - chi_e = chi_bf + chi_ff
+    # alternatively, this could again just be selected in
+    # trace packet, it seems pretty straightforward to do
+
+    if zrand < chi_bf / chi_continuum:
         bound_free_absorption(r_packet, time_explosion, plasma)
     else:
         free_free_absorption(r_packet, time_explosion)

From 954deb7e79cc173378bcb385e015aa375f362b44 Mon Sep 17 00:00:00 2001
From: rodot- <jto33@drexel.edu>
Date: Wed, 4 Aug 2021 10:15:40 -0400
Subject: [PATCH 02/17] started addition of free-free emission handler

---
 .../montecarlo_numba/interaction.py           | 25 +++++++++++++++++++
 1 file changed, 25 insertions(+)

diff --git a/tardis/montecarlo/montecarlo_numba/interaction.py b/tardis/montecarlo/montecarlo_numba/interaction.py
index b11eb5169d3..c67d6eff0d7 100644
--- a/tardis/montecarlo/montecarlo_numba/interaction.py
+++ b/tardis/montecarlo/montecarlo_numba/interaction.py
@@ -18,6 +18,7 @@
 from tardis.montecarlo.montecarlo_numba.utils import get_random_mu
 from tardis.montecarlo.montecarlo_numba.macro_atom import macro_atom
 
+
 def scatter(r_packet, time_explosion):
 
     old_doppler_factor = get_doppler_factor(
@@ -67,6 +68,30 @@ def continuum_event(r_packet, time_explosion, continuum):
         free_free_absorption(r_packet, time_explosion)
 
 
+def free_free_emission(r_packet, time_explosion, numba_plasma):
+    
+    inverse_doppler_factor = get_inverse_doppler_factor(
+        r_packet.r, r_packet.mu, time_explosion
+    )
+    comov_nu = numba_plasma.nu_ff_sampler(r_packet.current_shell_id)
+    r_packet.nu = comov_nu * inverse_doppler_factor
+
+    current_line_id = len(
+            numba_plasma.line_list_nu
+            ) - np.searchsorted(
+                    numba_plasma.line_list_nu[::-1], 
+                    r_packet.nu
+                    )
+
+    r_packet.next_line_id = current_line_id
+    if montecarlo_configuration.full_relativity:
+        r_packet.mu = angle_aberration_CMF_to_LF(
+            r_packet, time_explosion, r_packet.mu
+        )
+
+
+
+
 @njit(**njit_dict_no_parallel)
 def thomson_scatter(r_packet, time_explosion):
     """

From 2eae9de130eeb74e956999953fa3369dc34cdf01 Mon Sep 17 00:00:00 2001
From: rodot- <jto33@drexel.edu>
Date: Mon, 9 Aug 2021 11:51:15 -0400
Subject: [PATCH 03/17] updated interactions to work with new macroatom

---
 .../montecarlo_numba/interaction.py           | 114 ++++++++++++++----
 1 file changed, 92 insertions(+), 22 deletions(-)

diff --git a/tardis/montecarlo/montecarlo_numba/interaction.py b/tardis/montecarlo/montecarlo_numba/interaction.py
index c67d6eff0d7..0bb6633fb5b 100644
--- a/tardis/montecarlo/montecarlo_numba/interaction.py
+++ b/tardis/montecarlo/montecarlo_numba/interaction.py
@@ -16,7 +16,10 @@
     InteractionType,
 )
 from tardis.montecarlo.montecarlo_numba.utils import get_random_mu
-from tardis.montecarlo.montecarlo_numba.macro_atom import macro_atom
+from tardis.montecarlo.montecarlo_numba.macro_atom import (
+        macro_atom, MacroAtomTransitionType
+)
+
 
 
 def scatter(r_packet, time_explosion):
@@ -35,7 +38,8 @@ def scatter(r_packet, time_explosion):
  
     return comov_nu, inverse_new_doppler_factor
 
-def continuum_event(r_packet, time_explosion, continuum):
+# Maybe make the continuum selection a method of the continuum?
+def continuum_event(r_packet, time_explosion, continuum, plasma):
 
     comov_nu, inverse_new_doppler_factor = scatter(r_packet, time_explosion)
 
@@ -46,33 +50,73 @@ def continuum_event(r_packet, time_explosion, continuum):
 
     # Does this need to be re-calculated?
     continuum.calculate(comov_nu, r_packet.current_shell_id)
-    chi_continuum = continuum.chi_bf + continuum.chi_ff
-
-    # Since trace_packet differentiates between thomson scattering
-    # and other continuum processes, we need to renormalize
-    # our odds of selecting each continuum process given that
-    # we are not thomson scattering
-    # P(bf|~e_scat) = P(~e_scat|bf) P(bf) / P(~e_scat)
-    # P(~e_scat) = 1 - P(e_scat) = 1 - chi_e / chi_nu
-    # P(bf) = chi_bf / chi_nu
-    # P(bf|~e_scat) = chi_bf / chi_nu / (1 - chi_e / chi_nu)
-    # P(bf|~e_scat) = chi_bf / (chi_nu - chi_e)
-    # Since chi_nu is the sum of ff, bf, and e_scat opacities
-    # we can just set chi_continuum = chi_nu - chi_e = chi_bf + chi_ff
-    # alternatively, this could again just be selected in
-    # trace packet, it seems pretty straightforward to do
-
-    if zrand < chi_bf / chi_continuum:
-        bound_free_absorption(r_packet, time_explosion, plasma)
-    else:
-        free_free_absorption(r_packet, time_explosion)
 
+    # Need to determine if a collisional process or not.  If not:
+
+    # This somehow determines the transition type and continuum id needed
+    # but does not sample new frequencies
+    destination_level_idx = plasma.determine_continuum_macro_activation_idx(
+            comov_nu,
+            continuum.chi_bf_tot,
+            continuum.chi_ff,
+            continuum.chi_bf_contributions,
+            continuum.current_continua
+            )
+
+    transition_id, transition_type = macro_atom(
+            destination_level_idx, 
+            r_packet.current_shell_id, 
+            numba_plasma
+            )
+
+    # Then use this to get a transition id from the macroatom 
+    if transition_type == MacroAtomTransitionType.FF_EMISSION: 
+        free_free_emission(r_packet, time_explosion, numba_plasma)
+    elif transition_type = MacroAtomTransitionType.BF_EMISSION:
+        bound_free_emission(r_packet, 
+                time_explosion, 
+                numba_plasma, 
+                continuum, 
+                transition_id)
+
+    macro_atom(destination_level_id, current_shell_id, numba_plasma)
+
+
+# TODO: numbafy | Add cooling rates to numba plasma
+def get_emission_probabilities(plasma, shell):
+    
+    C_fb = sim.plasma.cool_rate_fb_tot.iloc[0, shell]
+    C_ff = sim.plasma.cool_rate_ff.iloc[0, shell]
+    C_cl = sim.plasma.cool_rate_adiabatic.iloc[0, shell]
+    C_cl += sim.plasma.cool_rate_coll_ion.iloc[0, shell]
+    
+    pi_fb = C_fb / (C_fb + C_ff + C_cl)
+    pi_ff = C_ff / (C_fb + C_ff + C_cl)
+    
+    return pi_fb, pi_ff
+
+# TODO: numbafy
+def free_free_absorption(r_packet, numba_plasma):
+    # Do the k_packet thing, but don't actually make one
+    
+    pi_fb, pi_ff = get_emission_probabilities(numba_plasma, shell)
+    
+    # a little hack
+    z = np.random.random() - pi_fb
+    if z < 0:
+        free_bound_emission(r_packet, time_explosion, numba_plasma) 
+    elif z < pi_ff:
+        free_free_emission(r_packet, time_explosion, numba_plasma)
+    else:
+        macroatom()
 
 def free_free_emission(r_packet, time_explosion, numba_plasma):
     
     inverse_doppler_factor = get_inverse_doppler_factor(
         r_packet.r, r_packet.mu, time_explosion
     )
+    # Need to get the sampler into numba somehow
+    # maybe I'll update the numba_plasma?
     comov_nu = numba_plasma.nu_ff_sampler(r_packet.current_shell_id)
     r_packet.nu = comov_nu * inverse_doppler_factor
 
@@ -89,7 +133,33 @@ def free_free_emission(r_packet, time_explosion, numba_plasma):
             r_packet, time_explosion, r_packet.mu
         )
 
+def get_current_continuum_id(nu, continuum, plasma):
+
+    # TODO:
+    continuum_id = 0
+    return continuum_id
 
+def free_bound_emission(r_packet, time_explosion, numba_plasma, continuum, continuum_d):
+
+    inverse_doppler_factor = get_inverse_doppler_factor(
+        r_packet.r, r_packet.mu, time_explosion
+    )
+
+    comov_nu = numba_plasma.nu_fb_sampler(
+            r_packet.current_shell_id, 
+            continuum_id
+            )
+    
+    r_packet.nu = comov_nu * inverse_doppler_factor
+
+    current_line_id = len(
+            numba_plasma.line_list_nu
+            ) - np.searchsorted(
+                    numba_plasma.line_list_nu[::-1], 
+                    r_packet.nu
+                    )
+
+    r_packet.next_line_id = current_line_id
 
 
 @njit(**njit_dict_no_parallel)

From 8762f0b3a2820243554ef64431be029c1f7a0172 Mon Sep 17 00:00:00 2001
From: rodot- <jto33@drexel.edu>
Date: Mon, 9 Aug 2021 15:17:21 -0400
Subject: [PATCH 04/17] Added continuum process frquency samplers to the
 continuum jitclass constructor

---
 tardis/montecarlo/montecarlo_numba/base.py        |  5 ++++-
 .../montecarlo_numba/numba_interface.py           | 15 ++++++++++++++-
 2 files changed, 18 insertions(+), 2 deletions(-)

diff --git a/tardis/montecarlo/montecarlo_numba/base.py b/tardis/montecarlo/montecarlo_numba/base.py
index 885d7805a14..1cf3122890d 100644
--- a/tardis/montecarlo/montecarlo_numba/base.py
+++ b/tardis/montecarlo/montecarlo_numba/base.py
@@ -58,7 +58,10 @@ def montecarlo_radial1d(model, plasma, runner):
     packet_seeds = montecarlo_configuration.packet_seeds
 
     number_of_vpackets = montecarlo_configuration.number_of_vpackets
-    ContinuumObject = create_continuum_class(plasma.chi_continuum_calculator)
+    ContinuumObject = create_continuum_class(
+            plasma.chi_continuum_calculator,
+            plasma.nu_fb_sampler,
+            plasma.nu_ff_sampler)
 
     (
         v_packets_energy_hist,
diff --git a/tardis/montecarlo/montecarlo_numba/numba_interface.py b/tardis/montecarlo/montecarlo_numba/numba_interface.py
index 88b8406974b..8acb70ff0ef 100644
--- a/tardis/montecarlo/montecarlo_numba/numba_interface.py
+++ b/tardis/montecarlo/montecarlo_numba/numba_interface.py
@@ -318,7 +318,10 @@ def set_properties(
 
 
 
-def create_continuum_class(chi_continuum_calculator):
+def create_continuum_class(
+        chi_continuum_calculator, 
+        nu_fb_sampler, 
+        nu_ff_sampler):
     """called before mainloop"""
     continuum_spec = [
             ("chi_bf_tot", float64),
@@ -348,6 +351,16 @@ def calculate(self, nu, shell):
             self.chi_ff,
             ) = chi_continuum_calculator(nu, shell)
 
+        def sample_nu_free_bound(self, shell, continuum_id):
+
+            return nu_fb_sampler(shell, continuum_id)
+
+        def sample_nu_free_free(self, shell):
+
+            return nu_ff_sampler(shell)
+
+
+
     @njit(**njit_dict_no_parallel)
     def continuum_constructor():
         return Continuum()

From 5f42d66f2b56b13ca15b80fa0fddf9a7b6d6588c Mon Sep 17 00:00:00 2001
From: rodot- <jto33@drexel.edu>
Date: Wed, 11 Aug 2021 08:30:31 -0400
Subject: [PATCH 05/17] Fixed typo from merge conflict

---
 tardis/montecarlo/montecarlo_numba/base.py | 4 ----
 1 file changed, 4 deletions(-)

diff --git a/tardis/montecarlo/montecarlo_numba/base.py b/tardis/montecarlo/montecarlo_numba/base.py
index da863b76e4b..52215e49835 100644
--- a/tardis/montecarlo/montecarlo_numba/base.py
+++ b/tardis/montecarlo/montecarlo_numba/base.py
@@ -208,15 +208,11 @@ def montecarlo_main_loop(
         vpacket_collection = vpacket_collections[i]
 
         loop = single_packet_loop(
-
             r_packet,
             numba_model,
             numba_plasma,
             estimators,
             vpacket_collection,
-            continuum,
-
-            r_packet, numba_model, numba_plasma, estimators, vpacket_collection,
             continuum
         )
         # if loop and 'stop' in loop:

From bd1483d5d45e9b1b0933a02e3327636c4159981d Mon Sep 17 00:00:00 2001
From: rodot- <jto33@drexel.edu>
Date: Wed, 11 Aug 2021 08:32:35 -0400
Subject: [PATCH 06/17] Fixed spacing

---
 tardis/plasma/properties/continuum_processes.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/tardis/plasma/properties/continuum_processes.py b/tardis/plasma/properties/continuum_processes.py
index 51daf9bc44c..584c169e212 100644
--- a/tardis/plasma/properties/continuum_processes.py
+++ b/tardis/plasma/properties/continuum_processes.py
@@ -1029,7 +1029,7 @@ def nu_fb(shell, continuum_id):
             zrand = np.random.random()
             idx = np.searchsorted(em[start:end], zrand, side="right")
 
-            return phot_nus[idx] - (em[idx] - zrand) / (em[idx] - em[idx-1]) * (phot_nus[idx] - phot_nus[idx-1])
+            return phot_nus[idx] - (em[idx] - zrand) / (em[idx] - em[idx - 1]) * (phot_nus[idx] - phot_nus[idx - 1])
 
         return nu_fb
 

From ecfdb5bf1abbd714b6e9c32dddc80e765c6c7bfa Mon Sep 17 00:00:00 2001
From: rodot- <jto33@drexel.edu>
Date: Wed, 11 Aug 2021 08:37:56 -0400
Subject: [PATCH 07/17] small formatting updates to get the changes to track

---
 tardis/montecarlo/montecarlo_numba/r_packet.py   | 1 +
 tardis/plasma/properties/property_collections.py | 2 +-
 2 files changed, 2 insertions(+), 1 deletion(-)

diff --git a/tardis/montecarlo/montecarlo_numba/r_packet.py b/tardis/montecarlo/montecarlo_numba/r_packet.py
index 8ff94d19c0c..ba466be3a0b 100644
--- a/tardis/montecarlo/montecarlo_numba/r_packet.py
+++ b/tardis/montecarlo/montecarlo_numba/r_packet.py
@@ -210,6 +210,7 @@ def trace_packet(
         distance = min(distance_trace, distance_boundary, distance_continuum)
 
         if distance_trace != 0:
+
             if distance == distance_boundary:
                 interaction_type = InteractionType.BOUNDARY  # BOUNDARY
                 r_packet.next_line_id = cur_line_id
diff --git a/tardis/plasma/properties/property_collections.py b/tardis/plasma/properties/property_collections.py
index dba6bbcebd7..c6d281b69fd 100644
--- a/tardis/plasma/properties/property_collections.py
+++ b/tardis/plasma/properties/property_collections.py
@@ -93,7 +93,7 @@ class PlasmaPropertyCollection(list):
 )
 continuum_interaction_properties = PlasmaPropertyCollection(
     [
-        MacroAtomData,  # TODO: Remove this once everything works
+        MacroAtomData, # TODO: Remove this once everything works
         PhotoIonizationData,
         SpontRecombRateCoeff,
         PhotoIonRateCoeff,

From b979fb78f5812b717a19db0794c9a5c3c195f1e3 Mon Sep 17 00:00:00 2001
From: rodot- <jto33@drexel.edu>
Date: Wed, 11 Aug 2021 08:38:51 -0400
Subject: [PATCH 08/17] propagated continuum through single packet loop

---
 tardis/montecarlo/montecarlo_numba/single_packet_loop.py | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/tardis/montecarlo/montecarlo_numba/single_packet_loop.py b/tardis/montecarlo/montecarlo_numba/single_packet_loop.py
index f308a1d91fe..a34f70e58e6 100644
--- a/tardis/montecarlo/montecarlo_numba/single_packet_loop.py
+++ b/tardis/montecarlo/montecarlo_numba/single_packet_loop.py
@@ -38,7 +38,7 @@
 @njit
 def single_packet_loop(
     r_packet, numba_model, numba_plasma, estimators, vpacket_collection,
-    chi_bf_interpolator
+    continuum
 ):
     """
     Parameters
@@ -78,7 +78,7 @@ def single_packet_loop(
     while r_packet.status == PacketStatus.IN_PROCESS:
         distance, interaction_type, delta_shell = trace_packet(
             r_packet, numba_model, numba_plasma,
-            estimators, chi_bf_interpolator
+            estimators, continuum
         )
 
         if interaction_type == InteractionType.BOUNDARY:

From bea3c39b23572ba97e4da99cc49c6fba8d3b763d Mon Sep 17 00:00:00 2001
From: rodot- <jto33@drexel.edu>
Date: Wed, 11 Aug 2021 08:49:09 -0400
Subject: [PATCH 09/17] Cleaned up file, bound-free and free-free emission
 pathways should be in the correct place now

---
 .../montecarlo_numba/interaction.py           | 49 +++++++++----------
 1 file changed, 24 insertions(+), 25 deletions(-)

diff --git a/tardis/montecarlo/montecarlo_numba/interaction.py b/tardis/montecarlo/montecarlo_numba/interaction.py
index 0bb6633fb5b..5775ddfa8c7 100644
--- a/tardis/montecarlo/montecarlo_numba/interaction.py
+++ b/tardis/montecarlo/montecarlo_numba/interaction.py
@@ -39,7 +39,8 @@ def scatter(r_packet, time_explosion):
     return comov_nu, inverse_new_doppler_factor
 
 # Maybe make the continuum selection a method of the continuum?
-def continuum_event(r_packet, time_explosion, continuum, plasma):
+@njit(**njit_dict_no_parallel)
+def continuum_event(r_packet, time_explosion, continuum, numba_plasma):
 
     comov_nu, inverse_new_doppler_factor = scatter(r_packet, time_explosion)
 
@@ -79,9 +80,6 @@ def continuum_event(r_packet, time_explosion, continuum, plasma):
                 continuum, 
                 transition_id)
 
-    macro_atom(destination_level_id, current_shell_id, numba_plasma)
-
-
 # TODO: numbafy | Add cooling rates to numba plasma
 def get_emission_probabilities(plasma, shell):
     
@@ -110,6 +108,19 @@ def free_free_absorption(r_packet, numba_plasma):
     else:
         macroatom()
 
+@njit(**njit_dict_no_parallel)
+def get_current_line_id(nu, numba_plasma):
+    '''
+    Get the next line id corresponding to a packet at frequency nu
+    '''
+
+    reverse_line_list = numba_plasma.line_list_nu[::-1]
+    number_of_lines = len(numba_plasma.line_list_nu)
+    line_id = number_of_lines - np.searchsorted(reverse_line_list, nu)
+    return line_id
+
+
+@njit(**njit_dict_no_parallel)
 def free_free_emission(r_packet, time_explosion, numba_plasma):
     
     inverse_doppler_factor = get_inverse_doppler_factor(
@@ -119,26 +130,15 @@ def free_free_emission(r_packet, time_explosion, numba_plasma):
     # maybe I'll update the numba_plasma?
     comov_nu = numba_plasma.nu_ff_sampler(r_packet.current_shell_id)
     r_packet.nu = comov_nu * inverse_doppler_factor
-
-    current_line_id = len(
-            numba_plasma.line_list_nu
-            ) - np.searchsorted(
-                    numba_plasma.line_list_nu[::-1], 
-                    r_packet.nu
-                    )
-
+    current_line_id = get_current_line_id(r_packet.nu, numba_plasma) 
     r_packet.next_line_id = current_line_id
+    
     if montecarlo_configuration.full_relativity:
         r_packet.mu = angle_aberration_CMF_to_LF(
             r_packet, time_explosion, r_packet.mu
         )
 
-def get_current_continuum_id(nu, continuum, plasma):
-
-    # TODO:
-    continuum_id = 0
-    return continuum_id
-
+@njit(**njit_dict_no_parallel)
 def free_bound_emission(r_packet, time_explosion, numba_plasma, continuum, continuum_d):
 
     inverse_doppler_factor = get_inverse_doppler_factor(
@@ -151,15 +151,14 @@ def free_bound_emission(r_packet, time_explosion, numba_plasma, continuum, conti
             )
     
     r_packet.nu = comov_nu * inverse_doppler_factor
+    current_line_id = get_current_line_id(r_packet.nu, numba_plasma)
+    r_packet.next_line_id = current_line_id
 
-    current_line_id = len(
-            numba_plasma.line_list_nu
-            ) - np.searchsorted(
-                    numba_plasma.line_list_nu[::-1], 
-                    r_packet.nu
-                    )
+    if montecarlo_configuration.full_relativity:
+        r_packet.mu = angle_aberration_CMF_to_LF(
+            r_packet, time_explosion, r_packet.mu
+        )
 
-    r_packet.next_line_id = current_line_id
 
 
 @njit(**njit_dict_no_parallel)

From 5d23b2c7b8966e2a79d5a5ef37e229fc9ba37fa2 Mon Sep 17 00:00:00 2001
From: rodot- <jto33@drexel.edu>
Date: Wed, 11 Aug 2021 08:59:00 -0400
Subject: [PATCH 10/17] added fixture for a new continuum object

---
 .../montecarlo/montecarlo_numba/tests/conftest.py | 15 +++++++++++++++
 1 file changed, 15 insertions(+)

diff --git a/tardis/montecarlo/montecarlo_numba/tests/conftest.py b/tardis/montecarlo/montecarlo_numba/tests/conftest.py
index 53cf34788d7..d6aa94120f8 100644
--- a/tardis/montecarlo/montecarlo_numba/tests/conftest.py
+++ b/tardis/montecarlo/montecarlo_numba/tests/conftest.py
@@ -24,6 +24,21 @@ def nb_simulation_verysimple(config_verysimple, atomic_dataset):
     return sim
 
 
+@pytest.fixture(scope="package")
+def verysimple_continuum(nu_simulation_verysimple):
+    plasma = nu_simulation_verysimple
+    numba_plasma = numba_interface.numba_plasma_initialize(plasma, "macroatom")
+
+    ContinuumClass = create_continuum_class(
+            plasma.chi_continuum_calculator,
+            plasma.nu_bf_sampler,
+            plasma.nu_ff_sampelr
+            )
+
+    continuum = ContinuumClass()
+    return continuum
+
+
 @pytest.fixture()
 def verysimple_collection(nb_simulation_verysimple):
     runner = nb_simulation_verysimple.runner

From af93d048def2f72c74d4cfecd5b02bd35cce2033 Mon Sep 17 00:00:00 2001
From: rodot- <jto33@drexel.edu>
Date: Wed, 11 Aug 2021 09:10:55 -0400
Subject: [PATCH 11/17] Added placeholder test for calculating the continuum
 opacities.  I'm weary to put in exact values until we're done

---
 .../tests/test_numba_interface.py              | 18 ++++++++++++++++++
 1 file changed, 18 insertions(+)

diff --git a/tardis/montecarlo/montecarlo_numba/tests/test_numba_interface.py b/tardis/montecarlo/montecarlo_numba/tests/test_numba_interface.py
index 75b287e0faa..7c3715ad9a0 100644
--- a/tardis/montecarlo/montecarlo_numba/tests/test_numba_interface.py
+++ b/tardis/montecarlo/montecarlo_numba/tests/test_numba_interface.py
@@ -4,6 +4,24 @@
 import numpy as np
 
 
+@pytest.mark.parametrize(
+       ["current_shell_id", "nu"],
+       [(0, 0.6), (1, 0.4)]
+       ) 
+def test_coninuum_opacities(
+    verysimple_continuum, current_shell_id, nu):
+
+    continuum = verysimple_continuum
+    continuum.calculate(nu, current_shell_id)
+    print(current_shell_id, nu)
+    print(continuum.chi_bf_tot)
+    print(continuum.chi_bf_contributions)
+    print(continuum.current_continuua)
+    print(continuum.x_sect_bfs)
+    print(continuum.chi_ff)
+
+    
+
 @pytest.mark.parametrize("input_params", ["scatter", "macroatom", "downbranch"])
 def test_numba_plasma_initialize(nb_simulation_verysimple, input_params):
     line_interaction_type = input_params

From fc5396c7397ce944b0ebd7bbe00fa322ce87d71c Mon Sep 17 00:00:00 2001
From: rodot- <jto33@drexel.edu>
Date: Wed, 11 Aug 2021 09:28:54 -0400
Subject: [PATCH 12/17] reworked the contruction of the continuum class so that
 only a plasma object needs to be passed

---
 tardis/montecarlo/montecarlo_numba/base.py    |  6 +----
 .../montecarlo_numba/numba_interface.py       | 15 ++++++++++-
 .../montecarlo_numba/tests/conftest.py        |  6 +----
 .../tests/test_numba_interface.py             | 25 ++++++++++++++++++-
 4 files changed, 40 insertions(+), 12 deletions(-)

diff --git a/tardis/montecarlo/montecarlo_numba/base.py b/tardis/montecarlo/montecarlo_numba/base.py
index 52215e49835..3f97ab125d5 100644
--- a/tardis/montecarlo/montecarlo_numba/base.py
+++ b/tardis/montecarlo/montecarlo_numba/base.py
@@ -59,11 +59,7 @@ def montecarlo_radial1d(model, plasma, runner):
     packet_seeds = montecarlo_configuration.packet_seeds
 
     number_of_vpackets = montecarlo_configuration.number_of_vpackets
-    ContinuumObject = create_continuum_class(
-            plasma.chi_continuum_calculator,
-            plasma.nu_fb_sampler,
-            plasma.nu_ff_sampler)
-
+    ContinuumObject = create_continuum_class(plasma)
     (
         v_packets_energy_hist,
         last_interaction_type,
diff --git a/tardis/montecarlo/montecarlo_numba/numba_interface.py b/tardis/montecarlo/montecarlo_numba/numba_interface.py
index 93ae9b00d23..7eabfa68cd4 100644
--- a/tardis/montecarlo/montecarlo_numba/numba_interface.py
+++ b/tardis/montecarlo/montecarlo_numba/numba_interface.py
@@ -318,11 +318,16 @@ def set_properties(
 
 
 
-def create_continuum_class(
+def create_continuum_class(plasma)
         chi_continuum_calculator, 
         nu_fb_sampler, 
         nu_ff_sampler):
     """called before mainloop"""
+
+    chi_continuum_calculator = plasma.chi_continuum_calculator
+    nu_fb_sampler = plasma.nu_fb_sampler
+    nu_ff_sampler = plasma.nu_ff_sampler
+    get_macro_activation_idx = plasma.determine_continuum_macro_activation_idx
     continuum_spec = [
             ("chi_bf_tot", float64),
             ("chi_bf_contributions", float64[:]),
@@ -359,6 +364,14 @@ def sample_nu_free_free(self, shell):
 
             return nu_ff_sampler(shell)
 
+        def determine_macro_activation_idx(self, nu, shell):
+
+            self.calculate(nu, shell) # maybe don't do this here
+            idx = get_macro_activation_idx(
+                    nu, self.chi_bf, self.chi_ff, 
+                    self.chi_bf_contributions, self.current_continua
+                    )
+            return idx
 
 
     @njit(**njit_dict_no_parallel)
diff --git a/tardis/montecarlo/montecarlo_numba/tests/conftest.py b/tardis/montecarlo/montecarlo_numba/tests/conftest.py
index d6aa94120f8..480ca530b72 100644
--- a/tardis/montecarlo/montecarlo_numba/tests/conftest.py
+++ b/tardis/montecarlo/montecarlo_numba/tests/conftest.py
@@ -29,11 +29,7 @@ def verysimple_continuum(nu_simulation_verysimple):
     plasma = nu_simulation_verysimple
     numba_plasma = numba_interface.numba_plasma_initialize(plasma, "macroatom")
 
-    ContinuumClass = create_continuum_class(
-            plasma.chi_continuum_calculator,
-            plasma.nu_bf_sampler,
-            plasma.nu_ff_sampelr
-            )
+    ContinuumClass = create_continuum_class(plasma)
 
     continuum = ContinuumClass()
     return continuum
diff --git a/tardis/montecarlo/montecarlo_numba/tests/test_numba_interface.py b/tardis/montecarlo/montecarlo_numba/tests/test_numba_interface.py
index 7c3715ad9a0..0311b677d13 100644
--- a/tardis/montecarlo/montecarlo_numba/tests/test_numba_interface.py
+++ b/tardis/montecarlo/montecarlo_numba/tests/test_numba_interface.py
@@ -20,7 +20,30 @@ def test_coninuum_opacities(
     print(continuum.x_sect_bfs)
     print(continuum.chi_ff)
 
-    
+ 
+@pytest.mark.parametrize(
+       "current_shell_id",
+       [0, 1, 2]
+       ) 
+def test_continuum_free_free_sampler(
+        verysimple_continuum, current_shell_id):
+
+    continuum = verysimple_continuum
+    nu = continuum.sample_nu_free_free(current_shell_id)
+    print(nu)
+
+@pytest.mark.parametrize(
+       ["current_shell_id", "continuum_id"],
+       [(0, 0), (1, 0)]
+       ) 
+def test_coninuum_opacities(
+    verysimple_continuum, current_shell_id, continuum_id):
+
+    continuum = verysimple_continuum
+    nu = continuum.sample_nu_free_bound(current_shell_id, continuum_id)
+    print(nu)
+
+
 
 @pytest.mark.parametrize("input_params", ["scatter", "macroatom", "downbranch"])
 def test_numba_plasma_initialize(nb_simulation_verysimple, input_params):

From 96ad962ca6b4efd639aabddceb0febffd2244722 Mon Sep 17 00:00:00 2001
From: rodot- <jto33@drexel.edu>
Date: Wed, 11 Aug 2021 09:32:24 -0400
Subject: [PATCH 13/17] Updated interaction to now get the macro activation
 level from the continuum object

---
 tardis/montecarlo/montecarlo_numba/interaction.py | 9 ++-------
 1 file changed, 2 insertions(+), 7 deletions(-)

diff --git a/tardis/montecarlo/montecarlo_numba/interaction.py b/tardis/montecarlo/montecarlo_numba/interaction.py
index 5775ddfa8c7..539f104c7c6 100644
--- a/tardis/montecarlo/montecarlo_numba/interaction.py
+++ b/tardis/montecarlo/montecarlo_numba/interaction.py
@@ -56,13 +56,8 @@ def continuum_event(r_packet, time_explosion, continuum, numba_plasma):
 
     # This somehow determines the transition type and continuum id needed
     # but does not sample new frequencies
-    destination_level_idx = plasma.determine_continuum_macro_activation_idx(
-            comov_nu,
-            continuum.chi_bf_tot,
-            continuum.chi_ff,
-            continuum.chi_bf_contributions,
-            continuum.current_continua
-            )
+    destination_level_idx = continuum.determine_macro_activation_idx(
+            comov_nu, r_packet.current_shell_id)
 
     transition_id, transition_type = macro_atom(
             destination_level_idx, 

From 86c885c85d0826baa4804964f4fcd83b03814f74 Mon Sep 17 00:00:00 2001
From: rodot- <jto33@drexel.edu>
Date: Wed, 11 Aug 2021 09:36:16 -0400
Subject: [PATCH 14/17] refer to previous commit

---
 tardis/montecarlo/montecarlo_numba/numba_interface.py | 3 ---
 1 file changed, 3 deletions(-)

diff --git a/tardis/montecarlo/montecarlo_numba/numba_interface.py b/tardis/montecarlo/montecarlo_numba/numba_interface.py
index 7eabfa68cd4..7cd3d86d159 100644
--- a/tardis/montecarlo/montecarlo_numba/numba_interface.py
+++ b/tardis/montecarlo/montecarlo_numba/numba_interface.py
@@ -319,9 +319,6 @@ def set_properties(
 
 
 def create_continuum_class(plasma)
-        chi_continuum_calculator, 
-        nu_fb_sampler, 
-        nu_ff_sampler):
     """called before mainloop"""
 
     chi_continuum_calculator = plasma.chi_continuum_calculator

From aa411459c14b6175b70a6a4602d146ad8c074c19 Mon Sep 17 00:00:00 2001
From: rodot- <jto33@drexel.edu>
Date: Wed, 11 Aug 2021 09:40:43 -0400
Subject: [PATCH 15/17] Added the continuum_process into the single packet loop
 so packets can now actually interact with the contniuum

---
 tardis/montecarlo/montecarlo_numba/interaction.py  |  3 ++-
 .../montecarlo_numba/single_packet_loop.py         | 14 ++++++++++++++
 2 files changed, 16 insertions(+), 1 deletion(-)

diff --git a/tardis/montecarlo/montecarlo_numba/interaction.py b/tardis/montecarlo/montecarlo_numba/interaction.py
index 539f104c7c6..59676f72931 100644
--- a/tardis/montecarlo/montecarlo_numba/interaction.py
+++ b/tardis/montecarlo/montecarlo_numba/interaction.py
@@ -50,12 +50,13 @@ def continuum_event(r_packet, time_explosion, continuum, numba_plasma):
     zrand = np.random.random()
 
     # Does this need to be re-calculated?
-    continuum.calculate(comov_nu, r_packet.current_shell_id)
+    # continuum.calculate(comov_nu, r_packet.current_shell_id)
 
     # Need to determine if a collisional process or not.  If not:
 
     # This somehow determines the transition type and continuum id needed
     # but does not sample new frequencies
+    # This reruns continuum.calculate
     destination_level_idx = continuum.determine_macro_activation_idx(
             comov_nu, r_packet.current_shell_id)
 
diff --git a/tardis/montecarlo/montecarlo_numba/single_packet_loop.py b/tardis/montecarlo/montecarlo_numba/single_packet_loop.py
index a34f70e58e6..1939bdba993 100644
--- a/tardis/montecarlo/montecarlo_numba/single_packet_loop.py
+++ b/tardis/montecarlo/montecarlo_numba/single_packet_loop.py
@@ -17,6 +17,7 @@
     InteractionType,
     thomson_scatter,
     line_scatter,
+    continuum_event,
 )
 from tardis.montecarlo.montecarlo_numba.numba_interface import (
     LineInteractionType,
@@ -116,6 +117,19 @@ def single_packet_loop(
             trace_vpacket_volley(
                 r_packet, vpacket_collection, numba_model, numba_plasma
             )
+        elif interaction_type = InteractionType.CONTINUUM_PROCESS:
+            r_packet.last_interaction_type = InteractionType.CONTINUUM_PROCESS
+             move_r_packet(
+                r_packet, distance, numba_model.time_explosion, estimators
+            )
+            continuum_event(r_packet, numba_model.time_explosion,
+                    continuum, numba_plasma)
+
+            trace_vpacket_volley(
+                r_packet, vpacket_collection, numba_model, numba_plasma
+            )
+
+
         if mc_tracker.DEBUG_MODE:
             r_packet_track_nu.append(r_packet.nu)
             r_packet_track_mu.append(r_packet.mu)

From 3d4cd82b8d5bf8d165586a5a15f3d3592cc08345 Mon Sep 17 00:00:00 2001
From: rodot- <jto33@drexel.edu>
Date: Wed, 11 Aug 2021 11:33:33 -0400
Subject: [PATCH 16/17] various small bugfixes to get the numba funcs to
 compile properly

---
 .../montecarlo_numba/interaction.py           | 19 ++++++++++---------
 .../montecarlo_numba/numba_interface.py       |  4 ++--
 .../montecarlo/montecarlo_numba/r_packet.py   | 11 -----------
 .../montecarlo_numba/single_packet_loop.py    |  4 ++--
 4 files changed, 14 insertions(+), 24 deletions(-)

diff --git a/tardis/montecarlo/montecarlo_numba/interaction.py b/tardis/montecarlo/montecarlo_numba/interaction.py
index 59676f72931..2adc5db42d4 100644
--- a/tardis/montecarlo/montecarlo_numba/interaction.py
+++ b/tardis/montecarlo/montecarlo_numba/interaction.py
@@ -1,4 +1,5 @@
 from numba import njit
+import numpy as np
 from tardis.montecarlo.montecarlo_numba import njit_dict, njit_dict_no_parallel
 from tardis.montecarlo.montecarlo_numba.numba_interface import (
     LineInteractionType,
@@ -22,6 +23,7 @@
 
 
 
+@njit(**njit_dict_no_parallel)
 def scatter(r_packet, time_explosion):
 
     old_doppler_factor = get_doppler_factor(
@@ -68,8 +70,8 @@ def continuum_event(r_packet, time_explosion, continuum, numba_plasma):
 
     # Then use this to get a transition id from the macroatom 
     if transition_type == MacroAtomTransitionType.FF_EMISSION: 
-        free_free_emission(r_packet, time_explosion, numba_plasma)
-    elif transition_type = MacroAtomTransitionType.BF_EMISSION:
+        free_free_emission(r_packet, time_explosion, numba_plasma, continuum)
+    elif transition_type == MacroAtomTransitionType.BF_EMISSION:
         bound_free_emission(r_packet, 
                 time_explosion, 
                 numba_plasma, 
@@ -98,7 +100,7 @@ def free_free_absorption(r_packet, numba_plasma):
     # a little hack
     z = np.random.random() - pi_fb
     if z < 0:
-        free_bound_emission(r_packet, time_explosion, numba_plasma) 
+        bound_free_emission(r_packet, time_explosion, numba_plasma) 
     elif z < pi_ff:
         free_free_emission(r_packet, time_explosion, numba_plasma)
     else:
@@ -117,14 +119,14 @@ def get_current_line_id(nu, numba_plasma):
 
 
 @njit(**njit_dict_no_parallel)
-def free_free_emission(r_packet, time_explosion, numba_plasma):
+def free_free_emission(r_packet, time_explosion, numba_plasma, continuum):
     
     inverse_doppler_factor = get_inverse_doppler_factor(
         r_packet.r, r_packet.mu, time_explosion
     )
     # Need to get the sampler into numba somehow
     # maybe I'll update the numba_plasma?
-    comov_nu = numba_plasma.nu_ff_sampler(r_packet.current_shell_id)
+    comov_nu = continuum.sample_nu_free_free(r_packet.current_shell_id)
     r_packet.nu = comov_nu * inverse_doppler_factor
     current_line_id = get_current_line_id(r_packet.nu, numba_plasma) 
     r_packet.next_line_id = current_line_id
@@ -135,16 +137,15 @@ def free_free_emission(r_packet, time_explosion, numba_plasma):
         )
 
 @njit(**njit_dict_no_parallel)
-def free_bound_emission(r_packet, time_explosion, numba_plasma, continuum, continuum_d):
+def bound_free_emission(r_packet, time_explosion, numba_plasma, continuum, continuum_id):
 
     inverse_doppler_factor = get_inverse_doppler_factor(
         r_packet.r, r_packet.mu, time_explosion
     )
 
-    comov_nu = numba_plasma.nu_fb_sampler(
+    comov_nu = continuum.sample_nu_free_bound(
             r_packet.current_shell_id, 
-            continuum_id
-            )
+            continuum_id)
     
     r_packet.nu = comov_nu * inverse_doppler_factor
     current_line_id = get_current_line_id(r_packet.nu, numba_plasma)
diff --git a/tardis/montecarlo/montecarlo_numba/numba_interface.py b/tardis/montecarlo/montecarlo_numba/numba_interface.py
index 7cd3d86d159..627d1fee660 100644
--- a/tardis/montecarlo/montecarlo_numba/numba_interface.py
+++ b/tardis/montecarlo/montecarlo_numba/numba_interface.py
@@ -318,7 +318,7 @@ def set_properties(
 
 
 
-def create_continuum_class(plasma)
+def create_continuum_class(plasma):
     """called before mainloop"""
 
     chi_continuum_calculator = plasma.chi_continuum_calculator
@@ -365,7 +365,7 @@ def determine_macro_activation_idx(self, nu, shell):
 
             self.calculate(nu, shell) # maybe don't do this here
             idx = get_macro_activation_idx(
-                    nu, self.chi_bf, self.chi_ff, 
+                    nu, self.chi_bf_tot, self.chi_ff, 
                     self.chi_bf_contributions, self.current_continua
                     )
             return idx
diff --git a/tardis/montecarlo/montecarlo_numba/r_packet.py b/tardis/montecarlo/montecarlo_numba/r_packet.py
index ba466be3a0b..eec2f19b341 100644
--- a/tardis/montecarlo/montecarlo_numba/r_packet.py
+++ b/tardis/montecarlo/montecarlo_numba/r_packet.py
@@ -161,17 +161,6 @@ def trace_packet(
     chi_continuum = chi_e + chi_bf + chi_ff
     distance_continuum = tau_event / chi_continuum
 
-    (
-        chi_bf,
-        chi_bf_contributions,
-        current_continua,
-        x_sect_bfs,
-        chi_ff,
-    ) = chi_continuum_calculator(comov_nu, r_packet.current_shell_id)
-
-    chi_continuum = chi_e + chi_bf + chi_ff
-    distance_continuum = tau_event / chi_continuum
-
     cur_line_id = start_line_id  # initializing varibale for Numba
     # - do not remove
     last_line_id = len(numba_plasma.line_list_nu) - 1
diff --git a/tardis/montecarlo/montecarlo_numba/single_packet_loop.py b/tardis/montecarlo/montecarlo_numba/single_packet_loop.py
index 1939bdba993..a269fab899b 100644
--- a/tardis/montecarlo/montecarlo_numba/single_packet_loop.py
+++ b/tardis/montecarlo/montecarlo_numba/single_packet_loop.py
@@ -117,9 +117,9 @@ def single_packet_loop(
             trace_vpacket_volley(
                 r_packet, vpacket_collection, numba_model, numba_plasma
             )
-        elif interaction_type = InteractionType.CONTINUUM_PROCESS:
+        elif interaction_type == InteractionType.CONTINUUM_PROCESS:
             r_packet.last_interaction_type = InteractionType.CONTINUUM_PROCESS
-             move_r_packet(
+            move_r_packet(
                 r_packet, distance, numba_model.time_explosion, estimators
             )
             continuum_event(r_packet, numba_model.time_explosion,

From 166c0e6d5a9bee944d5dd1a1cf8ffe2d86d2702c Mon Sep 17 00:00:00 2001
From: rodot- <jto33@drexel.edu>
Date: Wed, 11 Aug 2021 11:36:47 -0400
Subject: [PATCH 17/17] Removed extra functions

---
 .../montecarlo_numba/interaction.py           | 28 -------------------
 1 file changed, 28 deletions(-)

diff --git a/tardis/montecarlo/montecarlo_numba/interaction.py b/tardis/montecarlo/montecarlo_numba/interaction.py
index 2adc5db42d4..8e7c5ffdf9a 100644
--- a/tardis/montecarlo/montecarlo_numba/interaction.py
+++ b/tardis/montecarlo/montecarlo_numba/interaction.py
@@ -78,34 +78,6 @@ def continuum_event(r_packet, time_explosion, continuum, numba_plasma):
                 continuum, 
                 transition_id)
 
-# TODO: numbafy | Add cooling rates to numba plasma
-def get_emission_probabilities(plasma, shell):
-    
-    C_fb = sim.plasma.cool_rate_fb_tot.iloc[0, shell]
-    C_ff = sim.plasma.cool_rate_ff.iloc[0, shell]
-    C_cl = sim.plasma.cool_rate_adiabatic.iloc[0, shell]
-    C_cl += sim.plasma.cool_rate_coll_ion.iloc[0, shell]
-    
-    pi_fb = C_fb / (C_fb + C_ff + C_cl)
-    pi_ff = C_ff / (C_fb + C_ff + C_cl)
-    
-    return pi_fb, pi_ff
-
-# TODO: numbafy
-def free_free_absorption(r_packet, numba_plasma):
-    # Do the k_packet thing, but don't actually make one
-    
-    pi_fb, pi_ff = get_emission_probabilities(numba_plasma, shell)
-    
-    # a little hack
-    z = np.random.random() - pi_fb
-    if z < 0:
-        bound_free_emission(r_packet, time_explosion, numba_plasma) 
-    elif z < pi_ff:
-        free_free_emission(r_packet, time_explosion, numba_plasma)
-    else:
-        macroatom()
-
 @njit(**njit_dict_no_parallel)
 def get_current_line_id(nu, numba_plasma):
     '''