Vaguely working LSQ lattice refinement (doesn't quite work)

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

1 files changed, 59 insertions, 23 deletions

diff --git a/src/mapping.c b/src/mapping.c
index f9b47c0..e1e59c5 100644
--- a/src/mapping.c
+++ b/src/mapping.c
@@ -20,6 +20,35 @@
 #include "displaywindow.h"
 #include "cache.h"
 
+void mapping_rotate(double x, double y, double z, double *ddx, double *ddy, double *ddz, double omega, double tilt) {
+
+	double nx, ny, nz;
+	double x_temp, y_temp, z_temp;
+	
+	/* First: rotate image clockwise until tilt axis is aligned horizontally. */
+	nx = x*cos(omega) + y*sin(omega);
+	ny = -x*sin(omega) + y*cos(omega);
+	nz = z;
+	
+	/* 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;
+	
+	/* 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);
+	nz = z_temp;
+	
+	*ddx = nx;
+	*ddy = ny;
+	*ddz = nz;
+	
+}
+
 int mapping_map_to_space(ImageFeature *refl, double *ddx, double *ddy, double *ddz, double *twotheta) {
 
 	/* "Input" space */
@@ -34,7 +63,6 @@ int mapping_map_to_space(ImageFeature *refl, double *ddx, double *ddy, double *d
 	double omega;
 	
 	double x_temp, y_temp, z_temp;
-	double nx, ny, nz;
 	
 	imagerecord = refl->parent;
 	x = refl->x - imagerecord->x_centre;	y = refl->y - imagerecord->y_centre;
@@ -54,7 +82,7 @@ int mapping_map_to_space(ImageFeature *refl, double *ddx, double *ddy, double *d
 		d = sqrt((x*x) + (y*y));
 		theta = atan2(d, k);
 	} else {
-		fprintf(stderr, "Unrecognised formulation mode in reflection_add_from_dp\n");
+		fprintf(stderr, "Unrecognised formulation mode in mapping_map_to_space.\n");
 		return -1;
 	}
 	psi = atan2(y, x);
@@ -63,34 +91,42 @@ int mapping_map_to_space(ImageFeature *refl, double *ddx, double *ddy, double *d
 	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);
-	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;
-	
-	/* 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);
-	nz = z_temp;
+	mapping_rotate(x_temp, y_temp, z_temp, ddx, ddy, ddz, omega, tilt);
 	
-	*ddx = nx;
-	*ddy = ny;
-	*ddz = nz;
 	*twotheta = theta;	/* Radians.  I've used the "wrong" nomenclature above */
 	
 	return 0;
 
 }
 
+int mapping_scale(ImageFeature *refl, double *ddx, double *ddy) {
+
+	double x, y;
+	ImageRecord *imagerecord;
+	double k;
+	
+	imagerecord = refl->parent;
+	x = refl->x - imagerecord->x_centre;
+	y = refl->y - imagerecord->y_centre;
+	k = 1/imagerecord->lambda;
+	
+	if ( imagerecord->fmode == FORMULATION_CLEN ) {
+		x /= imagerecord->resolution;
+		y /= imagerecord->resolution;	/* Convert pixels to metres */
+		*ddx = x * k / imagerecord->camera_len;
+		*ddy = y * k / imagerecord->camera_len;	
+	} else if (imagerecord->fmode == FORMULATION_PIXELSIZE ) {
+		*ddx = x * imagerecord->pixel_size;
+		*ddy = y * imagerecord->pixel_size;	/* Convert pixels to metres^-1 */
+	} else {
+		fprintf(stderr, "Unrecognised formulation mode in mapping_scale.\n");
+		return -1;
+	}
+	
+	return 0;
+	
+}
+
 void mapping_map_features(ControlContext *ctx) {
 
 	int i;