aboutsummaryrefslogtreecommitdiff
diff options
context:
space:
mode:
-rw-r--r--src/get_hkl.c63
1 files changed, 35 insertions, 28 deletions
diff --git a/src/get_hkl.c b/src/get_hkl.c
index 0e556a4e..454737b0 100644
--- a/src/get_hkl.c
+++ b/src/get_hkl.c
@@ -3,11 +3,11 @@
*
* Small program to manipulate reflection lists
*
- * Copyright © 2012 Deutsches Elektronen-Synchrotron DESY,
+ * Copyright © 2013 Deutsches Elektronen-Synchrotron DESY,
* a research centre of the Helmholtz Association.
*
* Authors:
- * 2009-2012 Thomas White <taw@physics.org>
+ * 2009-2013 Thomas White <taw@physics.org>
*
* This file is part of CrystFEL.
*
@@ -173,43 +173,36 @@ static RefList *twin_reflections(RefList *in, int need_all_parts,
RefListIterator *iter;
RefList *out;
SymOpMask *m;
-
- if ( !is_subgroup(holo, mero) ) {
- ERROR("%s is not a subgroup of %s!\n", symmetry_name(mero),
- symmetry_name(holo));
- return NULL;
- }
+ int n;
out = reflist_new();
+
+ /* No need to free and reallocate this for every reflection */
m = new_symopmask(holo);
for ( refl = first_refl(in, &iter);
refl != NULL;
- refl = next_refl(refl, iter) ) {
-
+ refl = next_refl(refl, iter) )
+ {
double total, sigma;
- int multi, nbits;
+ int multi;
signed int h, k, l;
- int n, j;
+ int j;
int skip;
+ /* Figure out where to put the twinned version, and check it's
+ * not there already. */
get_indices(refl, &h, &k, &l);
-
- /* There is a many-to-one correspondence between reflections
- * in the merohedral and holohedral groups. Do the calculation
- * only once for each reflection in the holohedral group, which
- * contains fewer reflections.
- */
get_asymm(holo, h, k, l, &h, &k, &l);
if ( find_refl(out, h, k, l) != NULL ) continue;
+ special_position(holo, m, h, k, l);
+ n = num_equivs(holo, m);
+
total = 0.0;
sigma = 0.0;
multi = 0;
skip = 0;
- nbits = 0;
- special_position(holo, m, h, k, l);
- n = num_equivs(holo, m);
for ( j=0; j<n; j++ ) {
@@ -219,6 +212,7 @@ static RefList *twin_reflections(RefList *in, int need_all_parts,
int r;
get_equiv(holo, m, j, h, k, l, &he, &ke, &le);
+ get_asymm(mero, he, ke, le, &he, &ke, &le);
/* Do we have this reflection?
* We might not have the particular (merohedral)
@@ -240,20 +234,33 @@ static RefList *twin_reflections(RefList *in, int need_all_parts,
}
- part = find_refl(in, hu, ku, lu);
+ if ( r ) {
- total += get_intensity(part);
- sigma += pow(get_esd_intensity(part), 2.0);
- multi += get_redundancy(part);
- nbits++;
+ double i, sigi;
+ int mult;
+
+ part = find_refl(in, hu, ku, lu);
+
+ i = get_intensity(part);
+ sigi = get_esd_intensity(part);
+ mult = get_redundancy(part);
+
+ total += mult*i;
+ sigma += pow(sigi*mult, 2.0);
+ multi += mult;
+
+ set_intensity(part, 0.0);
+ set_esd_intensity(part, 0.0);
+ set_redundancy(part, 0);
+ }
}
if ( !skip ) {
Reflection *new = add_refl(out, h, k, l);
- set_intensity(new, total/nbits);
- set_esd_intensity(new, sqrt(sigma)/nbits);
+ set_intensity(new, total/multi);
+ set_esd_intensity(new, sqrt(sigma)/multi);
set_redundancy(new, multi);
}