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

3 files changed, 49 insertions, 29 deletions

diff --git a/src/displaywindow.c b/src/displaywindow.c
index fa66382..1fef6ce 100644
--- a/src/displaywindow.c
+++ b/src/displaywindow.c
@@ -285,9 +285,11 @@ static gint displaywindow_gl_motion_notify(GtkWidget *widget, GdkEventMotion *ev
 
 static void displaywindow_gl_create_list(DisplayWindow *dw) {
 
+	GLfloat bblue[] = { 0.0, 0.0, 1.0, 1.0 };
 	GLfloat blue[] = { 0.0, 0.0, 0.5, 1.0 };
 	GLfloat blue_spec[] = { 0.0, 0.0, 1.0, 1.0 };
 	GLfloat red[] = { 1.0, 0.0, 0.0, 1.0 };
+	GLfloat purple[] = { 0.8, 0.0, 1.0, 1.0 };
 	GLfloat green[] = { 0.0, 1.0, 0.0, 1.0 };
 	GLfloat yellow[] = { 1.0, 1.0, 0.0, 1.0 };
 	GLfloat glass[] = { 0.2, 0.0, 0.8, 000.1 };
@@ -405,7 +407,8 @@ static void displaywindow_gl_create_list(DisplayWindow *dw) {
 	glNewList(dw->gl_list_id, GL_COMPILE);
 	
 	/* Bounding cube: 100 nm^-1 side length */
-	glMaterialfv(GL_FRONT_AND_BACK, GL_AMBIENT_AND_DIFFUSE, red);
+	glMaterialfv(GL_FRONT_AND_BACK, GL_AMBIENT_AND_DIFFUSE, purple);
+	glMaterialfv(GL_FRONT_AND_BACK, GL_EMISSION, black);
 	glMaterialfv(GL_FRONT_AND_BACK, GL_SPECULAR, black);
 	glMaterialf(GL_FRONT_AND_BACK, GL_SHININESS, 0.0);
 	glBegin(GL_LINE_LOOP);
@@ -430,10 +433,29 @@ static void displaywindow_gl_create_list(DisplayWindow *dw) {
 		glVertex3f(50, -50, 50);
 		glVertex3f(50, -50, -50);
 	glEnd();
-
+	
+	/* x, y, z pointers */
+	glBegin(GL_LINES);
+		glMaterialfv(GL_FRONT_AND_BACK, GL_AMBIENT_AND_DIFFUSE, black);
+		glMaterialfv(GL_FRONT_AND_BACK, GL_EMISSION, red);
+		glVertex3f(0.0, 0.0, 0.0);
+		glVertex3f(5.0, 0.0, 0.0);
+		glMaterialfv(GL_FRONT_AND_BACK, GL_EMISSION, green);
+		glVertex3f(0.0, 0.0, 0.0);
+		glVertex3f(0.0, 5.0, 0.0);
+		glMaterialfv(GL_FRONT_AND_BACK, GL_EMISSION, bblue);
+		glVertex3f(0.0, 0.0, 0.0);
+		glVertex3f(0.0, 0.0, 5.0);
+		glMaterialfv(GL_FRONT_AND_BACK, GL_EMISSION, black);
+	glEnd();
+	
 	/* Tilt axis */
 	glPushMatrix();
-	glRotatef(-ctx->images[0].omega, 0.0, 0.0, 1.0);
+	/* Images rotate clockwise by omega to put tilt axis at +x,
+	 *	so rotate tilt axis anticlockwise by omega.
+	 *	Since the rotation is about +z, this is already anticlockwise
+	 *	when looking down z. */
+	glRotatef(ctx->images[0].omega, 0.0, 0.0, 1.0);
 	glPolygonOffset(1.0, 1.0);
 	glEnable(GL_POLYGON_OFFSET_LINE);
 	glBegin(GL_LINES);
@@ -581,7 +603,7 @@ static void displaywindow_gl_create_list(DisplayWindow *dw) {
 	
 	glEndList();
 	
-	printf("DW: Vertex counts: meas:%i, mark:%i, gen:%i\n", dw->gl_ref_num_vertices, dw->gl_marker_num_vertices, dw->gl_gen_num_vertices);
+	//printf("DW: Vertex counts: meas:%i, mark:%i, gen:%i\n", dw->gl_ref_num_vertices, dw->gl_marker_num_vertices, dw->gl_gen_num_vertices);
 
 }
 
diff --git a/src/reflections.c b/src/reflections.c
index a490364..dde8d84 100644
--- a/src/reflections.c
+++ b/src/reflections.c
@@ -162,26 +162,26 @@ void reflection_add_from_dp(ControlContext *ctx, double x, double y, ImageRecord
 	psi = atan2(y, x);
 	
 	x_temp = k*sin(theta)*cos(psi);
-	y_temp = -k*sin(theta)*sin(psi);	/* Minus sign to define axes as y going upwards */
-	z_temp = k- k*cos(theta);
+	y_temp = k*sin(theta)*sin(psi);
+	z_temp = k - k*cos(theta);
 	
 	/* Apply the rotations...
 		First: rotate image clockwise until tilt axis is aligned horizontally. */
-	nx = x_temp*cos(omega) + y_temp*-sin(omega);
-	ny = x_temp*sin(omega) + y_temp*cos(omega);
+	nx = x_temp*cos(omega) + y_temp*sin(omega);
+	ny = -x_temp*sin(omega) + y_temp*cos(omega);
 	nz = z_temp;
 	
 	/* Now, tilt about the x-axis ANTICLOCKWISE around +x, i.e. the "wrong" way.
 		This is because the crystal is rotated in the experiment, not the Ewald sphere. */
 	x_temp = nx; y_temp = ny; z_temp = nz;	
 	nx = x_temp;
-	ny = cos(tilt)*y_temp - sin(tilt)*z_temp;
-	nz = -sin(tilt)*y_temp - cos(tilt)*z_temp;
+	ny = cos(tilt)*y_temp + sin(tilt)*z_temp;
+	nz = -sin(tilt)*y_temp + cos(tilt)*z_temp;
 	
 	/* Finally, reverse the omega rotation to restore the location of the image in 3D space */
 	x_temp = nx; y_temp = ny; z_temp = nz;	
-	nx = x_temp*cos(-omega) + y_temp*-sin(-omega);
-	ny = x_temp*sin(-omega) + y_temp*cos(-omega);
+	nx = x_temp*cos(-omega) + y_temp*sin(-omega);
+	ny = -x_temp*sin(-omega) + y_temp*cos(-omega);
 	nz = z_temp;
 	
 	reflection_add(ctx->reflectionctx, nx, ny, nz, intensity, REFLECTION_NORMAL);
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);