BraggPrediction: Minor improvements - predict from -max_hkl to +max_hkl (incl. +max_hkl); small fixed in notation

This commit is contained in:
2026-01-03 15:02:26 +01:00
parent 68e99afc9e
commit fae05a5cfe
2 changed files with 16 additions and 16 deletions

View File

@@ -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,

View File

@@ -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;