Vaguely working refinement to work with

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

1 files changed, 94 insertions, 10 deletions

diff --git a/src/refine.c b/src/refine.c
index aaefdf9..c996691 100644
--- a/src/refine.c
+++ b/src/refine.c
@@ -17,6 +17,7 @@
 #include <math.h>
 #include <stdlib.h>
 #include <stdio.h>
+#include <assert.h>
 
 #include "displaywindow.h"
 #include "gtk-valuegraph.h"
@@ -56,6 +57,12 @@ static double refine_image_deviation(ImageRecord *image, ReflectionList *cell_la
 
 }
 
+typedef enum {
+	INDEX_A = 1<<0,
+	INDEX_B = 1<<1,
+	INDEX_C = 1<<2
+} RefinementIndex;
+
 /* Use the IPR algorithm to make "cell" fit the given image */
 static void refine_fit_image(Basis *cell, ImageRecord *image) {
 
@@ -63,6 +70,14 @@ static void refine_fit_image(Basis *cell, ImageRecord *image) {
 	ImageFeatureList *flist;
 	int i;
 	ReflectionList *cell_lattice;
+	Basis cd; /* Cell delta */
+	int n_a = 0;
+	int n_b = 0;
+	int n_c = 0;
+	
+	cd.a.x = 0.0;	cd.a.y = 0.0;	cd.a.z = 0.0;
+	cd.b.x = 0.0;	cd.b.y = 0.0;	cd.b.z = 0.0;
+	cd.c.x = 0.0;	cd.c.y = 0.0;	cd.c.z = 0.0;
 	
 	cell_lattice = reflection_list_from_cell(cell);
 	rflist = reproject_get_reflections(image, cell_lattice);
@@ -71,36 +86,105 @@ static void refine_fit_image(Basis *cell, ImageRecord *image) {
 	
 	for ( i=0; i<rflist->n_features; i++ ) {
 		
-		double dx, dy;
-		double x, y, z, twotheta;
+		double dix, diy;
+		double dx, dy, dz, twotheta;
 		double old_x, old_y;
+		RefinementIndex index, delta, shared;
+		signed int h, k, l;
+		double a11, a12, a13, a21, a22, a23, a31, a32, a33, det;
+		double dh, dk, dl;
 		
 		/* Skip if no partner */
 		if ( !rflist->features[i].partner ) continue;
 		
 		/* Determine the difference vector */
-		dx = rflist->features[i].partner->x - rflist->features[i].x;
-		dy = rflist->features[i].partner->y - rflist->features[i].y;
+		dix = rflist->features[i].partner->x - rflist->features[i].x;
+		diy = rflist->features[i].partner->y - rflist->features[i].y;
 		
 		/* Map the difference vector to the relevant tilted plane */
 		old_x = rflist->features[i].partner->x;
 		old_y = rflist->features[i].partner->y;
-		rflist->features[i].partner->x = dx;
-		rflist->features[i].partner->y = dy;
-		mapping_map_to_space(rflist->features[i].partner, &x, &y, &z, &twotheta);
+		rflist->features[i].partner->x = dix;
+		rflist->features[i].partner->y = diy;
+		mapping_map_to_space(rflist->features[i].partner, &dx, &dy, &dz, &twotheta);
 		rflist->features[i].partner->x = old_x;
 		rflist->features[i].partner->y = old_y;
 		
-		printf("Feature %3i: %3i %3i %3i dev=%8e %8e %8e (%5f mrad)\n",
-			i, rflist->features[i].h, rflist->features[i].k, rflist->features[i].l, x, y, z, twotheta*1e3);
+		h = rflist->features[i].h;
+		k = rflist->features[i].k;
+		l = rflist->features[i].l;
+		printf("Feature %3i: %3i %3i %3i dev=%8e %8e %8e (%5f mrad)\n", i, h, k, l, dx, dy, dz, twotheta*1e3);
+		
+		/* Select the basis vectors which are allowed to be altered */
+		index = 0;
+		if ( h ) index |= INDEX_A;
+		if ( k ) index |= INDEX_B;
+		if ( l ) index |= INDEX_C;
+		assert(index != 0);	/* Can't refine using the central beam! */
 		
-		/* Work out weightings for a, b and c */
+		/* Set up the coordinate transform from hkl to xyz */
+		a11 = cell->a.x;  a12 = cell->a.y;  a13 = cell->a.z;
+		a21 = cell->b.x;  a22 = cell->b.y;  a23 = cell->b.z;
+		a31 = cell->c.x;  a32 = cell->c.y;  a33 = cell->c.z;
 		
+		/* Invert the matrix to get dh,dk,dl from dx,dy,dz */
+		det = a11*(a22*a33 - a23*a32) - a12*(a21*a33 - a23*a31) + a13*(a21*a32 - a22*a31);
+		dh = ((a22*a33-a23*a32)*dx + (a23*a31-a21*a33)*dy + (a21*a32-a22*a31)*dz) / det;
+		dk = ((a13*a32-a12*a33)*dx + (a11*a33-a13*a31)*dy + (a12*a31-a11*a32)*dz) / det;
+		dl = ((a12*a23-a13*a22)*dx + (a13*a21-a11*a23)*dy + (a11*a22-a12*a21)*dz) / det;
+		printf("dev(hkl) = %f %f %f\n", dh, dk, dl);
 	
+		delta = 0;
+		if ( fabs(dh) < 0.001 ) delta |= INDEX_A;
+		if ( fabs(dk) < 0.001 ) delta |= INDEX_B;
+		if ( fabs(dl) < 0.001 ) delta |= INDEX_C;
+		/* Typically, 'delta' will have all three bits set */
+		
+		shared = index & delta;
+		
+		if ( shared == 0 ) {
+			/* No indices left - 'pure shear' (delta is perpendicular (in the abc basis) to index) */
+			shared = index;
+		}
+		
+		if ( shared & INDEX_A ) {
+			double w = h / (h+k+l);
+			cd.a.x += w*dx / h;
+			cd.a.y += w*dy / h;
+			cd.a.z += w*dz / h;
+			n_a++;
+		}
+		if ( shared & INDEX_B ) {
+			double w = k / (h+k+l);
+			cd.b.x += w*dx / k;
+			cd.b.y += w*dy / k;
+			cd.b.z += w*dz / k;
+			n_b++;
+		}
+		if ( shared & INDEX_C ) {
+			double w = l / (h+k+l);
+			cd.c.x += w*dx / l;
+			cd.c.y += w*dy / l;
+			cd.c.z += w*dz / l;
+			n_c++;
+		}
+
 	}
 	
 	image_feature_list_free(rflist);
 	
+	cd.a.x /= n_a;	cd.a.y /= n_a;	cd.a.z /= n_a;
+	cd.b.x /= n_b;	cd.b.y /= n_b;	cd.b.z /= n_b;
+	cd.c.x /= n_c;	cd.c.y /= n_c;	cd.c.z /= n_c;
+	
+	printf("Total distortion(a) = %+8e %+8e %+8e\n", cd.a.x, cd.a.y, cd.a.z);
+	printf("Total distortion(b) = %+8e %+8e %+8e\n", cd.b.x, cd.b.y, cd.b.z);
+	printf("Total distortion(c) = %+8e %+8e %+8e\n", cd.c.x, cd.c.y, cd.c.z);
+	
+	cell->a.x += cd.a.x;	cell->a.y += cd.a.y;	cell->a.z += cd.a.z;
+	cell->b.x += cd.b.x;	cell->b.y += cd.b.y;	cell->b.z += cd.b.z;
+	cell->c.x += cd.c.x;	cell->c.y += cd.c.y;	cell->c.z += cd.c.z;
+	
 }
 
 /* Display a graph of root sum squared deviation distance against some other parameter */