More geometry work (fixes a previous commit)

2 files changed, 33 insertions, 43 deletions

diff --git a/src/geometry.c b/src/geometry.c
index 2b94e7ba..1f21c96f 100644
--- a/src/geometry.c
+++ b/src/geometry.c
@@ -76,28 +76,19 @@ static signed int locate_peak(double x, double y, double z, double lambda,
 
 
 static double excitation_error(double xl, double yl, double zl,
-                               double ds_sq, double lambda)
+                               double ds, double lambda)
 {
-	double tt;
-	double a, b, c;
-	double r1, r2;
-	int n;
+	double tt, al;
+	double r;
 
 	tt = angle_between(0.0, 0.0, 1.0,  xl, yl, zl+1.0/lambda);
+	al = M_PI_2 - asin(-zl/ds);
 
-	a = 1.0;
-	b = - cos(tt) * 2.0/lambda;
-	c = pow(lambda, -2.0) - ds_sq;
+	r = ds * sin(al) / sin(tt);
 
-	/* FIXME: I don't trust GSL's quadratic solver */
-	n = gsl_poly_solve_quadratic(a, b, c, &r1, &r2);
-	assert(n == 2);
+	r -= 1.0/lambda;
 
-	r1 -= 1.0/lambda;
-	r2 -= 1.0/lambda;
-
-	if ( r1 > r2 ) return r2;
-	return r1;
+	return r;
 }
 
 
@@ -115,7 +106,7 @@ struct cpeak *find_intersections(struct image *image, UnitCell *cell,
 	double bandwidth = image->bw;
 	double divergence = image->div;
 	double lambda = image->lambda;
-	const double profile_cutoff = 0.05e9;  /* 0.1 nm^-1 */
+	const double profile_cutoff = 0.005e9;  /* 0.1 nm^-1 */
 	double llow, lhigh;    /* Wavelength */
 
 	cpeaks = malloc(sizeof(struct cpeak)*MAX_CPEAKS);
@@ -148,10 +139,10 @@ struct cpeak *find_intersections(struct image *image, UnitCell *cell,
 		double xl, yl, zl;
 		double ds, ds_sq;
 		double rlow, rhigh;    /* "Excitation error" */
-		signed int plow, phigh;
-		double xdalow, xdahigh;
-		double ydalow, ydahigh;
+		signed int p;
+		double lcen;
 		double xda, yda;
+		int close, straddled;
 
 		/* Ignore central beam */
 		if ( (h==0) && (k==0) && (l==0) ) continue;
@@ -168,32 +159,32 @@ struct cpeak *find_intersections(struct image *image, UnitCell *cell,
 		if ( ds > mres ) continue;  /* Outside resolution range */
 
 		/* Calculate excitation errors */
-		rlow = excitation_error(xl, yl, zl, ds_sq, llow);
-		rhigh = excitation_error(xl, yl, zl, ds_sq, lhigh);
-
-		if ( (fabs(rlow) > profile_cutoff)
-		  && (fabs(rhigh) > profile_cutoff) ) {
-			/* Reflection is not close to Bragg at either of
-			 * the two extremes of wavelength, so skip it. */
-			continue;
-		}
+		rlow = excitation_error(xl, yl, zl, ds, llow);
+		rhigh = excitation_error(xl, yl, zl, ds, lhigh);
+
+		/* 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"? */
+		straddled = signbit(rlow) ^ signbit(rhigh);
+
+		/* Neither?  Skip it. */
+		if ( !(close || straddled) ) continue;
+
+		lcen = -2.0*zl / ds_sq;
 
 		/* Locate peak on detector, and check it doesn't span panels */
-		plow = locate_peak(xl, yl, zl, llow, image->det,
-		                   &xdalow, &ydalow);
-		if ( plow == -1 ) continue;
-		phigh = locate_peak(xl, yl, zl, lhigh, image->det,
-		                    &xdahigh, &ydahigh);
-		if ( phigh == -1 ) continue;
-		if ( plow != phigh ) continue;
-
-		xda = xdahigh;
-		yda = ydahigh;
+		p = locate_peak(xl, yl, zl, lcen, image->det, &xda, &yda);
+		if ( p == -1 ) continue;
+
 		cpeaks[np].h = h;
 		cpeaks[np].k = k;
 		cpeaks[np].l = l;
-		cpeaks[np].x = xdahigh;
-		cpeaks[np].y = ydahigh;
+		cpeaks[np].x = xda;
+		cpeaks[np].y = yda;
 		np++;
 
 		if ( output ) {
diff --git a/src/templates.c b/src/templates.c
index e9bd7b72..4f47c572 100644
--- a/src/templates.c
+++ b/src/templates.c
@@ -173,8 +173,7 @@ IndexingPrivate *generate_templates(UnitCell *cell, const char *filename,
 
 		cell_rot = rotate_cell(cell, omega, phi, 0.0);
 
-		cpeaks = find_intersections(&image, cell_rot, 5.0e-3,
-		                          3.0/100.0, &n, 0);
+		cpeaks = find_intersections(&image, cell_rot, &n, 0);
 		if ( cpeaks == NULL ) {
 			ERROR("Template calculation failed.\n");
 			return NULL;