Don't change which reflections are used for refinement until the next macrocycle

5 files changed, 125 insertions, 113 deletions

diff --git a/src/facetron.c b/src/facetron.c
index 2bac7511..663a3236 100644
--- a/src/facetron.c
+++ b/src/facetron.c
@@ -101,7 +101,7 @@ static void refine_image(int mytask, void *tasks)
 		return;
 	}
 
-	spots = find_intersections(image, image->indexed_cell, &n, 0);
+	spots = find_intersections(image, image->indexed_cell, &n, 0, NULL);
 	dev = +INFINITY;
 	i = 0;
 	do {
@@ -162,7 +162,7 @@ static void integrate_image(int mytask, void *tasks)
 	}
 
 	/* Figure out which spots should appear in this pattern */
-	spots = find_intersections(image, image->indexed_cell, &n, 0);
+	spots = find_intersections(image, image->indexed_cell, &n, 0, NULL);
 
 	/* For each reflection, estimate the partiality */
 	for ( j=0; j<n; j++ ) {
diff --git a/src/geometry.c b/src/geometry.c
index 664817e0..09544120 100644
--- a/src/geometry.c
+++ b/src/geometry.c
@@ -117,8 +117,122 @@ static double partiality(double r1, double r2, double r)
 }
 
 
+static int check_reflection(struct image *image, double mres, int output,
+                            struct cpeak *cpeaks, int np,
+                            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)
+{
+	double xl, yl, zl;
+	double ds, ds_sq;
+	double rlow, rhigh;     /* "Excitation error" */
+	signed int p;           /* Panel number */
+	double xda, yda;        /* Position on detector */
+	int close, inside;
+	double part;            /* Partiality */
+	int clamp_low = 0;
+	int clamp_high = 0;
+	double bandwidth = image->bw;
+	double divergence = image->div;
+	double lambda = image->lambda;
+	double klow, kcen, khigh;    /* Wavenumber */
+	/* Bounding sphere for the shape transform approximation */
+	const double profile_cutoff = 0.02e9;  /* 0.02 nm^-1 */
+
+	/* "low" gives the largest Ewald sphere,
+	 * "high" gives the smallest Ewald sphere. */
+	klow = 1.0/(lambda - lambda*bandwidth/2.0);
+	kcen = 1.0/lambda;
+	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" */
+	if ( zl > profile_cutoff ) return 0;
+	xl = h*asx + k*bsx + l*csx;
+	yl = h*asy + k*bsy + l*csy;
+
+	/* Calculate reciprocal lattice point modulus (and square) */
+	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);
+	rhigh = excitation_error(xl, yl, zl, ds, khigh, +divergence);
+
+	/* Is the reciprocal lattice point close to either extreme of
+	 * the sphere, maybe just outside the "Ewald volume"? */
+	close = (fabs(rlow) < profile_cutoff)
+	     || (fabs(rhigh) < profile_cutoff);
+
+	/* Is the reciprocal lattice point somewhere between the
+	 * extremes of the sphere, i.e. inside the "Ewald volume"? */
+	inside = signbit(rlow) ^ signbit(rhigh);
+
+	/* Can't be both inside and close */
+	if ( inside ) close = 0;
+
+	/* Neither?  Skip it. */
+	if ( !(close || inside) ) return 0;
+
+	/* If the "lower" Ewald sphere is a long way away, use the
+	 * position at which the Ewald sphere would just touch the
+	 * reflection. */
+	if ( rlow < -profile_cutoff ) {
+		rlow = -profile_cutoff;
+		clamp_low = -1;
+	}
+	if ( rlow > +profile_cutoff ) {
+		rlow = +profile_cutoff;
+		clamp_low = +1;
+	}
+	/* Likewise the "higher" Ewald sphere */
+	if ( rhigh < -profile_cutoff ) {
+		rhigh = -profile_cutoff;
+		clamp_high = -1;
+	}
+	if ( rhigh > +profile_cutoff ) {
+		rhigh = +profile_cutoff;
+		clamp_high = +1;
+	}
+	/* The six possible combinations of clamp_{low,high} (including
+	 * zero) correspond to the six situations in Table 3 of Rossmann
+	 * et al. (1979). */
+
+	/* Calculate partiality and reject if too small */
+	part = partiality(rlow, rhigh, image->profile_radius);
+	if ( part < 0.1 ) return 0;
+
+	/* Locate peak on detector. */
+	p = locate_peak(xl, yl, zl, kcen, image->det, &xda, &yda);
+	if ( p == -1 ) return 0;
+
+	/* Add peak to list */
+	cpeaks[np].h = h;
+	cpeaks[np].k = k;
+	cpeaks[np].l = l;
+	cpeaks[np].x = xda;
+	cpeaks[np].y = yda;
+	cpeaks[np].r1 = rlow;
+	cpeaks[np].r2 = rhigh;
+	cpeaks[np].p = part;
+	cpeaks[np].clamp1 = clamp_low;
+	cpeaks[np].clamp2 = clamp_high;
+	np++;
+
+	if ( output ) {
+		printf("%3i %3i %3i %6f (at %5.2f,%5.2f) %5.2f\n",
+		       h, k, l, 0.0, xda, yda, part);
+	}
+
+	return 1;
+}
+
+
 struct cpeak *find_intersections(struct image *image, UnitCell *cell,
-                                 int *n, int output)
+                                 int *n, int output, struct cpeak *t)
 {
 	double asx, asy, asz;
 	double bsx, bsy, bsz;
@@ -128,12 +242,6 @@ struct cpeak *find_intersections(struct image *image, UnitCell *cell,
 	int hmax, kmax, lmax;
 	double mres;
 	signed int h, k, l;
-	double bandwidth = image->bw;
-	double divergence = image->div;
-	double lambda = image->lambda;
-	double klow, kcen, khigh;    /* Wavenumber */
-	/* Bounding sphere for the shape transform approximation */
-	const double profile_cutoff = 0.02e9;  /* 0.02 nm^-1 */
 
 	cpeaks = malloc(sizeof(struct cpeak)*MAX_CPEAKS);
 	if ( cpeaks == NULL ) {
@@ -150,112 +258,14 @@ struct cpeak *find_intersections(struct image *image, UnitCell *cell,
 	kmax = mres / modulus(bsx, bsy, bsz);
 	lmax = mres / modulus(csx, csy, csz);
 
-	/* "low" gives the largest Ewald sphere,
-	 * "high" gives the smallest Ewald sphere. */
-	klow = 1.0/(lambda - lambda*bandwidth/2.0);
-	kcen = 1.0/lambda;
-	khigh = 1.0/(lambda + lambda*bandwidth/2.0);
-
 	for ( h=-hmax; h<hmax; h++ ) {
 	for ( k=-kmax; k<kmax; k++ ) {
 	for ( l=-lmax; l<lmax; l++ ) {
-
-		double xl, yl, zl;
-		double ds, ds_sq;
-		double rlow, rhigh;     /* "Excitation error" */
-		signed int p;           /* Panel number */
-		double xda, yda;        /* Position on detector */
-		int close, inside;
-		double part;            /* Partiality */
-		int clamp_low = 0;
-		int clamp_high = 0;
-
 		/* Ignore central beam */
 		if ( (h==0) && (k==0) && (l==0) ) continue;
-
-		/* Get the coordinates of the reciprocal lattice point */
-		zl = h*asz + k*bsz + l*csz;
-		/* Throw out if it's "in front" */
-		if ( zl > profile_cutoff ) continue;
-		xl = h*asx + k*bsx + l*csx;
-		yl = h*asy + k*bsy + l*csy;
-
-		/* Calculate reciprocal lattice point modulus (and square) */
-		ds_sq = modulus_squared(xl, yl, zl);  /* d*^2 */
-		ds = sqrt(ds_sq);
-		if ( ds > mres ) continue;  /* Outside resolution range */
-
-		/* Calculate excitation errors */
-		rlow = excitation_error(xl, yl, zl, ds, klow, -divergence);
-		rhigh = excitation_error(xl, yl, zl, ds, khigh, +divergence);
-
-		/* Is the reciprocal lattice point close to either extreme of
-		 * the sphere, maybe just outside the "Ewald volume"? */
-		close = (fabs(rlow) < profile_cutoff)
-		     || (fabs(rhigh) < profile_cutoff);
-
-		/* Is the reciprocal lattice point somewhere between the
-		 * extremes of the sphere, i.e. inside the "Ewald volume"? */
-		inside = signbit(rlow) ^ signbit(rhigh);
-
-		/* Can't be both inside and close */
-		if ( inside ) close = 0;
-
-		/* Neither?  Skip it. */
-		if ( !(close || inside) ) continue;
-
-		/* If the "lower" Ewald sphere is a long way away, use the
-		 * position at which the Ewald sphere would just touch the
-		 * reflection. */
-		if ( rlow < -profile_cutoff ) {
-			rlow = -profile_cutoff;
-			clamp_low = -1;
-		}
-		if ( rlow > +profile_cutoff ) {
-			rlow = +profile_cutoff;
-			clamp_low = +1;
-		}
-		/* Likewise the "higher" Ewald sphere */
-		if ( rhigh < -profile_cutoff ) {
-			rhigh = -profile_cutoff;
-			clamp_high = -1;
-		}
-		if ( rhigh > +profile_cutoff ) {
-			rhigh = +profile_cutoff;
-			clamp_high = +1;
-		}
-		/* The six possible combinations of clamp_{low,high} (including
-		 * zero) correspond to the six situations in Table 3 of Rossmann
-		 * et al. (1979). */
-
-		/* Calculate partiality and reject if too small */
-		part = partiality(rlow, rhigh, image->profile_radius);
-		if ( part < 0.1 ) continue;
-
-		/* Locate peak on detector. */
-		p = locate_peak(xl, yl, zl, kcen, image->det, &xda, &yda);
-		if ( p == -1 ) continue;
-
-		/* Add peak to list */
-		cpeaks[np].h = h;
-		cpeaks[np].k = k;
-		cpeaks[np].l = l;
-		cpeaks[np].x = xda;
-		cpeaks[np].y = yda;
-		cpeaks[np].r1 = rlow;
-		cpeaks[np].r2 = rhigh;
-		cpeaks[np].p = part;
-		cpeaks[np].clamp1 = clamp_low;
-		cpeaks[np].clamp2 = clamp_high;
-		np++;
-
-		if ( output ) {
-			printf("%3i %3i %3i %6f (at %5.2f,%5.2f) %5.2f\n",
-			       h, k, l, 0.0, xda, yda, part);
-		}
-
+		np += check_reflection(image, mres, output, cpeaks, np, h, k, l,
+		                       asx,asy,asz,bsx,bsy,bsz,csx,csy,csz);
 		if ( np == MAX_CPEAKS ) goto out;
-
 	}
 	}
 	}
diff --git a/src/geometry.h b/src/geometry.h
index 0782a3e1..03fde52f 100644
--- a/src/geometry.h
+++ b/src/geometry.h
@@ -18,7 +18,7 @@
 #endif
 
 extern struct cpeak *find_intersections(struct image *image, UnitCell *cell,
-                                        int *n, int output);
+                                        int *n, int output, struct cpeak *t);
 
 extern double integrate_all(struct image *image, struct cpeak *cpeaks, int n);
 
diff --git a/src/post-refinement.c b/src/post-refinement.c
index e2fa6d50..c4ca03cc 100644
--- a/src/post-refinement.c
+++ b/src/post-refinement.c
@@ -165,6 +165,7 @@ double pr_iterate(struct image *image, double *i_full, const char *sym,
 	gsl_vector *shifts;
 	int h, param;
 	struct cpeak *spots = *pspots;
+	struct cpeak *spots_old;
 
 	M = gsl_matrix_calloc(NUM_PARAMS, NUM_PARAMS);
 	v = gsl_vector_calloc(NUM_PARAMS);
@@ -241,8 +242,9 @@ double pr_iterate(struct image *image, double *i_full, const char *sym,
 	gsl_vector_free(v);
 	gsl_vector_free(shifts);
 
-	free(spots);
-	spots = find_intersections(image, image->indexed_cell, n, 0);
+	spots_old = spots;
+	spots = find_intersections(image, image->indexed_cell, n, 0, spots_old);
+	free(spots_old);
 	*pspots = spots;
 	return mean_partial_dev(image, spots, *n, sym, i_full, NULL);
 }
diff --git a/src/templates.c b/src/templates.c
index 4bd85a76..57434473 100644
--- a/src/templates.c
+++ b/src/templates.c
@@ -172,7 +172,7 @@ IndexingPrivate *generate_templates(UnitCell *cell, const char *filename,
 
 		cell_rot = rotate_cell(cell, omega, phi, 0.0);
 
-		cpeaks = find_intersections(&image, cell_rot, &n, 0);
+		cpeaks = find_intersections(&image, cell_rot, &n, 0, NULL);
 		if ( cpeaks == NULL ) {
 			ERROR("Template calculation failed.\n");
 			return NULL;