diff --git a/image_analysis/indexing/AnalyzeIndexing.cpp b/image_analysis/indexing/AnalyzeIndexing.cpp
index 02181045..8d02357b 100644
--- a/image_analysis/indexing/AnalyzeIndexing.cpp
+++ b/image_analysis/indexing/AnalyzeIndexing.cpp
@@ -29,9 +29,13 @@ namespace {
 
     // Wrap to [-180, 180] (useful for residuals)
     inline float wrap_deg_pm180(float deg) {
-        while (deg > 180.0f) deg -= 360.0f;
-        while (deg < -180.0f) deg += 360.0f;
-        return deg;
+        if (!std::isfinite(deg)) return std::numeric_limits<float>::quiet_NaN();; // or std::nullopt upstream
+
+        deg = std::fmod(deg + 180.0f, 360.0f);
+        if (!std::isfinite(deg))
+            return std::numeric_limits<float>::quiet_NaN();
+        if (deg < 0) deg += 360.0f;
+        return deg - 180.0f;
     }
 
      // XDS convention: zeta = |m2 · e1| where e1 = (S × S0) / |S × S0|
@@ -74,7 +78,7 @@ namespace {
         double ll = 0.0;
         for (size_t i = 0; i < tau_values.size(); ++i) {
             float R = R_fraction(tau_values[i], delta_phi, sigma_M, zeta_values[i]);
-            if (R > 1e-30f) {
+            if (std::isfinite(R) && R > 1e-30f) {
                 ll += std::log(static_cast<double>(R));
             } else {
                 ll += -70.0;  // Large penalty for zero probability
@@ -98,8 +102,8 @@ namespace {
         float d = a + golden * (b - a);
 
         const float tol = 1e-6f;
-
-        while (std::fabs(b - a) > tol) {
+        int iter = 0;
+        while (std::fabs(b - a) > tol && iter++ < 100) {
             double fc = log_likelihood(tau_values, zeta_values, delta_phi, c);
             double fd = log_likelihood(tau_values, zeta_values, delta_phi, d);
 
@@ -126,7 +130,7 @@ namespace {
             return 0;
 
         const float rhs = -D / R;
-        if (rhs < -1.0f || rhs > 1.0f)
+        if (!std::isfinite(rhs) || rhs < -1.0f || rhs > 1.0f)
             return 0;
 
         const float phi_ref = std::atan2(B, A);
@@ -137,9 +141,24 @@ namespace {
 
         const float two_pi = 2.0f * static_cast<float>(M_PI);
         auto shift_near = [&](float x) {
-            while (x < phi0 - two_pi) x += two_pi;
-            while (x > phi1 + two_pi) x -= two_pi;
-            return x;
+            if (!std::isfinite(x))
+                return std::numeric_limits<float>::quiet_NaN();
+
+            const float span_center = 0.5f * (phi0 + phi1);
+
+            // Bring x close to the interval center using modulo 2π
+            float shifted = x - span_center;
+            shifted = std::fmod(shifted, two_pi);
+            if (!std::isfinite(shifted))
+                return std::numeric_limits<float>::quiet_NaN();
+
+            // fmod can return negative values; normalize to [-π, π]
+            if (shifted < -static_cast<float>(M_PI))
+                shifted += two_pi;
+            else if (shifted > static_cast<float>(M_PI))
+                shifted -= two_pi;
+
+            return shifted + span_center;
         };
 
         s1 = shift_near(s1);
@@ -232,7 +251,7 @@ namespace {
 
             // Find predicted phi angle
             const auto phi_pred_opt = predict_phi_deg_local(pstar, S0, m2, 0.0f, axis.GetWedge_deg());
-            if (!phi_pred_opt.has_value())
+            if (!phi_pred_opt.has_value() || !std::isfinite(phi_pred_opt.value()))
                 continue;
 
             // τ (tau) = angular deviation from Bragg position to center of oscillation range
@@ -257,7 +276,7 @@ namespace {
             float zeta = calc_zeta(S, S0, m2);
 
             // Filter out reflections with very small zeta (poorly determined)
-            if (zeta < 0.1f)
+            if (!std::isfinite(zeta) || !std::isfinite(tau_rad) || zeta < 0.1f)
                 continue;
 
             tau_values.push_back(tau_rad);