avoid divide by 0 warnings

qpv-research-group · Sep 26, 2024 · f6eb284 · f6eb284
1 parent 8c9bead
commit f6eb284
Show file tree

Hide file tree

Showing 2 changed files with 6 additions and 3 deletions.
diff --git a/examples/pvk_Si_analytical_RT.py b/examples/pvk_Si_analytical_RT.py
@@ -55,7 +55,6 @@
 options.maximum_passes = 100
 options.n_rays = n_rays
 options.pol = 'u'
-options.use_numba = False
 
 front_text = regular_pyramids(52, True, 1,
                                 interface_layers=[Layer(100e-9, MgF2), Layer(1000e-9, Pvk)],

diff --git a/rayflare/ray_tracing/analytical_rt.py b/rayflare/ray_tracing/analytical_rt.py
@@ -597,7 +597,10 @@ def analytical_per_face(current_surf,
             refracted_rays[:, 2] / np.linalg.norm(refracted_rays,
                                                   axis=1))  # cos (global) of refracted ray
 
-        T_pol_per_it[:, N_interaction] = T_stack/T_prob[:,:,None]
+        # T_pol_per_it[:, N_interaction] = T_stack/T_prob[:,:,None]
+        T_pol_per_it[:, N_interaction] = np.divide(T_stack, T_prob[:, :, None], out=T_pol_per_it[:, N_interaction],
+                  where=T_prob[:, :, None] > 1e-15)
+
         cos_inc[reflected_direction[:, 2] > 0] = 0
         stop_it[
             np.all((np.all(reflected_direction[:, 2] > 0, axis=1), stop_it > N_interaction),
@@ -651,7 +654,8 @@ def analytical_per_face(current_surf,
     final_T_weights[final_T_weights < 0] = 0
     final_T_directions = np.array(final_T_directions)
     final_T_pol = np.array(final_T_pol)
-    final_T_pol = final_T_pol / np.sum(final_T_pol, -1)[:, :, None]
+    final_T_pol = np.divide(final_T_pol, np.sum(final_T_pol, -1)[:, :, None],
+                            out=final_T_pol, where=np.sum(final_T_pol, -1)[:, :, None] > 1e-15)
 
     # A_per_it has dimensions: (face, layer, interaction, wavelength)
     # sum over layers: