aboutsummaryrefslogtreecommitdiff
diff options
context:
space:
mode:
authorThomas White <taw@physics.org>2011-07-07 16:10:41 +0200
committerThomas White <taw@physics.org>2012-02-22 15:27:32 +0100
commit8ac0fea56a5ef1c0405af178d5da0038052fa045 (patch)
tree8e2664e7d7518c6e83eae49a369fb3097953f129
parent2c932f6cef1d273ef7f71fddd1582f9fb8bbe8d8 (diff)
Don't waste time predicting reflections that don't appear in the pattern
-rw-r--r--src/geometry.c76
1 files changed, 30 insertions, 46 deletions
diff --git a/src/geometry.c b/src/geometry.c
index b338da56..f94980d3 100644
--- a/src/geometry.c
+++ b/src/geometry.c
@@ -122,12 +122,11 @@ static double partiality(double r1, double r2, double r)
}
-static int check_reflection(struct image *image, double mres,
- RefList *reflections,
- signed int h, signed int k, signed int l,
- double asx, double asy, double asz,
- double bsx, double bsy, double bsz,
- double csx, double csy, double csz)
+static Reflection *check_reflection(struct image *image, RefList *reflections,
+ signed int h, signed int k, signed int l,
+ double asx, double asy, double asz,
+ double bsx, double bsy, double bsz,
+ double csx, double csy, double csz)
{
const int output = 0;
double xl, yl, zl;
@@ -154,13 +153,12 @@ static int check_reflection(struct image *image, double mres,
/* 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 0;
+ if ( zl > image->profile_radius ) return NULL;
xl = h*asx + k*bsx + l*csx;
yl = h*asy + k*bsy + l*csy;
ds_sq = modulus_squared(xl, yl, zl); /* d*^2 */
ds = sqrt(ds_sq);
- if ( ds > mres ) return 0; /* Outside resolution range */
/* Calculate excitation errors */
rlow = excitation_error(xl, yl, zl, ds, klow, -divergence);
@@ -179,7 +177,7 @@ static int check_reflection(struct image *image, double mres,
if ( inside ) close = 0;
/* Neither? Skip it. */
- if ( !(close || inside) ) return 0;
+ if ( !(close || inside) ) return NULL;
/* If the "lower" Ewald sphere is a long way away, use the
* position at which the Ewald sphere would just touch the
@@ -211,7 +209,7 @@ static int check_reflection(struct image *image, double mres,
/* Locate peak on detector. */
p = locate_peak(xl, yl, zl, kcen, image->det, &xda, &yda);
- if ( p == -1 ) return 0;
+ if ( p == -1 ) return NULL;
/* Add peak to list */
refl = add_refl(reflections, h, k, l);
@@ -225,7 +223,7 @@ static int check_reflection(struct image *image, double mres,
h, k, l, 0.0, xda, yda, part);
}
- return 1;
+ return refl;
}
@@ -264,7 +262,7 @@ RefList *find_intersections(struct image *image, UnitCell *cell)
for ( h=-hmax; h<=hmax; h++ ) {
for ( k=-kmax; k<=kmax; k++ ) {
for ( l=-lmax; l<=lmax; l++ ) {
- check_reflection(image, mres, reflections, h, k, l,
+ check_reflection(image, reflections, h, k, l,
asx,asy,asz,bsx,bsy,bsz,csx,csy,csz);
}
}
@@ -282,63 +280,49 @@ void predict_corresponding_reflections(struct image *image, const char *sym,
Reflection *refl;
RefListIterator *iter;
RefList *predicted;
+ double asx, asy, asz;
+ double bsx, bsy, bsz;
+ double csx, csy, csz;
*n_expected = 0;
*n_found = 0;
*n_notfound = 0;
- predicted = find_intersections(image, image->indexed_cell);
+ cell_get_reciprocal(image->indexed_cell, &asx, &asy, &asz,
+ &bsx, &bsy, &bsz, &csx, &csy, &csz);
- for ( refl = first_refl(predicted, &iter);
+ /* Scratch list to give check_reflection() something to add to */
+ predicted = reflist_new();
+
+ for ( refl = first_refl(image->reflections, &iter);
refl != NULL;
refl = next_refl(refl, iter) )
{
-
- Reflection *p_peak;
+ Reflection *vals;
double r1, r2, p, x, y;
signed int h, k, l;
- signed int ha, ka, la;
int clamp1, clamp2;
- int found = 0;
- /* Get predicted indices and location */
- get_indices(refl, &h, &k, &l);
- get_detector_pos(refl, &x, &y);
if ( n_expected != NULL ) (*n_expected)++;
- /* Get the asymmetric indices, with which the reflection in the
- * image's list will be indexed. */
- get_asymm(h, k, l, &ha, &ka, &la, sym);
-
- /* Look for this reflection in the pattern */
- p_peak = find_refl(image->reflections, ha, ka, la);
- do {
-
- signed int hs, ks, ls;
-
- if ( p_peak != NULL ) {
- get_symmetric_indices(p_peak, &hs, &ks, &ls);
- if ( (hs==h) && (ks==k) && (ls==l) ) found = 1;
- }
+ get_symmetric_indices(refl, &h, &k, &l);
- if ( !found && (p_peak != NULL ) ) {
- p_peak = next_found_refl(p_peak);
- }
+ vals = check_reflection(image, predicted, h, k, l,
+ asx,asy,asz,bsx,bsy,bsz,csx,csy,csz);
- } while ( !found && (p_peak != NULL) );
- if ( !found ) {
- if (n_notfound != NULL) (*n_notfound)++;
+ if ( vals == NULL ) {
+ (*n_notfound)++;
continue;
}
- if ( n_found != NULL ) (*n_found)++;
+ (*n_found)++;
/* Transfer partiality stuff */
- get_partial(refl, &r1, &r2, &p, &clamp1, &clamp2);
- set_partial(p_peak, r1, r2, p, clamp1, clamp2);
+ get_partial(vals, &r1, &r2, &p, &clamp1, &clamp2);
+ set_partial(refl, r1, r2, p, clamp1, clamp2);
/* Transfer detector location */
- get_detector_pos(refl, &x, &y);
- set_detector_pos(p_peak, 0.0, x, y);
+ get_detector_pos(vals, &x, &y);
+ set_detector_pos(refl, 0.0, x, y);
}