Merge branch 'master' into simple-search

2 files changed, 17 insertions, 16 deletions

diff --git a/src/basis.c b/src/basis.c
index 37fedba..2736d0e 100644
--- a/src/basis.c
+++ b/src/basis.c
@@ -22,6 +22,7 @@
 #include "basis.h"
 #include "utils.h"
 #include "displaywindow.h"
+#include "reproject.h"
 
 static void basis_print(Basis *cell)
 {
@@ -224,6 +225,7 @@ static gint basis_load_response(GtkWidget *widget, gint response,
 						ctx->dw);
 		} else {
 			displaywindow_update(ctx->dw);
+			reproject_lattice_changed(ctx);
 		}
 		g_free(filename);
 	}
diff --git a/src/reproject.c b/src/reproject.c
index c41aa3f..df85f4f 100644
--- a/src/reproject.c
+++ b/src/reproject.c
@@ -49,9 +49,9 @@ ImageFeatureList *reproject_get_reflections(ImageRecord *image, ReflectionList *
 
 	ImageFeatureList *flist;
 	Reflection *reflection;
-	double smax = 0.1e9;
-	double tilt, omega;
-	double nx, ny, nz, nxt, nyt, nzt;	/* "normal" vector (and calculation intermediates) */
+	double smax = 0.4e9;
+	double tilt, omega, k;
+	double nx, ny, nz;			/* "normal" vector */
 	double kx, ky, kz;			/* Electron wavevector ("normal" times 1/lambda */
 	double ux, uy, uz, uxt, uyt, uzt;	/* "up" vector (and calculation intermediates) */
 	//int done_debug = 0;
@@ -60,13 +60,10 @@ ImageFeatureList *reproject_get_reflections(ImageRecord *image, ReflectionList *
 
 	tilt = image->tilt;
 	omega = image->omega;
+	k = 1.0/image->lambda;
 
-	/* Calculate the (normalised) incident electron wavevector
-	 *	n is rotated "into" the reconstruction, so only one omega step. */
-	nxt = 0.0;  nyt = 0.0;  nzt = 1.0;
-	nx = nxt;  ny = cos(tilt)*nyt + sin(tilt)*nzt;  nz = -sin(tilt)*nyt + cos(tilt)*nzt;
-	nxt = nx;  nyt = ny;  nzt = nz;
-	nx = nxt*cos(-omega) + nyt*sin(-omega);  ny = -nxt*sin(-omega) + nyt*cos(-omega);  nz = nzt;
+	/* Calculate the (normalised) incident electron wavevector */
+	mapping_rotate(0.0, 0.0, 1.0, &nx, &ny, &nz, omega, tilt);
 	kx = nx / image->lambda;
 	ky = ny / image->lambda;
 	kz = nz / image->lambda;	/* This is the centre of the Ewald sphere */
@@ -85,21 +82,23 @@ ImageFeatureList *reproject_get_reflections(ImageRecord *image, ReflectionList *
 
 		double xl, yl, zl;
 		double a, b, c;
-		double A1, A2, s1, s2, s;
+		double s1, s2, s, t;
+		double g_sq, gn;
 
 		/* Get the coordinates of the reciprocal lattice point */
 		xl = reflection->x;
 		yl = reflection->y;
 		zl = reflection->z;
+		g_sq = modulus_squared(xl, yl, zl);
+		gn = xl*nx + yl*ny + zl*nz;
 
 		/* Next, solve the relrod equation to calculate the excitation error */
 		a = 1.0;
-		b = -2.0*(xl*nx + yl*ny + zl*nz);
-		c = xl*xl + yl*yl + zl*zl - 1.0/(image->lambda*image->lambda);
-		A1 = (-b + sqrt(b*b-4.0*a*c))/(2.0*a);
-		A2 = (-b - sqrt(b*b-4.0*a*c))/(2.0*a);
-		s1 = 1.0/image->lambda - A1;
-		s2 = 1.0/image->lambda - A2;
+		b = 2.0*(gn - k);
+		c = -2.0*gn*k + g_sq;
+		t = -0.5*(b + sign(b)*sqrt(b*b - 4.0*a*c));
+		s1 = t/a;
+		s2 = c/t;
 		if ( fabs(s1) < fabs(s2) ) s = s1; else s = s2;
 
 		/* Skip this reflection if s is large */