From fae05a5cfebec279f47a7463a9813643ac08b11a Mon Sep 17 00:00:00 2001
From: Filip Leonarski <filip.leonarski@psi.ch>
Date: Sat, 3 Jan 2026 15:02:26 +0100
Subject: [PATCH] BraggPrediction: Minor improvements - predict from -max_hkl
 to +max_hkl (incl. +max_hkl); small fixed in notation

---
 .../bragg_prediction/BraggPrediction.cpp      | 12 +++++------
 .../bragg_prediction/BraggPredictionGPU.cu    | 20 +++++++++----------
 2 files changed, 16 insertions(+), 16 deletions(-)

diff --git a/image_analysis/bragg_prediction/BraggPrediction.cpp b/image_analysis/bragg_prediction/BraggPrediction.cpp
index 00f2389b..ab7eb892 100644
--- a/image_analysis/bragg_prediction/BraggPrediction.cpp
+++ b/image_analysis/bragg_prediction/BraggPrediction.cpp
@@ -39,19 +39,19 @@ int BraggPrediction::Calc(const DiffractionExperiment &experiment, const Crystal
 
     int i = 0;
 
-    for (int h = -settings.max_hkl; h < settings.max_hkl; h++) {
+    for (int h = -settings.max_hkl; h <= settings.max_hkl; h++) {
         // Precompute A* h contribution
         const float Ah_x = Astar.x * h;
         const float Ah_y = Astar.y * h;
         const float Ah_z = Astar.z * h;
 
-        for (int k = -settings.max_hkl; k < settings.max_hkl; k++) {
+        for (int k = -settings.max_hkl; k <= settings.max_hkl; k++) {
             // Accumulate B* k contribution
             const float AhBk_x = Ah_x + Bstar.x * k;
             const float AhBk_y = Ah_y + Bstar.y * k;
             const float AhBk_z = Ah_z + Bstar.z * k;
 
-            for (int l = -settings.max_hkl; l < settings.max_hkl; l++) {
+            for (int l = -settings.max_hkl; l <= settings.max_hkl; l++) {
                 if (systematic_absence(h, k, l, settings.centering))
                     continue;
 
@@ -62,8 +62,8 @@ int BraggPrediction::Calc(const DiffractionExperiment &experiment, const Crystal
                 float recip_y = AhBk_y + Cstar.y * l;
                 float recip_z = AhBk_z + Cstar.z * l;
 
-                float dot = recip_x * recip_x + recip_y * recip_y + recip_z * recip_z;
-                if (dot > one_over_dmax_sq)
+                float recip_sq = recip_x * recip_x + recip_y * recip_y + recip_z * recip_z;
+                if (recip_sq > one_over_dmax_sq)
                     continue;
 
                 float S_x = recip_x + S0.x;
@@ -90,7 +90,7 @@ int BraggPrediction::Calc(const DiffractionExperiment &experiment, const Crystal
                     if ((x < 0) || (x >= det_width_pxl) || (y < 0) || (y >= det_height_pxl))
                         continue;
 
-                    float d = 1.0f / sqrtf(dot);
+                    float d = 1.0f / sqrtf(recip_sq);
                     reflections[i] = Reflection{
                         .h = h,
                         .k = k,
diff --git a/image_analysis/bragg_prediction/BraggPredictionGPU.cu b/image_analysis/bragg_prediction/BraggPredictionGPU.cu
index 7ac8f1fb..b80563d4 100644
--- a/image_analysis/bragg_prediction/BraggPredictionGPU.cu
+++ b/image_analysis/bragg_prediction/BraggPredictionGPU.cu
@@ -55,14 +55,14 @@ namespace {
         float AhBk_x = Ah_x + C.Bstar.x * k;
         float AhBk_y = Ah_y + C.Bstar.y * k;
         float AhBk_z = Ah_z + C.Bstar.z * k;
-        float rx = AhBk_x + C.Cstar.x * l;
-        float ry = AhBk_y + C.Cstar.y * l;
-        float rz = AhBk_z + C.Cstar.z * l;
-        float dot = rx * rx + ry * ry + rz * rz;
-        if (dot > C.one_over_dmax_sq) return false;
-        float Sx = rx + C.S0.x;
-        float Sy = ry + C.S0.y;
-        float Sz = rz + C.S0.z;
+        float recip_x = AhBk_x + C.Cstar.x * l;
+        float recip_y = AhBk_y + C.Cstar.y * l;
+        float recip_z = AhBk_z + C.Cstar.z * l;
+        float recip_sq = recip_x * recip_x + recip_y * recip_y + recip_z * recip_z;
+        if (recip_sq > C.one_over_dmax_sq) return false;
+        float Sx = recip_x + C.S0.x;
+        float Sy = recip_y + C.S0.y;
+        float Sz = recip_z + C.S0.z;
         float S_len = sqrtf(Sx * Sx + Sy * Sy + Sz * Sz);
         float dist_ewald = fabsf(S_len - C.one_over_wavelength);
         if (dist_ewald > C.ewald_cutoff) return false;
@@ -79,7 +79,7 @@ namespace {
         out.l = l;
         out.predicted_x = x;
         out.predicted_y = y;
-        out.d = 1.0f / sqrtf(dot);
+        out.d = 1.0f / sqrtf(recip_sq);
         out.dist_ewald = dist_ewald;
         return true;
     }
@@ -89,7 +89,7 @@ namespace {
                                     int max_reflections,
                                     Reflection *__restrict__ out,
                                     int *__restrict__ counter) {
-        int range = 2 * max_hkl;
+        int range = 2 * max_hkl + 1;
         int hi = blockIdx.x * blockDim.x + threadIdx.x;
         int ki = blockIdx.y * blockDim.y + threadIdx.y;
         int li = blockIdx.z * blockDim.z + threadIdx.z;