theta/psi formulation for reprojection

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

1 files changed, 34 insertions, 17 deletions

diff --git a/src/reproject.c b/src/reproject.c
index bec9a31..1deb044 100644
--- a/src/reproject.c
+++ b/src/reproject.c
@@ -14,6 +14,7 @@
 
 #include "control.h"
 #include "reflections.h"
+#include "utils.h"
 
 #define MAX_IMAGE_REFLECTIONS 8*1024
 
@@ -65,36 +66,52 @@ ImageReflection *reproject_get_reflections(ImageRecord image, size_t *n, Reflect
 		/* Skip this reflection if s is large */
 		if ( fabs(s) <= smax ) {
 		
-			double xddd, yddd, zddd;
-			double xdd, ydd, zdd;
-			double xd, yd, zd;
+			double xi, yi, zi;
+			double gx, gy, gz;
+			double kx, ky, kz;
+			double cx, cy, cz;
+			double ux, uy, uz;
+			double uxt, uyt, uzt;
 			double theta, psi;
 			double x, y;
 			
 			/* Determine the intersection point */
-			xddd = xl + s*nx;  yddd = yl + s*ny;  zddd = zl + s*nz;
+			xi = xl + s*nx;  yi = yl + s*ny;  zi = zl + s*nz;
 			reflection_add(ctx->reflectionctx, xl, yl, zl, 1, REFLECTION_CENTRAL);
-			reflection_add(ctx->reflectionctx, xddd, yddd, zddd, 1, REFLECTION_MARKER);
+			reflection_add(ctx->reflectionctx, xi, yi, zi, 1, REFLECTION_MARKER);
 			
-			/* Invert the image->3D mapping to get the image coordinates */
-			xdd = xddd;
-			ydd = (yddd/cos(tilt) - zddd*tan(tilt)/cos(tilt))/(1+tan(tilt)*tan(tilt));
-			zdd = (-zddd-ydd*sin(tilt))/cos(tilt);
-			
-			yd = (ydd-xdd*tan(omega))/(sin(omega)*tan(omega)+cos(omega));	/* Needs to be done first */
-			xd = (xdd+yd*sin(omega))/cos(omega);
-			zd = zdd;
+			/* Calculate Bragg angle and azimuth of point in image */
+			kx = nx / image.lambda;
+			ky = ny / image.lambda;
+			kz = nz / image.lambda;	/* This is the centre of the Ewald sphere */
+			gx = xi - kx;
+			gy = yi - ky;
+			gz = zi - kz;	/* This is the vector from the centre of the sphere to the intersection */
+			theta = angle_between(-kx, -ky, -kz, gx, gy, gz);
+			printf("theta=%f   ", theta);
+			theta = deg2rad(theta);
+			cx = nz*gy - ny*gz;
+			cy = nz*gx - nx*gz;
+			cz = ny*gx - nx*gy;
+			uxt = 0;
+			uyt = cos(tilt);
+			uzt = -sin(tilt);
+			ux = uyt*-sin(omega);
+			uy = uyt*cos(omega);
+			uz = uzt;
+			psi = angle_between(cx, cy, cz, ux, uy, uz);
+			psi -= 90;
+			printf("psi=%f\n", psi);
+			psi = deg2rad(psi);
 			
 			if ( image.fmode == FORMULATION_CLEN ) {
-				psi = atan2(-yd, xd);
-				theta = acos(1+zd*image.lambda);
 				x = image.camera_len*sin(theta)*cos(psi);
 				y = image.camera_len*sin(theta)*sin(psi);
 				x *= image.resolution;
 				y *= image.resolution;
 			} else if ( image.fmode == FORMULATION_PIXELSIZE ) {
-				x = xd / image.pixel_size;
-				y = yd / image.pixel_size;
+				x = 0;
+				y = 0;
 			} else {
 				fprintf(stderr, "Unrecognised formulation mode in reproject_get_reflections()\n");
 				return NULL;