Match the unit cell to a model cell after indexing

10 files changed, 259 insertions, 35 deletions

diff --git a/src/cell.c b/src/cell.c
index 86be7cc4..a3fdc72a 100644
--- a/src/cell.c
+++ b/src/cell.c
@@ -22,6 +22,7 @@
 
 #include "cell.h"
 #include "utils.h"
+#include "image.h"
 
 
 /* Update the cartesian representation from the crystallographic one */
@@ -190,6 +191,55 @@ UnitCell *cell_new_from_parameters(double a, double b, double c,
 }
 
 
+static UnitCell *cell_new_from_axes(struct rvec as, struct rvec bs,
+                                    struct rvec cs)
+{
+	UnitCell *cell;
+	int s;
+	gsl_matrix *m;
+	gsl_matrix *inv;
+	gsl_permutation *perm;
+
+	cell = cell_new();
+	if ( !cell ) return NULL;
+
+	m = gsl_matrix_alloc(3, 3);
+	gsl_matrix_set(m, 0, 0, as.u);
+	gsl_matrix_set(m, 0, 1, as.v);
+	gsl_matrix_set(m, 0, 2, as.w);
+	gsl_matrix_set(m, 1, 0, bs.u);
+	gsl_matrix_set(m, 1, 1, bs.v);
+	gsl_matrix_set(m, 1, 2, bs.w);
+	gsl_matrix_set(m, 2, 0, cs.u);
+	gsl_matrix_set(m, 2, 1, cs.v);
+	gsl_matrix_set(m, 2, 2, cs.w);
+
+	/* Invert */
+	perm = gsl_permutation_alloc(m->size1);
+	inv = gsl_matrix_alloc(m->size1, m->size2);
+	gsl_linalg_LU_decomp(m, perm, &s);
+	gsl_linalg_LU_invert(m, perm, inv);
+	gsl_permutation_free(perm);
+	gsl_matrix_free(m);
+
+	/* Transpose */
+	gsl_matrix_transpose(inv);
+
+	cell->ax = gsl_matrix_get(inv, 0, 0);
+	cell->ay = gsl_matrix_get(inv, 0, 1);
+	cell->az = gsl_matrix_get(inv, 0, 2);
+	cell->bx = gsl_matrix_get(inv, 1, 0);
+	cell->by = gsl_matrix_get(inv, 1, 1);
+	cell->bz = gsl_matrix_get(inv, 1, 2);
+	cell->cx = gsl_matrix_get(inv, 2, 0);
+	cell->cy = gsl_matrix_get(inv, 2, 1);
+	cell->cz = gsl_matrix_get(inv, 2, 2);
+
+	cell_update_crystallographic(cell);
+	return cell;
+}
+
+
 void cell_get_cartesian(UnitCell *cell,
                         double *ax, double *ay, double *az,
                         double *bx, double *by, double *bz,
@@ -247,6 +297,187 @@ void cell_get_reciprocal(UnitCell *cell,
 }
 
 
+static void cell_print(UnitCell *cell)
+{
+	STATUS("  a     b     c         alpha   beta  gamma\n");
+	STATUS("%5.2f %5.2f %5.2f nm    %6.2f %6.2f %6.2f deg\n",
+	       cell->a*1e9, cell->b*1e9, cell->c*1e9,
+	       rad2deg(cell->alpha), rad2deg(cell->beta), rad2deg(cell->gamma));
+}
+
+
+#define MAX_CAND (1024)
+
+static int within_tolerance(double a, double b, double percent)
+{
+	double tol = a * (percent/100.0);
+	if ( fabs(b-a) < tol ) return 1;
+	return 0;
+}
+
+
+struct cvec {
+	struct rvec vec;
+	float na;
+	float nb;
+	float nc;
+};
+
+
+/* Attempt to make 'cell' fit into 'template' somehow */
+UnitCell *match_cell(UnitCell *cell, UnitCell *template)
+{
+	signed int n1l, n2l, n3l;
+	double asx, asy, asz;
+	double bsx, bsy, bsz;
+	double csx, csy, csz;
+	int i, j;
+	double lengths[3];
+	double angles[3];
+	struct cvec *cand[3];
+
+	int ncand[3] = {0,0,0};
+	float ltl = 5.0;
+	float angtol = 5.0;
+	UnitCell *new_cell = NULL;
+
+	cell_get_reciprocal(template, &asx, &asy, &asz,
+	                              &bsx, &bsy, &bsz,
+	                              &csx, &csy, &csz);
+
+	lengths[0] = modulus(asx, asy, asz);
+	lengths[1] = modulus(bsx, bsy, bsz);
+	lengths[2] = modulus(csx, csy, csz);
+
+	angles[0] = angle_between(bsx, bsy, bsz, csx, csy, csz);
+	angles[1] = angle_between(asx, asy, asz, csx, csy, csz);
+	angles[2] = angle_between(asx, asy, asz, bsx, bsy, bsz);
+
+	cand[0] = malloc(MAX_CAND*sizeof(struct cvec));
+	cand[1] = malloc(MAX_CAND*sizeof(struct cvec));
+	cand[2] = malloc(MAX_CAND*sizeof(struct cvec));
+
+	cell_get_reciprocal(cell, &asx, &asy, &asz,
+	                          &bsx, &bsy, &bsz,
+	                          &csx, &csy, &csz);
+
+	STATUS("Performing unit cell magic...\n");
+
+	/* Negative values mean 1/n, positive means n, zero means zero */
+	for ( n1l=-4; n1l<=4; n1l++ ) {
+	for ( n2l=-4; n2l<=4; n2l++ ) {
+	for ( n3l=-4; n3l<=4; n3l++ ) {
+
+		float n1, n2, n3;
+		signed int b1, b2, b3;
+
+		n1 = (n1l>=0) ? (n1l) : (1.0/n1l);
+		n2 = (n2l>=0) ? (n2l) : (1.0/n2l);
+		n3 = (n3l>=0) ? (n3l) : (1.0/n3l);
+
+		/* 'bit' values can be +1 or -1 */
+		for ( b1=-1; b1<=1; b1+=2 ) {
+		for ( b2=-1; b2<=1; b2+=2 ) {
+		for ( b3=-1; b3<=1; b3+=2 ) {
+
+			double tx, ty, tz;
+			double tlen;
+			int i;
+
+			n1 *= b1;  n2 *= b2;  n3 *= b3;
+
+			tx = n1*asx + n2*asy + n3*asz;
+			ty = n1*bsx + n2*bsy + n3*bsz;
+			tz = n1*csx + n2*csy + n3*csz;
+			tlen = modulus(tx, ty, tz);
+
+			/* Test modulus for agreement with moduli of template */
+			for ( i=0; i<3; i++ ) {
+				if ( within_tolerance(lengths[i], tlen, ltl) ) {
+					STATUS("sought %e, found %e (%e %e) for %i\n",
+					       lengths[i], tlen, 1.0/lengths[i],
+					       1.0/tlen, i);
+					cand[i][ncand[i]].vec.u = tx;
+					cand[i][ncand[i]].vec.v = ty;
+					cand[i][ncand[i]].vec.w = tz;
+					cand[i][ncand[i]].na = n1;
+					cand[i][ncand[i]].nb = n2;
+					cand[i][ncand[i]].nc = n3;
+					ncand[i]++;
+				}
+			}
+
+		}
+		}
+		}
+
+	}
+	}
+	}
+
+	for ( i=0; i<ncand[0]; i++ ) {
+	for ( j=0; j<ncand[1]; j++ ) {
+
+		double ang;
+		int k;
+
+		/* Measure the angle between the ith candidate for axis 0
+		 * and the jth candidate for axis 1 */
+		ang = angle_between(cand[0][i].vec.u, cand[0][i].vec.v,
+		                    cand[0][i].vec.w, cand[1][j].vec.u,
+		                    cand[1][j].vec.v, cand[1][j].vec.w);
+
+		/* Angle between axes 0 and 1 should be angle 2 */
+		if ( !within_tolerance(ang, angles[2], angtol) ) continue;
+
+		for ( k=0; k<ncand[2]; k++ ) {
+
+			/* Measure the angle between the current candidate for
+			 * axis 0 and the kth candidate for axis 2 */
+			ang = angle_between(cand[0][i].vec.u, cand[0][i].vec.v,
+			                    cand[0][i].vec.w, cand[2][k].vec.u,
+			                    cand[2][k].vec.v, cand[2][k].vec.w);
+
+			/* ... it should be angle 1 ... */
+			if ( !within_tolerance(ang, angles[1], angtol) ) continue;
+
+			/* Finally, the angle between the current candidate for
+			 * axis 1 and the kth candidate for axis 2 */
+			ang = angle_between(cand[1][j].vec.u, cand[1][j].vec.v,
+			                    cand[1][j].vec.w, cand[2][k].vec.u,
+			                    cand[2][k].vec.v, cand[2][k].vec.w);
+
+			/* ... it should be angle 0 ... */
+			if ( !within_tolerance(ang, angles[0], angtol) ) continue;
+
+			new_cell = cell_new_from_axes(cand[0][i].vec,
+			                              cand[1][j].vec,
+			                              cand[2][k].vec);
+
+			STATUS("Success! --------------- \n");
+			cell_print(new_cell);
+			STATUS("%f %f %f\n", cand[0][i].na, cand[0][i].nb,
+			                      cand[0][i].nc);
+			STATUS("%f %f %f\n", cand[1][j].na, cand[1][j].nb,
+			                      cand[1][j].nc);
+			STATUS("%f %f %f\n", cand[2][k].na, cand[2][k].nb,
+			                      cand[2][k].nc);
+			goto done;
+
+		}
+
+	}
+	}
+
+done:
+	free(cand[0]);
+	free(cand[1]);
+	free(cand[2]);
+
+	return new_cell;
+}
+
+
 double resolution(UnitCell *cell, signed int h, signed int k, signed int l)
 {
 	const double a = cell->a;
diff --git a/src/cell.h b/src/cell.h
index d1b49d22..687c7ef6 100644
--- a/src/cell.h
+++ b/src/cell.h
@@ -72,4 +72,6 @@ extern void cell_get_reciprocal(UnitCell *cell,
 extern double resolution(UnitCell *cell,
                          signed int h, signed int k, signed int l);
 
+extern UnitCell *match_cell(UnitCell *cell, UnitCell *template);
+
 #endif	/* CELL_H */
diff --git a/src/diffraction.c b/src/diffraction.c
index eb67e4a8..84206908 100644
--- a/src/diffraction.c
+++ b/src/diffraction.c
@@ -140,24 +140,12 @@ void get_diffraction(struct image *image, int na, int nb, int nc)
 	double bx, by, bz;
 	double cx, cy, cz;
 	double a, b, c, d;
-	struct molecule *mtmp;
 
 	/* Generate the array of reciprocal space vectors in image->qvecs */
 	if ( image->qvecs == NULL ) {
 		get_ewald(image);
 	}
 
-	/* FIXME: Nasty */
-	mtmp = load_molecule();
-	if ( image->molecule == NULL ) {
-		image->molecule = mtmp;
-	} else {
-		int i;
-		for ( i=0; i<32; i++ ) {
-			image->molecule->species[i] = mtmp->species[i];
-		}
-		image->molecule->n_species = mtmp->n_species;
-	}
 	if ( image->molecule == NULL ) return;
 
 	cell_get_cartesian(image->molecule->cell, &ax, &ay, &az,
diff --git a/src/dirax.c b/src/dirax.c
index 2c0675da..4161c79c 100644
--- a/src/dirax.c
+++ b/src/dirax.c
@@ -69,13 +69,7 @@ static void dirax_parseline(const char *line, struct image *image)
 		if ( line[i] == 'R' ) rf = 1;
 		if ( (line[i] == 'D') && rf ) {
 			image->dirax_read_cell = 1;
-			if ( image->molecule == NULL ) {
-				image->molecule = malloc(sizeof(struct molecule));
-			} else if ( image->molecule->cell ) {
-				free(image->molecule->cell);
-				free(image->molecule);
-			}
-			image->molecule->cell = cell_new();
+			image->indexed_cell = cell_new();
 			return;
 		}
 		i++;
@@ -86,7 +80,7 @@ static void dirax_parseline(const char *line, struct image *image)
 		/* First row of unit cell values */
 		float ax, ay, az, d;
 		sscanf(line, "%f %f %f %f %f %f", &d, &d, &d, &ax, &ay, &az);
-		cell_set_cartesian_a(image->molecule->cell,
+		cell_set_cartesian_a(image->indexed_cell,
 		                     ax*1e-10, ay*1e-10, az*1e-10);
 		image->dirax_read_cell++;
 		return;
@@ -94,7 +88,7 @@ static void dirax_parseline(const char *line, struct image *image)
 		/* First row of unit cell values */
 		float bx, by, bz, d;
 		sscanf(line, "%f %f %f %f %f %f", &d, &d, &d, &bx, &by, &bz);
-		cell_set_cartesian_b(image->molecule->cell,
+		cell_set_cartesian_b(image->indexed_cell,
 		                     bx*1e-10, by*1e-10, bz*1e-10);
 		image->dirax_read_cell++;
 		return;
@@ -102,7 +96,7 @@ static void dirax_parseline(const char *line, struct image *image)
 		/* First row of unit cell values */
 		float cx, cy, cz, d;
 		sscanf(line, "%f %f %f %f %f %f", &d, &d, &d, &cx, &cy, &cz);
-		cell_set_cartesian_c(image->molecule->cell,
+		cell_set_cartesian_c(image->indexed_cell,
 		                     cx*1e-10, cy*1e-10, cz*1e-10);
 		STATUS("Read a direct space unit cell from DirAx\n");
 		/* FIXME: Do something */
diff --git a/src/image.h b/src/image.h
index 0f64d248..59435b4e 100644
--- a/src/image.h
+++ b/src/image.h
@@ -78,6 +78,7 @@ struct image {
 	struct rvec		*qvecs;
 	double			*twotheta;
 	struct molecule		*molecule;
+	UnitCell		*indexed_cell;
 
 	struct quaternion	orientation;
 
diff --git a/src/index.c b/src/index.c
index b3c0b7e7..a20ad3a5 100644
--- a/src/index.c
+++ b/src/index.c
@@ -96,7 +96,16 @@ void index_pattern(struct image *image, int no_index, int use_dirax)
 	}
 
 	if ( use_dirax ) {
+
+		UnitCell *new_cell;
+
 		run_dirax(image, no_index);
+
+		new_cell = match_cell(image->indexed_cell,
+		                      image->molecule->cell);
+
+		if ( new_cell != NULL ) image->indexed_cell = new_cell;
+
 		return;
 	}
 }
diff --git a/src/indexamajig.c b/src/indexamajig.c
index 1a94882b..245d1a81 100644
--- a/src/indexamajig.c
+++ b/src/indexamajig.c
@@ -30,6 +30,7 @@
 #include "peaks.h"
 #include "diffraction.h"
 #include "detector.h"
+#include "sfac.h"
 
 
 static void show_help(const char *s)
@@ -132,8 +133,9 @@ int main(int argc, char *argv[])
 		chomp(line);
 
 		image.features = NULL;
-		image.molecule = NULL;
+		image.molecule = load_molecule();
 		image.data = NULL;
+		image.indexed_cell = NULL;
 
 		STATUS("Processing '%s'\n", line);
 
@@ -162,10 +164,8 @@ int main(int argc, char *argv[])
 			if ( config_noindex ) goto done;
 
 			/* Calculate orientation matrix (by magic) */
-			index_pattern(&image, config_noindex,
-			              config_dirax);
-
-			if ( image.molecule == NULL ) goto done;
+			index_pattern(&image, config_noindex, config_dirax);
+			if ( image.indexed_cell == NULL ) goto done;
 
 			if ( config_nearbragg || config_simulate ) {
 
@@ -186,7 +186,7 @@ int main(int argc, char *argv[])
 
 			if ( config_nearbragg ) {
 				/* Read h,k,l,I */
-				output_intensities(&image);
+				output_intensities(&image, image.indexed_cell);
 			}
 
 			if ( config_simulate ) {
diff --git a/src/intensities.c b/src/intensities.c
index 8b100d9f..d5b3604c 100644
--- a/src/intensities.c
+++ b/src/intensities.c
@@ -49,7 +49,7 @@ static int sum_nearby_points(int16_t *data, int width, int x, int y)
 }
 
 
-void output_intensities(struct image *image)
+void output_intensities(struct image *image, UnitCell *cell)
 {
 	int x, y;
 	double ax, ay, az;
@@ -59,9 +59,7 @@ void output_intensities(struct image *image)
 	int n_hits = 0;
 	int i;
 
-	cell_get_cartesian(image->molecule->cell, &ax, &ay, &az,
-		                                  &bx, &by, &bz,
-		                                  &cx, &cy, &cz);
+	cell_get_cartesian(cell, &ax, &ay, &az, &bx, &by, &bz, &cx, &cy, &cz);
 
 	for ( x=0; x<image->width; x++ ) {
 	for ( y=0; y<image->height; y++ ) {
diff --git a/src/intensities.h b/src/intensities.h
index 8f863287..10664ab2 100644
--- a/src/intensities.h
+++ b/src/intensities.h
@@ -17,7 +17,8 @@
 #define INTENSITIES_H
 
 #include "image.h"
+#include "cell.h"
 
-extern void output_intensities(struct image *image);
+extern void output_intensities(struct image *image, UnitCell *cell);
 
 #endif	/* INTENSITIES_H */
diff --git a/src/pattern_sim.c b/src/pattern_sim.c
index 9e4abf90..7a27ceb6 100644
--- a/src/pattern_sim.c
+++ b/src/pattern_sim.c
@@ -218,7 +218,7 @@ int main(int argc, char *argv[])
 	image.camera_len = 0.05;  /* 5 cm (front CCD can move from 5cm-20cm) */
 	image.resolution = 13333.3; /* 75 micron pixel size */
 	image.lambda = ph_en_to_lambda(eV_to_J(2.0e3));  /* Wavelength */
-	image.molecule = NULL;
+	image.molecule = load_molecule();
 
 	/* Splurge a few useful numbers */
 	STATUS("Wavelength is %f nm\n", image.lambda/1.0e-9);
@@ -263,7 +263,7 @@ int main(int argc, char *argv[])
 		             !config_nobloom);
 
 		if ( config_nearbragg ) {
-			output_intensities(&image);
+			output_intensities(&image, image.molecule->cell);
 		}
 
 		if ( !config_noimages ) {