diff options
author | taw27 <taw27@bf6ca9ba-c028-0410-8290-897cf20841d1> | 2007-10-31 16:56:05 +0000 |
---|---|---|
committer | taw27 <taw27@bf6ca9ba-c028-0410-8290-897cf20841d1> | 2007-10-31 16:56:05 +0000 |
commit | 68aaef4f7d3a73153f6052b731aae559da45a4d2 (patch) | |
tree | 01c86af1cb8a8e8462cd9f544c901f11170ca094 | |
parent | 4676d2c74e1a29d8aa5c0da56cb67ba7a6bb7e0f (diff) |
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
-rw-r--r-- | src/refine.c | 104 |
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 */ |