diff --git a/image_analysis/bragg_prediction/BraggPrediction.cpp b/image_analysis/bragg_prediction/BraggPrediction.cpp
index fe6c2b6d..0e982ad6 100644
--- a/image_analysis/bragg_prediction/BraggPrediction.cpp
+++ b/image_analysis/bragg_prediction/BraggPrediction.cpp
@@ -37,8 +37,11 @@ int BraggPrediction::Calc(const DiffractionExperiment &experiment, const Crystal
     float pixel_size = geom.GetPixelSize_mm();
     float F = det_distance / pixel_size;
 
-    int i = 0;
+    const float epsilon = 1e-5f;
+    const float s0_sq = S0 * S0;
+    const float rad_to_deg = 180.0f / static_cast<float>(M_PI);
 
+    int i = 0;
 
     for (int h = -settings.max_hkl; h <= settings.max_hkl; h++) {
         // Precompute A* h contribution
@@ -74,6 +77,31 @@ int BraggPrediction::Calc(const DiffractionExperiment &experiment, const Crystal
                 float dist_ewald_sphere = std::fabs(S_len - one_over_wavelength);
 
                 if (dist_ewald_sphere <= settings.ewald_dist_cutoff ) {
+                    const float s0_p0 = S0.x * recip_x + S0.y * recip_y + S0.z * recip_z;
+                    const float val = s0_sq * recip_sq - s0_p0 * s0_p0;
+
+                    float delta_phi_deg = NAN;
+                    if (std::fabs(val) >= epsilon && s0_sq > epsilon) {
+                        const float a_num = (s0_sq - 0.25f * recip_sq) * recip_sq;
+                        if (a_num >= 0.0f) {
+                            const float A = std::sqrt(a_num / val);
+                            const float B = (A * s0_p0 + 0.5f * recip_sq) / s0_sq;
+
+                            const float p_star_x = A * recip_x - B * S0.x;
+                            const float p_star_y = A * recip_y - B * S0.y;
+                            const float p_star_z = A * recip_z - B * S0.z;
+
+                            const float p_star_sq = p_star_x * p_star_x + p_star_y * p_star_y + p_star_z * p_star_z;
+                            const float denom = std::sqrt(p_star_sq * recip_sq);
+
+                            if (denom >= epsilon) {
+                                float c = (p_star_x * recip_x + p_star_y * recip_y + p_star_z * recip_z) / denom;
+                                c = std::fmax(-1.0f, std::fmin(1.0f, c));
+                                delta_phi_deg = std::acos(c) * rad_to_deg;
+                            }
+                        }
+                    }
+
                     // Inlined RecipToDector with rot1 and rot2 (rot3 = 0)
                     // Apply rotation matrix transpose
                     float S_rot_x = rot[0] * S_x + rot[1] * S_y + rot[2] * S_z;
@@ -96,7 +124,7 @@ int BraggPrediction::Calc(const DiffractionExperiment &experiment, const Crystal
                         .h = h,
                         .k = k,
                         .l = l,
-                        .delta_phi_deg = NAN,
+                        .delta_phi_deg = delta_phi_deg,
                         .predicted_x = x,
                         .predicted_y = y,
                         .d = d,
diff --git a/image_analysis/bragg_prediction/BraggPredictionGPU.cu b/image_analysis/bragg_prediction/BraggPredictionGPU.cu
index ea8080d2..e561542e 100644
--- a/image_analysis/bragg_prediction/BraggPredictionGPU.cu
+++ b/image_analysis/bragg_prediction/BraggPredictionGPU.cu
@@ -10,6 +10,36 @@
 namespace {
     __device__ inline bool is_odd(int v) { return (v & 1) != 0; }
 
+    __device__ inline float angle_from_ewald_sphere_deg(const Coord &S0, float recip_x, float recip_y, float recip_z, float recip_sq) {
+        const float epsilon = 1e-5f;
+        const float rad_to_deg = 180.0f / static_cast<float>(M_PI);
+
+        const float s0_sq = S0.x * S0.x + S0.y * S0.y + S0.z * S0.z;
+        const float s0_p0 = S0.x * recip_x + S0.y * recip_y + S0.z * recip_z;
+        const float val = s0_sq * recip_sq - s0_p0 * s0_p0;
+
+        if (fabsf(val) < epsilon || s0_sq < epsilon) return NAN;
+
+        const float a_num = (s0_sq - 0.25f * recip_sq) * recip_sq;
+        if (a_num < 0.0f) return NAN;
+
+        const float A = sqrtf(a_num / val);
+        const float B = (A * s0_p0 + 0.5f * recip_sq) / s0_sq;
+
+        const float p_star_x = A * recip_x - B * S0.x;
+        const float p_star_y = A * recip_y - B * S0.y;
+        const float p_star_z = A * recip_z - B * S0.z;
+
+        const float p_star_sq = p_star_x * p_star_x + p_star_y * p_star_y + p_star_z * p_star_z;
+        const float denom = sqrtf(p_star_sq * recip_sq);
+        if (denom < epsilon) return NAN;
+
+        float c = (p_star_x * recip_x + p_star_y * recip_y + p_star_z * recip_z) / denom;
+        c = fmaxf(-1.0f, fminf(1.0f, c));
+
+        return acosf(c) * rad_to_deg;
+    }
+
     __device__ inline bool compute_reflection(const KernelConsts &C, int h, int k, int l, Reflection &out) {
         if (h == 0 && k == 0 && l == 0)
             return false;
@@ -77,7 +107,7 @@ namespace {
         out.h = h;
         out.k = k;
         out.l = l;
-        out.delta_phi_deg = NAN;
+        out.delta_phi_deg = angle_from_ewald_sphere_deg(C.S0, recip_x, recip_y, recip_z, recip_sq);
         out.predicted_x = x;
         out.predicted_y = y;
         out.d = 1.0f / sqrtf(recip_sq);