Only calculate molecular transform at Bragg positions

6 files changed, 203 insertions, 49 deletions

diff --git a/src/diffraction.c b/src/diffraction.c
index 4ed13709..e112e7ce 100644
--- a/src/diffraction.c
+++ b/src/diffraction.c
@@ -76,55 +76,78 @@ static double lattice_factor(struct threevec q, double ax, double ay, double az,
 }
 
 
-/* Return structure factor for molecule 'mol' at energy 'en' (J/photon) at
- * scattering vector 'q' */
+/* Look up the structure factor for the nearest Bragg condition */
 static double complex molecule_factor(struct molecule *mol, struct threevec q,
-                                      double en)
+                                      double ax, double ay, double az,
+                                      double bx, double by, double bz,
+                                      double cx, double cy, double cz)
 {
-	int i;
-	double complex F = 0.0;
-	double s;
+	double hd, kd, ld;
+	signed int h, k, l;
 
-	/* s = sin(theta)/lambda = 1/2d = (1/d)/2.0 */
-	s = modulus(q.u, q.v, q.w) / 2.0;
+	hd = q.u * ax + q.v * ay + q.w * az;
+	kd = q.u * bx + q.v * by + q.w * bz;
+	ld = q.u * cx + q.v * cy + q.w * cz;
+	h = (signed int)nearbyint(hd);
+	k = (signed int)nearbyint(kd);
+	l = (signed int)nearbyint(ld);
 
-	/* Atoms are grouped by species for faster calculation */
-	for ( i=0; i<mol->n_species; i++ ) {
-
-		double complex sfac;
-		double complex contrib = 0.0;
-		struct mol_species *spec;
-		int j;
-
-		spec = mol->species[i];
+	return get_integral(mol->reflections, h, k, l);
+}
 
-		for ( j=0; j<spec->n_atoms; j++ ) {
 
-			double ph;
+static double complex water_factor(struct threevec q, double en)
+{
+	return 0.0;
+}
 
-			ph = q.u*spec->x[j] + q.v*spec->y[j] + q.w*spec->z[j];
 
-			/* Conversion from revolutions to radians is required */
-			contrib += cos(2.0*M_PI*ph) + I*sin(2.0*M_PI*ph);
+static void get_reflections_cached(struct molecule *mol, double en)
+{
+	char s[1024];
+	FILE *fh;
+	size_t r;
+
+	snprintf(s, 1023, "reflections-%ieV.cache", (int)J_to_eV(en));
+	fh = fopen(s, "rb");
+	if ( fh == NULL ) {
+		printf("No cache file found (looked for %s)\n", s);
+		goto calc;
+	}
 
-		}
+	mol->reflections = reflist_new();
+	r = fread(mol->reflections, sizeof(double complex), IDIM*IDIM*IDIM, fh);
+	if ( r <  IDIM*IDIM*IDIM ) {
+		printf("Found cache file (%s), but failed to read.\n", s);
+		goto calc;
+	}
 
-		sfac = get_sfac(spec->species, s, en);
-		F += sfac * contrib * exp(-2.0 * spec->B[j] * s);
+	printf("Read structure factors (at Bragg positions) from %s\n", s);
+	return;
 
+calc:
+	printf("Calculating structure factors at Bragg positions...\n");
+	mol->reflections = get_reflections(mol, en);
+	fh = fopen(s, "wb");
+	if ( fh == NULL ) {
+		printf("Failed to write cache file (%s)\n", s);
+		return;
 	}
 
-	return F;
-}
-
+	r = fwrite(mol->reflections, sizeof(double complex),
+	           IDIM*IDIM*IDIM, fh);
+	if ( r <  IDIM*IDIM*IDIM ) {
+		printf("Failed to write cache file (%s)\n", s);
+		return;
+	}
+	fclose(fh);
 
-static double complex water_factor(struct threevec q, double en)
-{
-	return 0.0;
+	printf("Successfully saved structure factors at Bragg positions to"
+	       " file %s\n", s);
 }
 
 
-void get_diffraction(struct image *image, UnitCell *cell)
+void get_diffraction(struct image *image)
 {
 	int x, y;
 	double ax, ay, az;
@@ -139,13 +162,17 @@ void get_diffraction(struct image *image, UnitCell *cell)
 		if ( image->molecule == NULL ) return;
 	}
 
-	cell_get_cartesian(cell, &ax, &ay, &az,
-		                 &bx, &by, &bz,
-		                 &cx, &cy, &cz);
+	cell_get_cartesian(image->molecule->cell, &ax, &ay, &az,
+		                                  &bx, &by, &bz,
+		                                  &cx, &cy, &cz);
 
 	image->sfacs = malloc(image->width * image->height
 	                      * sizeof(double complex));
 
+	if ( image->molecule->reflections == NULL ) {
+		get_reflections_cached(image->molecule, image->xray_energy);
+	}
+
 	progress_bar(0, image->width-1);
 	for ( x=0; x<image->width; x++ ) {
 	for ( y=0; y<image->height; y++ ) {
@@ -154,12 +181,13 @@ void get_diffraction(struct image *image, UnitCell *cell)
 		double complex f_molecule;
 		double complex f_water;
 		struct threevec q;
+		double complex val;
 
 		q = image->qvecs[x + image->width*y];
 
 		f_lattice = lattice_factor(q, ax,ay,az,bx,by,bz,cx,cy,cz);
 		f_molecule = molecule_factor(image->molecule, q,
-		                             image->xray_energy);
+		                             ax,ay,az,bx,by,bz,cx,cy,cz);
 
 		/* Nasty approximation follows */
 		if ( y == image->height/2 ) {
@@ -168,8 +196,8 @@ void get_diffraction(struct image *image, UnitCell *cell)
 			f_water = 0.0;
 		}
 
-		image->sfacs[x + image->width*y] = (f_lattice * f_molecule)
-		                                                  + f_water;
+		val = (f_molecule) + f_water;
+		image->sfacs[x + image->width*y] = val;
 
 	}
 	progress_bar(x, image->width-1);
diff --git a/src/diffraction.h b/src/diffraction.h
index c32e62c2..112818bc 100644
--- a/src/diffraction.h
+++ b/src/diffraction.h
@@ -19,6 +19,6 @@
 #include "image.h"
 #include "cell.h"
 
-extern void get_diffraction(struct image *image, UnitCell *cell);
+extern void get_diffraction(struct image *image);
 
 #endif	/* DIFFRACTION_H */
diff --git a/src/main.c b/src/main.c
index b7755ce6..f3aeadfc 100644
--- a/src/main.c
+++ b/src/main.c
@@ -41,7 +41,6 @@ static void main_show_help(const char *s)
 int main(int argc, char *argv[])
 {
 	int c, done;
-	UnitCell *cell;
 	struct image image;
 	char filename[1024];
 	int number = 1;
@@ -59,14 +58,6 @@ int main(int argc, char *argv[])
 
 	}
 
-	/* Define unit cell */
-	cell = cell_new_from_parameters(28.10e-9,
-	                                28.10e-9,
-	                                16.52e-9,
-	                                deg2rad(90.0),
-	                                deg2rad(90.0),
-	                                deg2rad(120.0));
-
 	/* Define image parameters */
 	image.width = 1024;
 	image.height = 1024;
@@ -124,7 +115,7 @@ again:
 	image.twotheta = NULL;
 	image.hdr = NULL;
 
-	get_diffraction(&image, cell);
+	get_diffraction(&image);
 	record_image(&image);
 
 	snprintf(filename, 1023, "results/sim-%i.h5", number);
diff --git a/src/sfac.c b/src/sfac.c
index 00d6c69d..c13ed3b8 100644
--- a/src/sfac.c
+++ b/src/sfac.c
@@ -272,6 +272,15 @@ struct molecule *load_molecule()
 	mol = malloc(sizeof(struct molecule));
 	if ( mol == NULL ) return NULL;
 	mol->n_species = 0;
+	mol->reflections = NULL;
+
+	/* FIXME: Read cell from file */
+	mol->cell = cell_new_from_parameters(28.10e-9,
+	                                     28.10e-9,
+	                                     16.52e-9,
+	                                     deg2rad(90.0),
+	                                     deg2rad(90.0),
+	                                     deg2rad(120.0));
 
 	fh = fopen("molecule.pdb", "r");
 	if ( fh == NULL ) {
@@ -370,3 +379,68 @@ struct molecule *load_molecule()
 
 	return mol;
 }
+
+
+double complex *get_reflections(struct molecule *mol, double en)
+{
+	double complex *reflections;
+	double asx, asy, asz;
+	double bsx, bsy, bsz;
+	double csx, csy, csz;
+	signed int h, k, l;
+
+	cell_get_reciprocal(mol->cell, &asx, &asy, &asz,
+	                               &bsx, &bsy, &bsz,
+	                               &csx, &csy, &csz);
+
+	reflections = reflist_new();
+
+	for ( h=-INDMAX; h<=INDMAX; h++ ) {
+	for ( k=-INDMAX; k<=INDMAX; k++ ) {
+	for ( l=-INDMAX; l<=INDMAX; l++ ) {
+
+		double complex F;
+		int i;
+		double s;
+
+		/* Atoms are grouped by species for faster calculation */
+		for ( i=0; i<mol->n_species; i++ ) {
+
+			double complex sfac;
+			double complex contrib = 0.0;
+			struct mol_species *spec;
+			int j;
+
+			spec = mol->species[i];
+
+			for ( j=0; j<spec->n_atoms; j++ ) {
+
+				double ph, u, v, w;
+
+				u = h*asx + k*bsx + l*csx;
+				v = h*asy + k*bsy + l*csy;
+				w = h*asz + k*bsz + l*csz;
+
+				ph = u*spec->x[j] + v*spec->y[j] + w*spec->z[j];
+
+				/* Conversion from revolutions to radians */
+				contrib += cos(-2.0*M_PI*ph)
+				                          + I*sin(-2.0*M_PI*ph);
+
+			}
+
+			sfac = get_sfac(spec->species, s, en);
+			F += sfac * contrib * exp(-2.0 * spec->B[j] * s);
+
+		}
+
+		integrate_reflection(reflections, h, k, l, F);
+
+	}
+	}
+	progress_bar((h+INDMAX+1), 2*INDMAX);
+	}
+	printf("\n");
+
+	return reflections;
+}
diff --git a/src/sfac.h b/src/sfac.h
index f88729f4..fc4b2bd9 100644
--- a/src/sfac.h
+++ b/src/sfac.h
@@ -18,6 +18,8 @@
 
 #include <complex.h>
 
+#include "cell.h"
+
 
 struct mol_species
 {
@@ -37,6 +39,13 @@ struct molecule
 	int n_species;
 	struct mol_species *species[32];
 
+	/* Unit cell */
+	UnitCell *cell;
+
+	/* Reflection intensities at Bragg positions */
+	double complex *reflections;
+
+	/* Offset to molecule's centre of scattering power */
 	double xc;
 	double yc;
 	double zc;
@@ -45,5 +54,6 @@ struct molecule
 
 extern double complex get_sfac(const char *n, double s, double en);
 extern struct molecule *load_molecule(void);
+extern double complex *get_reflections(struct molecule *mol, double en);
 
 #endif	/* SFAC_H */
diff --git a/src/utils.h b/src/utils.h
index 7940d77f..ca653973 100644
--- a/src/utils.h
+++ b/src/utils.h
@@ -17,6 +17,9 @@
 #endif
 
 #include <math.h>
+#include <complex.h>
+#include <string.h>
+#include <stdlib.h>
 
 
 /* Electron charge in C */
@@ -31,6 +34,10 @@
 /* Thomson scattering length (m) */
 #define THOMSON_LENGTH (2.81794e-15)
 
+/* Maxmimum index to go up to */
+#define INDMAX 20
+#define IDIM (INDMAX*2 +1)
+
 
 extern unsigned int biggest(signed int a, signed int b);
 extern unsigned int smallest(signed int a, signed int b);
@@ -69,4 +76,48 @@ extern void progress_bar(int val, int total);
 /* Joules to eV */
 #define J_to_eV(a) ((a)/ELECTRON_CHARGE)
 
+
+static inline void integrate_reflection(double complex *ref, signed int h,
+                                        signed int k, signed int l,
+                                        double complex i)
+{
+	int idx;
+
+	/* Not interested in central beam */
+	if ( (h==0) && (k==0) && (l==0) ) return;
+
+	if ( h < 0 ) h += IDIM;
+	if ( k < 0 ) k += IDIM;
+	if ( l < 0 ) l += IDIM;
+
+	idx = h + (IDIM*k) + (IDIM*IDIM*l);
+	ref[idx] += i;
+}
+
+
+static inline double complex get_integral(double complex *ref, signed int h,
+                                          signed int k, signed int l)
+{
+	int idx;
+
+	if ( h < 0 ) h += IDIM;
+	if ( k < 0 ) k += IDIM;
+	if ( l < 0 ) l += IDIM;
+
+	idx = h + (IDIM*k) + (IDIM*IDIM*l);
+	return ref[idx];
+}
+
+
+static inline double complex *reflist_new(void)
+{
+	double complex *r;
+	size_t r_size;
+	r_size = IDIM*IDIM*IDIM*sizeof(double complex);
+	r = malloc(r_size);
+	memset(r, 0, r_size);
+	return r;
+}
+
+
 #endif	/* UTILS_H */