Change the formulation of the rotations

git-svn-id: svn://cook.msm.cam.ac.uk:745/diff-tomo/dtr@121 bf6ca9ba-c028-0410-8290-897cf20841d1

1 files changed, 15 insertions, 17 deletions

diff --git a/src/reproject.c b/src/reproject.c
index 4f0dd3c..9a54281 100644
--- a/src/reproject.c
+++ b/src/reproject.c
@@ -25,9 +25,9 @@ ImageReflection *reproject_get_reflections(ImageRecord image, size_t *n, Reflect
 	int i;
 	double smax = 0.5e9;
 	double tilt, omega;
-	double xt, yt, zt;
-	double nx, ny, nz;
-	double kx, ky, kz;
+	double nx, ny, nz, nxt, nyt, nzt;	/* "normal" vector (and calculation intermediates) */
+	double kx, ky, kz;			/* Electron wavevector ("normal" times 1/lambda */
+	double ux, uy, uz, uxt, uyt, uzt;	/* "up" vector (and calculation intermediates) */
 			
 	tilt = deg2rad(image.tilt);
 	omega = deg2rad(image.omega);
@@ -39,16 +39,22 @@ ImageReflection *reproject_get_reflections(ImageRecord image, size_t *n, Reflect
 	reflection = rctx->reflections;
 	
 	/* Calculate the (normalised) incident electron wavevector */
-	xt = 0;  yt = sin(tilt);  zt = cos(tilt);
-	nx = xt*cos(omega) + yt*-sin(omega);
-	ny = xt*sin(omega) + yt*cos(omega);
-	nz = zt;
+	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;
 	kx = nx / image.lambda;
 	ky = ny / image.lambda;
 	kz = nz / image.lambda;	/* This is the centre of the Ewald sphere */
-	
 	reflection_add(ctx->reflectionctx, kx, ky, kz, 1, REFLECTION_VECTOR_MARKER_1);
 	
+	/* Determine where "up" is */
+	uxt = sin(-omega);  uyt = cos(-omega);  uzt = 0;
+	ux = uxt;  uy = cos(tilt)*uyt + sin(tilt)*uzt;  uz = -sin(tilt)*uyt + cos(tilt)*uzt;
+	uxt = ux;  uyt = uy;  uzt = uz;
+	ux = uxt*cos(omega) + uyt*-sin(omega);  uy = -uxt*sin(omega) + uyt*cos(omega);  uz = uzt;
+	reflection_add(ctx->reflectionctx, ux*50, uy*50, uz*50, 1, REFLECTION_VECTOR_MARKER_2);
+	
 	do {
 	
 		double xl, yl, zl;
@@ -77,8 +83,6 @@ ImageReflection *reproject_get_reflections(ImageRecord image, size_t *n, Reflect
 			double xi, yi, zi;
 			double gx, gy, gz;
 			double cx, cy, cz;
-			double ux, uy, uz;
-			double uyt, uzt;
 			double theta;
 			double x, y;
 			double rx, ry, rz;
@@ -95,16 +99,10 @@ ImageReflection *reproject_get_reflections(ImageRecord image, size_t *n, Reflect
 			if ( theta > 0.0 ) {
 			
 				double psi, disc;
-			
+				
 				reflection_add(ctx->reflectionctx, xl, yl, zl, 1, REFLECTION_GENERATED);
 				reflection_add(ctx->reflectionctx, xi, yi, zi, 1, REFLECTION_MARKER);
 				
-				/* Determine where "up" is */
-				uyt = cos(tilt);	uzt = -sin(tilt);  /* tilt it (uxt not needed) */
-				ux = uyt*-sin(omega);	uy = uyt*cos(omega);	uz = uzt;  /* rotate it */
-				
-				if ( i == 0 ) reflection_add(ctx->reflectionctx, ux*50, uy*50, uz*50, 1, REFLECTION_VECTOR_MARKER_2);
-				
 				/* Calculate azimuth of point in image (clockwise from "up", will be changed later) */
 				cx = nz*yi - ny*zi;  cy = nx*zi - nz*xi;  cz = ny*xi - nx*yi;  /* c = (-n) x i */
 				psi = angle_between(cx, cy, cz, ux, uy, uz);