From 4b83b1cc6d3a38d8a8f16de5fe5f645b5e6498f0 Mon Sep 17 00:00:00 2001 From: Thomas White Date: Thu, 1 Dec 2011 11:08:22 +0100 Subject: Fix incorrect culling for reciprocal lattice peaks --- libcrystfel/src/geometry.c | 16 +++++++++------- 1 file changed, 9 insertions(+), 7 deletions(-) diff --git a/libcrystfel/src/geometry.c b/libcrystfel/src/geometry.c index 485abba3..1aef080f 100644 --- a/libcrystfel/src/geometry.c +++ b/libcrystfel/src/geometry.c @@ -143,7 +143,7 @@ static Reflection *check_reflection(struct image *image, double lambda = image->lambda; double klow, kcen, khigh; /* Wavenumber */ Reflection *refl; - double tt; + double tt, psi; /* "low" gives the largest Ewald sphere, * "high" gives the smallest Ewald sphere. */ @@ -152,18 +152,20 @@ static Reflection *check_reflection(struct image *image, khigh = 1.0/(lambda + lambda*bandwidth/2.0); /* Get the coordinates of the reciprocal lattice point */ - zl = h*asz + k*bsz + l*csz; - /* Throw out if it's "in front". A tiny bit "in front" is OK. */ - if ( zl > image->profile_radius ) return NULL; xl = h*asx + k*bsx + l*csx; yl = h*asy + k*bsy + l*csy; - - tt = angle_between(0.0, 0.0, 1.0, xl, yl, zl+kcen); - if ( tt > deg2rad(90.0) ) return NULL; + zl = h*asz + k*bsz + l*csz; ds_sq = modulus_squared(xl, yl, zl); /* d*^2 */ ds = sqrt(ds_sq); + /* Simple (fast) check to reject reflection if it's "in front" */ + psi = atan2(zl, sqrt(xl*xl + yl*yl)); + if ( psi - atan2(image->profile_radius, ds) > divergence ) return NULL; + + tt = angle_between(0.0, 0.0, 1.0, xl, yl, zl+kcen); + if ( tt > deg2rad(90.0) ) return NULL; + /* Calculate excitation errors */ rlow = excitation_error(xl, yl, zl, ds, klow, -divergence/2.0, tt); rhigh = excitation_error(xl, yl, zl, ds, khigh, +divergence/2.0, tt); -- cgit v1.2.3