Force Friedel's Law when measuring intensities

1 files changed, 97 insertions, 55 deletions

diff --git a/src/intensities.c b/src/intensities.c
index 86f3de1..c933d3e 100644
--- a/src/intensities.c
+++ b/src/intensities.c
@@ -21,13 +21,99 @@
 #include "utils.h"
 #include "basis.h"
 
+/* Add a Friedel equivalent if missing, otherwise set the
+ *  intensities of both equivalents to the larger intensity */
+static void intensity_friedelise(Reflection *ref, ReflectionList *list,
+				 int *n_fried)
+{
+	Reflection *equiv;
+
+	equiv = reflectionlist_find(list, -ref->h, -ref->k, -ref->l);
+
+	/* Friedel equivalent exists? */
+	if ( equiv == NULL ) {
+
+		/* No, add it */
+		equiv = reflection_add(list, -ref->x, -ref->y, -ref->z,
+					ref->intensity, REFLECTION_GENERATED);
+
+		equiv->h = -ref->h;
+		equiv->k = -ref->k;
+		equiv->l = -ref->l;
+
+		(*n_fried)++;
+
+	} else {
+
+		/* Yes, use largest intensity */
+		if ( equiv->intensity > ref->intensity )
+			ref->intensity = equiv->intensity;
+		else
+			equiv->intensity = ref->intensity;
+
+	}
+}
+
+static void intensity_measure(signed int h, signed int k, signed int l,
+				ImageFeature *feature, ReflectionList *list,
+				double *max,
+				int *n_meas, int *n_dupl, int *n_fried)
+{
+	if ( (h!=0) || (k!=0) || (l!=0) ) {
+
+		double intensity;
+		Reflection *ref;
+
+		intensity = feature->partner->intensity;
+
+		ref = reflectionlist_find(list, h, k, l);
+
+		/* Do we already have a value for this reflection? */
+		if ( ref == NULL ) {
+
+			/* No, add it */
+			ref = reflection_add(list,
+					     feature->reflection->x,
+					     feature->reflection->y,
+					     feature->reflection->z,
+					     intensity, REFLECTION_GENERATED);
+
+			ref->h = h;
+			ref->k = k;
+			ref->l = l;
+
+			if ( intensity > *max ) *max = intensity;
+
+			(*n_meas)++;
+
+		} else {
+
+			/* Yes, record only if the new value is larger */
+			if ( intensity > ref->intensity ) {
+
+				ref->x = feature->reflection->x;
+				ref->y = feature->reflection->y;
+				ref->z = feature->reflection->z;
+				ref->intensity = intensity;
+
+			}
+
+			(*n_dupl)++;
+
+		}
+
+		intensity_friedelise(ref, list, n_fried);
+
+	}
+}
+
 /* Extract integrated reflection intensities by estimating the spike function
  * based on the observed intensity and the calculated excitation error from
  * the lattice refinement.  Easy.	*/
 void intensities_extract(ControlContext *ctx)
 {
 	int i, j;
-	int n_meas, n_dupl, n_notf;
+	int n_meas, n_dupl, n_notf, n_fried;
 	double max;
 	Reflection *reflection;
 
@@ -40,16 +126,19 @@ void intensities_extract(ControlContext *ctx)
 	n_meas = 0;
 	n_dupl = 0;
 	n_notf = 0;
+	n_fried = 0;
 	max = 0;
 	for ( i=0; i<ctx->images->n_images; i++ ) {
 
 		ImageRecord *image;
 
+		/* Reproject this pattern, if necessary */
 		image = &ctx->images->images[i];
 		if ( image->rflist == NULL )
 			image->rflist = reproject_get_reflections(image,
 							ctx->cell_lattice);
 
+		/* For each reprojected feature */
 		for ( j=0; j<image->rflist->n_features; j++ ) {
 
 			ImageFeature *feature;
@@ -61,61 +150,13 @@ void intensities_extract(ControlContext *ctx)
 			k = feature->reflection->k;
 			l = feature->reflection->l;
 
-			if ( feature->partner != NULL ) {
-
-				if ( (h!=0) || (k!=0) || (l!=0) ) {
-
-					double intensity;
-					Reflection *ref;
-
-					intensity = feature->partner->intensity;
-
-					ref = reflectionlist_find(
-						      ctx->integrated, h, k, l);
-
-					if ( ref == NULL ) {
-
-						Reflection *new;
-
-						new = reflection_add(
-							ctx->integrated,
-							feature->reflection->x,
-							feature->reflection->y,
-							feature->reflection->z,
-							intensity,
-							REFLECTION_GENERATED);
-
-						new->h = h;
-						new->k = k;
-						new->l = l;
-
-						if ( intensity > max )
-							max = intensity;
-
-						n_meas++;
-
-					} else {
-
-						if ( intensity > ref->intensity ) {
-
-							ref->x = feature->reflection->x;
-							ref->y = feature->reflection->y;
-							ref->z = feature->reflection->z;
-							ref->intensity = intensity;
-
-						}
-
-						n_dupl++;
-
-					}
-
-				}
-
-			} else {
-				//printf("IN: %3i %3i %3i not found\n", h, k, l);
+			/* Corresponds to a measured reflection? */
+			if ( feature->partner != NULL )
+				intensity_measure(h, k, l, feature,
+						  ctx->integrated, &max,
+						  &n_meas, &n_dupl, &n_fried);
+			else
 				n_notf++;
-			}
-
 		}
 
 	}
@@ -129,6 +170,7 @@ void intensities_extract(ControlContext *ctx)
 
 	printf("IN: %5i intensities measured\n", n_meas);
 	printf("IN: %5i duplicated measurements\n", n_dupl);
+	printf("IN: %5i reflections generated by Friedel's law\n", n_fried);
 	printf("IN: %5i predicted reflections not found\n", n_notf);
 }