Fix formal integral spikes bug (see tardis-sn#970) (tardis-sn#971)

tardis/montecarlo/formal_integral.py

@@ -199,9 +199,11 @@ def interpolate_integrator_quantities(self, att_S_ul, Jredlu,
         runner.electron_densities_integ = interp1d(
                 r_middle, plasma.electron_densities,
                 fill_value='extrapolate', kind='nearest')(r_middle_integ)
+        # Assume tau_sobolevs to be constant within a shell
+        # (as in the MC simulation)
         runner.tau_sobolevs_integ = interp1d(
                 r_middle, plasma.tau_sobolevs,
-                fill_value='extrapolate')(r_middle_integ)
+                fill_value='extrapolate', kind='nearest')(r_middle_integ)
         att_S_ul = interp1d(
                 r_middle, att_S_ul, fill_value='extrapolate')(r_middle_integ)
         Jredlu = interp1d(
@@ -210,4 +212,10 @@ def interpolate_integrator_quantities(self, att_S_ul, Jredlu,
                 r_middle, Jbluelu, fill_value='extrapolate')(r_middle_integ)
         e_dot_u = interp1d(
                 r_middle, e_dot_u, fill_value='extrapolate')(r_middle_integ)
+
+        # Set negative values from the extrapolation to zero
+        att_S_ul = att_S_ul.clip(0.)
+        Jbluelu = Jbluelu.clip(0.)
+        Jredlu = Jredlu.clip(0.)
+        e_dot_u = e_dot_u.clip(0.)
         return att_S_ul, Jredlu, Jbluelu, e_dot_u