Calculate reflections using reciprocal cell

7 files changed, 140 insertions, 95 deletions

diff --git a/src/Makefile.am b/src/Makefile.am
index b480cedd..46a5e9d9 100644
--- a/src/Makefile.am
+++ b/src/Makefile.am
@@ -1,8 +1,8 @@
 bin_PROGRAMS = template_index simulate_patterns
 
 template_index_SOURCES = main.c relrod.c utils.c image.c cell.c
-template_index_LDADD = @LIBS@ -lm
+template_index_LDADD = @LIBS@ -lm -lgsl -lgslcblas
 AM_CFLAGS = -Wall -g
 
 simulate_patterns_SOURCES = sim-main.c relrod.c utils.c image.c cell.c
-simulate_patterns_LDADD = @LIBS@ -lm
+simulate_patterns_LDADD = @LIBS@ -lm -lgsl -lgslcblas
diff --git a/src/cell.c b/src/cell.c
index 6fa91cd3..81638016 100644
--- a/src/cell.c
+++ b/src/cell.c
@@ -16,6 +16,9 @@
 #include <math.h>
 #include <stdlib.h>
 #include <stdio.h>
+#include <gsl/gsl_matrix.h>
+#include <gsl/gsl_blas.h>
+#include <gsl/gsl_linalg.h>
 
 #include "cell.h"
 #include "utils.h"
@@ -74,9 +77,9 @@ static void cell_update_crystallographic(UnitCell *cell)
 	cell->gamma = angle_between(cell->ax, cell->ay, cell->az,
 					cell->bx, cell->by, cell->bz);
 
-	printf("a=%f nm\n", cell->a);
-	printf("b=%f nm\n", cell->b);
-	printf("c=%f nm\n", cell->c);
+	printf("a=%f nm\n", cell->a/1e9);
+	printf("b=%f nm\n", cell->b/1e9);
+	printf("c=%f nm\n", cell->c/1e9);
 	printf("alpha = %f deg\n", rad2deg(cell->alpha));
 	printf(" beta = %f deg\n", rad2deg(cell->beta));
 	printf("gamma = %f deg\n", rad2deg(cell->gamma));
@@ -104,7 +107,7 @@ UnitCell *cell_new()
 
 
 void cell_set_parameters(UnitCell *cell, double a, double b, double c,
-				double alpha, double beta, double gamma)
+                         double alpha, double beta, double gamma)
 {
 	if ( !cell ) return;
 
@@ -135,7 +138,7 @@ void cell_set_cartesian(UnitCell *cell,
 
 
 UnitCell *cell_new_from_parameters(double a, double b, double c,
-				double alpha, double beta, double gamma)
+                                   double alpha, double beta, double gamma)
 {
 	UnitCell *cell;
 
@@ -149,9 +152,9 @@ UnitCell *cell_new_from_parameters(double a, double b, double c,
 
 
 void cell_get_cartesian(UnitCell *cell,
-			double *ax, double *ay, double *az,
-			double *bx, double *by, double *bz,
-			double *cx, double *cy, double *cz)
+                        double *ax, double *ay, double *az,
+                        double *bx, double *by, double *bz,
+                        double *cx, double *cy, double *cz)
 {
 	if ( !cell ) return;
 
@@ -159,3 +162,51 @@ void cell_get_cartesian(UnitCell *cell,
 	*bx = cell->bx;  *by = cell->by;  *bz = cell->bz;
 	*cx = cell->cx;  *cy = cell->cy;  *cz = cell->cz;
 }
+
+
+void cell_get_reciprocal(UnitCell *cell,
+                         double *asx, double *asy, double *asz,
+                         double *bsx, double *bsy, double *bsz,
+                         double *csx, double *csy, double *csz)
+{
+	int s;
+	gsl_matrix *m;
+	gsl_matrix *inv;
+	gsl_permutation *perm;
+
+	m = gsl_matrix_alloc(3, 3);
+	gsl_matrix_set(m, 0, 0, cell->ax);
+	gsl_matrix_set(m, 0, 1, cell->bx);
+	gsl_matrix_set(m, 0, 2, cell->cx);
+	gsl_matrix_set(m, 1, 0, cell->ay);
+	gsl_matrix_set(m, 1, 1, cell->by);
+	gsl_matrix_set(m, 1, 2, cell->cy);
+	gsl_matrix_set(m, 2, 0, cell->az);
+	gsl_matrix_set(m, 2, 1, cell->bz);
+	gsl_matrix_set(m, 2, 2, cell->cz);
+
+	/* 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);
+
+	*asx = gsl_matrix_get(inv, 0, 0);
+	*bsx = gsl_matrix_get(inv, 0, 1);
+	*csx = gsl_matrix_get(inv, 0, 2);
+	*asy = gsl_matrix_get(inv, 1, 0);
+	*bsy = gsl_matrix_get(inv, 1, 1);
+	*csy = gsl_matrix_get(inv, 1, 2);
+	*asz = gsl_matrix_get(inv, 2, 0);
+	*bsz = gsl_matrix_get(inv, 2, 1);
+	*csz = gsl_matrix_get(inv, 2, 2);
+
+	printf("a* = %f %f %f\n", *asx/1e9, *asy/1e9, *asz/1e9);
+	printf("b* = %f %f %f\n", *bsx/1e9, *bsy/1e9, *bsz/1e9);
+	printf("c* = %f %f %f\n", *csx/1e9, *csy/1e9, *csz/1e9);
+}
diff --git a/src/cell.h b/src/cell.h
index cd31aa66..196b9521 100644
--- a/src/cell.h
+++ b/src/cell.h
@@ -19,9 +19,9 @@
 typedef struct {
 
 	/* Crystallographic representation */
-	double a;	/* nm */
-	double b;	/* nm */
-	double c;	/* nm */
+	double a;	/* m */
+	double b;	/* m */
+	double c;	/* m */
 	double alpha;	/* Radians */
 	double beta;	/* Radians */
 	double gamma;	/* Radians */
@@ -40,7 +40,7 @@ typedef struct {
 
 extern UnitCell *cell_new(void);
 
-/* Lengths in nm, angles in radians */
+/* Lengths in m, angles in radians */
 extern UnitCell *cell_new_from_parameters(double a, double b, double c,
 				double alpha, double beta, double gamma);
 
@@ -57,4 +57,9 @@ extern void cell_get_cartesian(UnitCell *cell,
                                double *bx, double *by, double *bz,
                                double *cx, double *cy, double *cz);
 
+extern void cell_get_reciprocal(UnitCell *cell,
+                               double *asx, double *asy, double *asz,
+                               double *bsx, double *bsy, double *bsz,
+                               double *csx, double *csy, double *csz);
+
 #endif	/* CELL_H */
diff --git a/src/image.c b/src/image.c
index 8780c0dd..d3e65685 100644
--- a/src/image.c
+++ b/src/image.c
@@ -69,8 +69,8 @@ ImageList *image_list_new()
 }
 
 
-void image_add_feature_reflection(ImageFeatureList *flist, double x, double y,
-                                  struct image *parent, double intensity)
+void image_add_feature(ImageFeatureList *flist, double x, double y,
+                       struct image *parent, double intensity)
 {
 	if ( flist->features ) {
 		flist->features = realloc(flist->features,
@@ -93,13 +93,6 @@ void image_add_feature_reflection(ImageFeatureList *flist, double x, double y,
 }
 
 
-void image_add_feature(ImageFeatureList *flist, double x, double y,
-                       struct image *parent, double intensity)
-{
-	image_add_feature_reflection(flist, x, y, parent, intensity);
-}
-
-
 ImageFeatureList *image_feature_list_new()
 {
 	ImageFeatureList *flist;
diff --git a/src/image.h b/src/image.h
index 69b8cc79..fc8a1f2f 100644
--- a/src/image.h
+++ b/src/image.h
@@ -82,10 +82,13 @@ struct image {
 typedef struct _imagelist ImageList;
 
 
+/* Image lists */
 extern ImageList *image_list_new(void);
 
 extern int image_add(ImageList *list, struct image *image);
 
+
+/* Feature lists */
 extern ImageFeatureList *image_feature_list_new(void);
 
 extern void image_feature_list_free(ImageFeatureList *flist);
@@ -93,11 +96,6 @@ extern void image_feature_list_free(ImageFeatureList *flist);
 extern void image_add_feature(ImageFeatureList *flist, double x, double y,
                               struct image *parent, double intensity);
 
-extern void image_add_feature_reflection(ImageFeatureList *flist,
-                                         double x, double y,
-                                         struct image *parent,
-                                         double intensity);
-
 extern struct imagefeature *image_feature_closest(ImageFeatureList *flist,
                                                   double x, double y, double *d,
                                                   int *idx);
diff --git a/src/relrod.c b/src/relrod.c
index df4c4df3..e344133a 100644
--- a/src/relrod.c
+++ b/src/relrod.c
@@ -62,14 +62,18 @@ void get_reflections(struct image *image, UnitCell *cell)
 	double kx, ky, kz; /* Electron wavevector ("normal" times 1/lambda) */
 	double ux, uy, uz; /* "up" vector */
 	double rx, ry, rz; /* "right" vector */
-	double ax, ay, az; /* "a" lattice parameter */
-	double bx, by, bz; /* "b" lattice parameter */
-	double cx, cy, cz; /* "c" lattice parameter */
+	double asx, asy, asz; /* "a*" lattice parameter */
+	double bsx, bsy, bsz; /* "b*" lattice parameter */
+	double csx, csy, csz; /* "c*" lattice parameter */
 	signed int h, k, l;
 
-	flist = image_feature_list_new();
-	cell_get_cartesian(cell, &ax, &ay, &az, &bx, &by, &bz, &cx, &cy, &cz);
+	/* Get the reciprocal unit cell */
+	cell_get_reciprocal(cell, &asx, &asy, &asz,
+	                          &bsx, &bsy, &bsz,
+	                          &csx, &csy, &csz);
 
+	/* Prepare list and some parameters */
+	flist = image_feature_list_new();
 	tilt = image->tilt;
 	omega = image->omega;
 	wavenumber = 1.0/image->lambda;
@@ -96,9 +100,9 @@ void get_reflections(struct image *image, UnitCell *cell)
 		double g_sq, gn;
 
 		/* Get the coordinates of the reciprocal lattice point */
-		xl = h*ax + k*bx + l*cx;
-		yl = h*ay + k*by + l*cy;
-		zl = h*az + k*bz + l*cz;
+		xl = h*asx + k*bsx + l*csx;
+		yl = h*asy + k*bsy + l*csy;
+		zl = h*asz + k*bsz + l*csz;
 		g_sq = modulus_squared(xl, yl, zl);
 		gn = xl*nx + yl*ny + zl*nz;
 
@@ -119,6 +123,7 @@ void get_reflections(struct image *image, UnitCell *cell)
 			double gx, gy, gz;
 			double theta;
 			double x, y;
+			double dx, dy, psi;
 
 			/* Determine the intersection point */
 			xi = xl + s*nx;  yi = yl + s*ny;  zi = zl + s*nz;
@@ -130,56 +135,49 @@ void get_reflections(struct image *image, UnitCell *cell)
 					*  the sphere to the intersection */
 			theta = angle_between(-kx, -ky, -kz, gx, gy, gz);
 
-			if ( theta > 0.0 ) {
-
-				double dx, dy, psi;
-
-				/* Calculate azimuth of point in image
-				 * (anticlockwise from +x) */
-				dx = xi*rx + yi*ry + zi*rz;
-				dy = xi*ux + yi*uy + zi*uz;
-				psi = atan2(dy, dx);
-
-				/* Get image coordinates from polar
-				 * representation */
-				if ( image->fmode == FORMULATION_CLEN ) {
-					x = image->camera_len*tan(theta)*cos(psi);
-					y = image->camera_len*tan(theta)*sin(psi);
-					x *= image->resolution;
-					y *= image->resolution;
-				} else if ( image->fmode==FORMULATION_PIXELSIZE ) {
-					x = tan(theta)*cos(psi) / image->lambda;
-					y = tan(theta)*sin(psi) / image->lambda;
-					x /= image->pixel_size;
-					y /= image->pixel_size;
-				} else {
-					fprintf(stderr,
-						"Unrecognised formulation mode "
-						"in get_reflections\n");
-					return;
-				}
-
-				x += image->x_centre;
-				y += image->y_centre;
-
-				/* Sanity check */
-				if ( (x>=0) && (x<image->width)
-				  && (y>=0) && (y<image->height) ) {
-
-					/* Record the reflection
-					 *  Intensity should be multiplied by
-					 *  relrod spike function except
-					 *  reprojected reflections aren't used
-					 *  quantitatively for anything
-					 */
-					image_add_feature_reflection(flist, x,
-						y, image, 1.0);
-
-				} /* else it's outside the picture somewhere */
-
-			} /* else this is the central beam */
-
-		}
+			/* Calculate azimuth of point in image
+			 * (anticlockwise from +x) */
+			dx = xi*rx + yi*ry + zi*rz;
+			dy = xi*ux + yi*uy + zi*uz;
+			psi = atan2(dy, dx);
+
+			/* Get image coordinates from polar
+			 * representation */
+			if ( image->fmode == FORMULATION_CLEN ) {
+				x = image->camera_len*tan(theta)*cos(psi);
+				y = image->camera_len*tan(theta)*sin(psi);
+				x *= image->resolution;
+				y *= image->resolution;
+			} else if ( image->fmode==FORMULATION_PIXELSIZE ) {
+				x = tan(theta)*cos(psi) / image->lambda;
+				y = tan(theta)*sin(psi) / image->lambda;
+				x /= image->pixel_size;
+				y /= image->pixel_size;
+			} else {
+				fprintf(stderr,
+					"Unrecognised formulation mode "
+					"in get_reflections\n");
+				return;
+			}
+
+			x += image->x_centre;
+			y += image->y_centre;
+
+			/* Sanity check */
+			if ( (x>=0) && (x<image->width)
+			  && (y>=0) && (y<image->height) ) {
+
+				/* Record the reflection.
+				 * Intensity should be multiplied by relrod
+				 * spike function, except reprojected
+				 * reflections aren't used quantitatively for
+				 * anything. */
+				 printf("%i %i %i\n", h,k,l);
+				image_add_feature(flist, x, y, image, 1.0);
+
+			} /* else it's outside the picture somewhere */
+
+		} /* else reflection is not excited in this orientation */
 
 	}
 	}
diff --git a/src/sim-main.c b/src/sim-main.c
index 38ad45bc..23ba3e99 100644
--- a/src/sim-main.c
+++ b/src/sim-main.c
@@ -61,20 +61,20 @@ int main(int argc, char *argv[])
 	image.tilt = 0.0;
 	image.omega = 0.0;
 	image.fmode = FORMULATION_CLEN;
-	image.x_centre = 128;
-	image.y_centre = 128;
-	image.camera_len = 1.0; /* one metre */
-	image.resolution = 5120;
-	image.lambda = 0.6e-9;
+	image.x_centre = 255.5;
+	image.y_centre = 255.5;
+	image.camera_len = 0.2;  /* 20 cm */
+	image.resolution = 5120; /* 512 pixels in 10 cm */
+	image.lambda = 0.2e-9;   /* LCLS wavelength */
 	image.data = malloc(512*512*2);
 	image_add(list, &image);
 
-	cell = cell_new_from_parameters(1.0,
-	                                1.0,
-	                                1.0,
+	cell = cell_new_from_parameters(28.10e-9,
+	                                28.10e-9,
+	                                16.52e-9,
 	                                deg2rad(90.0),
 	                                deg2rad(90.0),
-	                                deg2rad(90.0));
+	                                deg2rad(120.0));
 
 	get_reflections(&image, cell);