Loads of lattice stuff

11 files changed, 205 insertions, 60 deletions

diff --git a/src/Makefile.am b/src/Makefile.am
index eca48660..f5ce4e6c 100644
--- a/src/Makefile.am
+++ b/src/Makefile.am
@@ -2,5 +2,6 @@ bin_PROGRAMS = pattern_sim
 
 AM_CFLAGS = -Wall -g @CFLAGS@
 
-pattern_sim_SOURCES = main.c diffraction.c utils.c image.c cell.c hdf5-file.c
+pattern_sim_SOURCES = main.c diffraction.c utils.c image.c cell.c hdf5-file.c \
+                      ewald.c
 pattern_sim_LDADD = @LIBS@
diff --git a/src/diffraction.c b/src/diffraction.c
index 93c6a22c..b2da4622 100644
--- a/src/diffraction.c
+++ b/src/diffraction.c
@@ -17,43 +17,64 @@
 #include "image.h"
 #include "utils.h"
 #include "cell.h"
-
-
-static 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;
-}
+#include "ewald.h"
 
 
 void get_diffraction(struct image *image, UnitCell *cell)
 {
+	int x, y;
+	double ax, ay, az;
+	double bx, by, bz;
+	double cx, cy, cz;	
+	int na = 64;
+	int nb = 64;
+	int nc = 64;
+
+	/* Generate the array of reciprocal space vectors in image->qvecs */
+	get_ewald(image);
+	
+	cell_get_cartesian(cell, &ax, &ay, &az,
+		                 &bx, &by, &bz,
+		                 &cx, &cy, &cz);
+	
+	image->sfacs = malloc(image->width * image->height * sizeof(double));
+	
+	for ( x=0; x<image->width; x++ ) {
+	for ( y=0; y<image->height; y++ ) {
+	
+		struct threevec q;
+		struct threevec Udotq;
+		double f1, f2, f3;
+		                
+		q = image->qvecs[x + image->width*y];
+		
+		Udotq.u = (ax*q.u + ay*q.v + az*q.w)/2.0;
+		Udotq.v = (bx*q.u + by*q.v + bz*q.w)/2.0;
+		Udotq.w = (cx*q.u + cy*q.v + cz*q.w)/2.0;
+		
+		if ( na > 1 ) {
+			f1 = sin(2.0*M_PI*(double)na*Udotq.u)
+			       / sin(2.0*M_PI*Udotq.u);
+		} else {
+			f1 = 1.0;
+		}
+		
+		if ( nb > 1 ) {
+			f2 = sin(2.0*M_PI*(double)nb*Udotq.v)
+			       / sin(2.0*M_PI*Udotq.v);
+		} else {
+			f2 = 1.0;
+		}
+		
+		if ( nc > 1 ) {
+			f3 = sin(2.0*M_PI*(double)nc*Udotq.w)
+			       / sin(2.0*M_PI*Udotq.w);
+		} else {
+			f3 = 1.0;
+		}
+		
+		image->sfacs[x + image->width*y] = f1 * f2 * f3;
 
+	}
+	}
 }
diff --git a/src/ewald.c b/src/ewald.c
new file mode 100644
index 00000000..9d8a0450
--- /dev/null
+++ b/src/ewald.c
@@ -0,0 +1,57 @@
+/*
+ * ewald.c
+ *
+ * Calculate q-vector arrays
+ *
+ * (c) 2007-2009 Thomas White <thomas.white@desy.de>
+ *
+ * pattern_sim - Simulate diffraction patterns from small crystals
+ *
+ */
+
+
+#include <stdlib.h>
+#include <math.h>
+#include <stdio.h>
+
+#include "image.h"
+#include "utils.h"
+#include "cell.h"
+
+
+void get_ewald(struct image *image)
+{
+	int x, y;
+	double k;  /* Wavenumber */
+	
+	k = 1/image->lambda;
+	
+	image->qvecs = malloc(image->width * image->height
+	                       * sizeof(struct threevec));
+	
+	for ( x=0; x<image->width; x++ ) {
+	for ( y=0; y<image->height; y++ ) {
+	
+		double rx, ry, r;
+		double twothetax, twothetay, twotheta;
+		double qx, qy, qz;
+		
+		rx = ((double)x - image->x_centre) / image->resolution;
+		ry = ((double)y - image->y_centre) / image->resolution;
+		r = sqrt(pow(rx, 2.0) + pow(ry, 2.0));
+		
+		twothetax = atan2(rx, image->camera_len);
+		twothetay = atan2(ry, image->camera_len);	
+		twotheta = atan2(r, image->camera_len);
+		
+		qx = k * sin(twothetax);
+		qy = k * sin(twothetay);
+		qz = k - k * cos(twotheta);
+		
+		image->qvecs[x + image->width*y].u = qx;
+		image->qvecs[x + image->width*y].v = qy;
+		image->qvecs[x + image->width*y].w = qz;
+	
+	}
+	}
+}
diff --git a/src/ewald.h b/src/ewald.h
new file mode 100644
index 00000000..ffa6e162
--- /dev/null
+++ b/src/ewald.h
@@ -0,0 +1,23 @@
+/*
+ * ewald.h
+ *
+ * Calculate q-vector arrays
+ *
+ * (c) 2007-2009 Thomas White <thomas.white@desy.de>
+ *
+ * pattern_sim - Simulate diffraction patterns from small crystals
+ *
+ */
+
+#ifdef HAVE_CONFIG_H
+#include <config.h>
+#endif
+
+#ifndef EWALD_H
+#define EWALD_H
+
+#include "image.h"
+
+extern void get_ewald(struct image *image);
+
+#endif	/* EWALD_H */
diff --git a/src/hdf5-file.c b/src/hdf5-file.c
index ca1d5450..6eb84aeb 100644
--- a/src/hdf5-file.c
+++ b/src/hdf5-file.c
@@ -22,7 +22,7 @@
 #include "image.h"
 
 
-int hdf5_write(const char *filename, const uint16_t *data,
+int hdf5_write(const char *filename, const double *data,
                int width, int height)
 {
 	hid_t fh, gh, sh, dh;	/* File, group, dataspace and data handles */
@@ -49,7 +49,7 @@ int hdf5_write(const char *filename, const uint16_t *data,
 	max_size[1] = height;
 	sh = H5Screate_simple(2, size, max_size);
 
-	dh = H5Dcreate(gh, "data", H5T_NATIVE_UINT16, sh,
+	dh = H5Dcreate(gh, "data", H5T_NATIVE_FLOAT, sh,
 	               H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT);
 	if ( dh < 0 ) {
 		fprintf(stderr, "Couldn't create dataset\n");
@@ -63,7 +63,7 @@ int hdf5_write(const char *filename, const uint16_t *data,
 	                           (int)size[1], (int)size[0],
 	                           (int)max_size[1], (int)max_size[0]);
 
-	r = H5Dwrite(dh, H5T_NATIVE_UINT16, H5S_ALL,
+	r = H5Dwrite(dh, H5T_NATIVE_DOUBLE, H5S_ALL,
 	             H5S_ALL, H5P_DEFAULT, data);
 	if ( r < 0 ) {
 		fprintf(stderr, "Couldn't write data\n");
@@ -72,6 +72,8 @@ int hdf5_write(const char *filename, const uint16_t *data,
 		return 1;
 	}
 
+	H5Gclose(gh);
+	H5Dclose(dh);
 	H5Fclose(fh);
 
 	return 0;
diff --git a/src/hdf5-file.h b/src/hdf5-file.h
index a263368a..7f74efa3 100644
--- a/src/hdf5-file.h
+++ b/src/hdf5-file.h
@@ -18,7 +18,7 @@
 
 #include <stdint.h>
 
-extern int hdf5_write(const char *filename, const uint16_t *data,
+extern int hdf5_write(const char *filename, const double *data,
                       int width, int height);
 
 extern int hdf5_read(struct image *image, const char *filename);
diff --git a/src/image.c b/src/image.c
index 874d1f4a..04d0f35d 100644
--- a/src/image.c
+++ b/src/image.c
@@ -14,6 +14,7 @@
 #include <stdlib.h>
 #include <assert.h>
 #include <math.h>
+#include <stdio.h>
 
 #include "image.h"
 #include "utils.h"
diff --git a/src/image.h b/src/image.h
index c596abfe..35bd5da9 100644
--- a/src/image.h
+++ b/src/image.h
@@ -46,10 +46,22 @@ struct imagefeature {
 /* An opaque type representing a list of image features */
 typedef struct _imagefeaturelist ImageFeatureList;
 
+
+/* A 3D vector in reciprocal space */
+struct threevec
+{
+	double   u;
+	double   v;
+	double   w;
+};
+
+
 /* Structure describing an image */
 struct image {
 
 	uint16_t		*data;
+	double			*sfacs;
+	struct threevec		*qvecs;
 
 	/* Radians.  Defines where the pattern lies in reciprocal space */
 	double			tilt;
diff --git a/src/main.c b/src/main.c
index 8ab11c85..ee02d340 100644
--- a/src/main.c
+++ b/src/main.c
@@ -19,7 +19,7 @@
 #include <unistd.h>
 
 #include "image.h"
-#include "relrod.h"
+#include "diffraction.h"
 #include "cell.h"
 #include "utils.h"
 #include "hdf5-file.h"
@@ -39,11 +39,9 @@ static void main_show_help(const char *s)
 
 int main(int argc, char *argv[])
 {
-	int c, i;
+	int c;
 	UnitCell *cell;
 	struct image image;
-	int nrefl;
-	float t;
 
 	while ((c = getopt(argc, argv, "h")) != -1) {
 
@@ -70,28 +68,21 @@ int main(int argc, char *argv[])
 	image.width = 512;
 	image.height = 512;
 	image.omega = deg2rad(40.0);
+	image.tilt = deg2rad(0.0);
 	image.fmode = FORMULATION_CLEN;
 	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);
-
-	for ( t=0.0; t<180.0; t+=10.0 ) {
-
-		char filename[32];
-
-		memset(image.data, 0, 512*512*2);
-		image.tilt = deg2rad(t);
-
-		get_diffraction(&image, cell);
-
-		/* Write the output file */
-		snprintf(filename, 32, "simulated-%.0f.h5", t);
-		hdf5_write(filename, image.data, image.width, image.height);
-
-	}
+	image.qvecs = NULL;
+	image.sfacs = NULL;
+	image.data = NULL;
+	
+	get_diffraction(&image, cell);
+
+	/* Write the output file */
+	hdf5_write("results/sim.h5", image.sfacs, image.width, image.height);
 
 	return 0;
 }
diff --git a/src/utils.c b/src/utils.c
index 38855671..fa239f90 100644
--- a/src/utils.c
+++ b/src/utils.c
@@ -122,3 +122,37 @@ int sign(double a)
 	if ( a > 0 ) return +1;
 	return 0;
 }
+
+
+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;
+}
diff --git a/src/utils.h b/src/utils.h
index a60abe7e..f1b72af8 100644
--- a/src/utils.h
+++ b/src/utils.h
@@ -33,7 +33,10 @@ extern double distance3d(double x1, double y1, double z1,
 extern size_t skipspace(const char *s);
 extern void chomp(char *s);
 extern int sign(double a);
-
+extern void mapping_rotate(double x, double y, double z,
+                           double *ddx, double *ddy, double *ddz,
+                           double omega, double tilt);
+                           
 #define rad2deg(a) ((a)*180/M_PI)
 #define deg2rad(a) ((a)*M_PI/180)