Introduce integr_sim

9 files changed, 436 insertions, 59 deletions

diff --git a/src/Makefile.am b/src/Makefile.am
index 3a251dfb..a10e59c0 100644
--- a/src/Makefile.am
+++ b/src/Makefile.am
@@ -1,4 +1,4 @@
-bin_PROGRAMS = pattern_sim #integr_sim
+bin_PROGRAMS = pattern_sim integr_sim
 
 AM_CFLAGS = -Wall -g @CFLAGS@
 
@@ -6,5 +6,6 @@ pattern_sim_SOURCES = pattern_sim.c diffraction.c utils.c image.c cell.c \
                       hdf5-file.c ewald.c detector.c sfac.c
 pattern_sim_LDADD = @LIBS@
 
-integr_sim_SOURCES = integr_sim.c diffraction.c utils.c ewald.c sfac.c
+integr_sim_SOURCES = integr_sim.c diffraction.c utils.c ewald.c sfac.c \
+                     statistics.c cell.c
 integr_sim_LDADD = @LIBS@
diff --git a/src/diffraction.c b/src/diffraction.c
index e2d0061a..6351042d 100644
--- a/src/diffraction.c
+++ b/src/diffraction.c
@@ -93,7 +93,7 @@ static double complex molecule_factor(struct molecule *mol, struct threevec q,
 	k = (signed int)rint(kd);
 	l = (signed int)rint(ld);
 
-	r = get_integral(mol->reflections, h, k, l);
+	r = lookup_sfac(mol->reflections, h, k, l);
 
 	return r;
 }
diff --git a/src/image.h b/src/image.h
index 5206048a..e9db7cad 100644
--- a/src/image.h
+++ b/src/image.h
@@ -20,6 +20,8 @@
 #include <stdint.h>
 #include <complex.h>
 
+#include "utils.h"
+
 
 /* How is the scaling of the image described? */
 typedef enum {
@@ -57,15 +59,6 @@ struct threevec
 };
 
 
-struct quaternion
-{
-	double w;
-	double x;
-	double y;
-	double z;
-};
-
-
 /* Structure describing an image */
 struct image {
 
diff --git a/src/integr_sim.c b/src/integr_sim.c
new file mode 100644
index 00000000..08788277
--- /dev/null
+++ b/src/integr_sim.c
@@ -0,0 +1,232 @@
+/*
+ * integr_sim.c
+ *
+ * Test integration of intensities
+ *
+ * (c) 2006-2009 Thomas White <thomas.white@desy.de>
+ *
+ * Part of CrystFEL - crystallography with a FEL
+ *
+ */
+
+
+#ifdef HAVE_CONFIG_H
+#include <config.h>
+#endif
+
+#include <stdarg.h>
+#include <stdlib.h>
+#include <stdio.h>
+#include <string.h>
+#include <unistd.h>
+
+#include "cell.h"
+#include "image.h"
+#include "utils.h"
+#include "statistics.h"
+#include "sfac.h"
+
+
+static void main_show_help(const char *s)
+{
+	printf("Syntax: %s\n\n", s);
+	printf("Test relrod integration\n\n");
+	printf("  -h              Display this help message\n");
+}
+
+
+static void write_RvsQ(const char *name, double *ref, double *trueref,
+                       unsigned int *counts, double scale, UnitCell *cell)
+{
+	FILE *fh;
+	double smax, sbracket;
+	signed int h, k, l;
+
+	fh = fopen(name, "w");
+
+	smax = 0.0;
+	for ( h=-INDMAX; h<INDMAX; h++ ) {
+	for ( k=-INDMAX; k<INDMAX; k++ ) {
+	for ( l=-INDMAX; l<INDMAX; l++ ) {
+		double s = resolution(cell, h, k, l);
+		if ( (lookup_count(counts, h, k, l) > 0) && (s > smax) ) {
+			smax = s;
+		}
+	}
+	}
+	}
+
+	for ( sbracket=0.0; sbracket<smax; sbracket+=smax/10.0 ) {
+
+		double top = 0.0;
+		double bot = 0.0;
+		int n = 0;
+
+		for ( h=-INDMAX; h<INDMAX; h++ ) {
+		for ( k=-INDMAX; k<INDMAX; k++ ) {
+		for ( l=-INDMAX; l<INDMAX; l++ ) {
+
+			double s;
+			int c;
+			c = lookup_count(counts, h, k, l);
+			s = resolution(cell, h, k, l);
+			if ((s>=sbracket) && (s<sbracket+smax/10.0) && (c>0)) {
+
+				double obs, calc, obsi;
+
+				obs = lookup_intensity(ref, h, k, l);
+				calc = lookup_intensity(trueref, h, k, l);
+
+				obsi = obs / (double)c;
+				top += fabs(obsi/scale - calc);
+				bot += obsi/scale;
+				n++;
+			}
+
+		}
+		}
+		}
+
+		fprintf(fh, "%8.5f %8.5f %i\n", sbracket+smax/20.0, top/bot, n);
+
+	}
+	fclose(fh);
+}
+
+
+int main(int argc, char *argv[])
+{
+	int c;
+	int i;
+	struct image image;
+	double *ref, *trueref;
+	unsigned int *counts;
+	size_t ref_size;
+	signed int h, k, l;
+	double R, scale;
+	FILE *fh;
+
+	while ((c = getopt(argc, argv, "h")) != -1) {
+
+		switch ( c ) {
+
+			case 'h' : {
+				main_show_help(argv[0]);
+				return 0;
+			}
+
+		}
+
+	}
+
+	/* Define image parameters */
+	image.width = 1024;
+	image.height = 1024;
+	image.fmode = FORMULATION_CLEN;
+	image.x_centre = 512.5;
+	image.y_centre = 512.5;
+	image.camera_len = 0.05;  /* 5 cm (front CCD can move from 5cm-20cm) */
+	image.resolution = 13333.3; /* 75 micron pixel size */
+	image.xray_energy = eV_to_J(2.0e3); /* 2 keV energy */
+	image.lambda = ph_en_to_lambda(image.xray_energy);  /* Wavelength */
+	image.molecule = load_molecule();
+
+	/* Prepare array for integrated intensities */
+	ref = new_list_intensity();
+
+	/* Array for sample counts */
+	counts = new_list_count();
+
+	/* Calculate true intensities */
+	get_reflections_cached(image.molecule, image.xray_energy);
+	/* Complex structure factors now in image.molecule->reflections */
+
+	for ( i=1; i<=10e3; i++ ) {
+
+		#if 0
+		image.orientation = random_quaternion();
+
+		/* Calculate reflections using large smax
+		 * (rather than the actual value) */
+		//get_reflections(&image, cell, 1.0e9);
+
+		nrefl = image_feature_count(image.rflist);
+		for ( j=0; j<nrefl; j++ ) {
+
+			struct imagefeature *f;
+			double t, s, intensity, F;
+
+			f = image_get_feature(image.rflist, j);
+
+			t = 100.0e-9;  /* Thickness 100 nm */
+			s = f->s;      /* Get excitation error */
+			F = structure_factor(f->h, f->k, f->l);
+
+			/* Calculate intensity from this reflection */
+			intensity = pow( F * SINC(M_PI*t*s), 2.0);
+
+			if ( intensity < 0.1 ) continue;
+
+			if ( (f->h == 2) && (f->k == 2) && (f->l == 2) ) {
+				fprintf(fh1, "%f %f\n", s, intensity);
+			}
+
+			if ( (f->h == 15) && (f->k == 15) && (f->l == 15) ) {
+				fprintf(fh2, "%f %f\n", s, intensity);
+			}
+
+			integrate_reflection(ref, f->h, f->k, f->l, intensity);
+			add_count(counts, f->h, f->k, f->l, 1);
+		}
+		#endif
+
+		if ( i % 1000 == 0 ) {
+
+			int j;
+			double mean_counts;
+			int ctot = 0;
+			int nmeas = 0;
+			double ff;
+			char name[64];
+
+			for ( j=0; j<ref_size; j++ ) {
+				ctot += counts[j];
+				if ( counts[j] > 0 ) nmeas++;
+			}
+			mean_counts = (double)ctot/nmeas;
+
+			ff = lookup_intensity(ref, 2, 2, 2)
+			                        / lookup_count(counts, 2, 2, 2);
+
+			R = stat_r2(ref, trueref, counts, ref_size, &scale);
+			printf("%8i: R=%5.2f%% (sf=%7.4f)"
+			       " mean meas/refl=%5.2f,"
+			       " %i reflections measured. %f\n",
+			       i, R*100.0, scale, mean_counts, nmeas,
+			       ff);
+
+			/* Record graph of R against q for this N */
+			snprintf(name, 63, "R_vs_q-%i.dat", i);
+			write_RvsQ(name, ref, trueref, counts,
+			           scale, image.molecule->cell);
+
+		}
+
+	}
+
+	fh = fopen("reflections.dat", "w");
+	for ( h=-INDMAX; h<INDMAX; h++ ) {
+	for ( k=-INDMAX; k<INDMAX; k++ ) {
+	for ( l=-INDMAX; l<INDMAX; l++ ) {
+		int N;
+		N = lookup_count(counts, h, k, l);
+		if ( N == 0 ) continue;
+		double F = lookup_intensity(ref, h, k, l) / N;
+		fprintf(fh, "%3i %3i %3i %f\n", h, k, l, F);
+	}
+	}
+	}
+	fclose(fh);
+
+	return 0;
+}
diff --git a/src/list_tmp.h b/src/list_tmp.h
new file mode 100644
index 00000000..ffb5f393
--- /dev/null
+++ b/src/list_tmp.h
@@ -0,0 +1,74 @@
+static inline void LABEL(integrate)(TYPE *ref, signed int h,
+                                    signed int k, signed int l,
+                                    TYPE i)
+{
+	int idx;
+
+	if ( (abs(h) > INDMAX) || (abs(k) > INDMAX) || (abs(l) > INDMAX) ) {
+		printf("\nReflection %i %i %i is out of range!\n", h, k, l);
+		printf("You need to re-configure INDMAX, delete the reflection"
+		       " cache file and re-run.\n");
+		exit(1);
+	}
+
+	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 void LABEL(set)(TYPE *ref, signed int h,
+                              signed int k, signed int l,
+                              TYPE i)
+{
+	int idx;
+
+	if ( (abs(h) > INDMAX) || (abs(k) > INDMAX) || (abs(l) > INDMAX) ) {
+		printf("\nReflection %i %i %i is out of range!\n", h, k, l);
+		printf("You need to re-configure INDMAX, delete the reflection"
+		       " cache file and re-run.\n");
+		exit(1);
+	}
+
+	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 TYPE LABEL(lookup)(TYPE *ref, signed int h,
+                                 signed int k, signed int l)
+{
+	int idx;
+
+	if ( (abs(h) > INDMAX) || (abs(k) > INDMAX) || (abs(l) > INDMAX) ) {
+		printf("\nReflection %i %i %i is out of range!\n", h, k, l);
+		printf("You need to re-configure INDMAX, delete the reflection"
+		       " cache file and re-run.\n");
+		exit(1);
+	}
+
+	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 TYPE *LABEL(new_list)(void)
+{
+	TYPE *r;
+	size_t r_size;
+	r_size = IDIM*IDIM*IDIM*sizeof(TYPE);
+	r = malloc(r_size);
+	memset(r, 0, r_size);
+	return r;
+}
diff --git a/src/sfac.c b/src/sfac.c
index e70bc750..e25ee783 100644
--- a/src/sfac.c
+++ b/src/sfac.c
@@ -400,7 +400,7 @@ double complex *get_reflections(struct molecule *mol, double en)
 	                               &bsx, &bsy, &bsz,
 	                               &csx, &csy, &csz);
 
-	reflections = reflist_new();
+	reflections = new_list_sfac();
 
 	for ( h=-INDMAX; h<=INDMAX; h++ ) {
 	for ( k=-INDMAX; k<=INDMAX; k++ ) {
@@ -446,7 +446,7 @@ double complex *get_reflections(struct molecule *mol, double en)
 
 		}
 
-		integrate_reflection(reflections, h, k, l, F);
+		set_sfac(reflections, h, k, l, F);
 		//if ( (h!=0) || (k!=0) || (l!=0) ) {
 		//	tscat += cabs(F);
 		//} else {
@@ -477,7 +477,7 @@ void get_reflections_cached(struct molecule *mol, double en)
 		goto calc;
 	}
 
-	mol->reflections = reflist_new();
+	mol->reflections = new_list_sfac();
 	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);
diff --git a/src/statistics.c b/src/statistics.c
new file mode 100644
index 00000000..adadf907
--- /dev/null
+++ b/src/statistics.c
@@ -0,0 +1,67 @@
+/*
+ * statistics.c
+ *
+ * Structure-factor statistics
+ *
+ * (c) 2007-2009 Thomas White <thomas.white@desy.de>
+ *
+ * integr_sim - Test relrod integration
+ *
+ */
+
+
+#ifdef HAVE_CONFIG_H
+#include <config.h>
+#endif
+
+#include <math.h>
+#include <stdlib.h>
+
+
+/* By what (best-fitted) factor must the list "list" be multiplied by,
+ * if it were to be merged with "target"? */
+double stat_scale_intensity(double *obs, double *calc, unsigned int *c,
+                            int size)
+{
+	double top = 0.0;
+	double bot = 0.0;
+	int i;
+
+	for ( i=0; i<size; i++ ) {
+
+		if ( c[i] > 0 ) {
+			double obsi;
+			obsi = obs[i] / (double)c[i];
+			top += obsi * calc[i];
+			bot += calc[i] * calc[i];
+		} /* else reflection not measured so don't worry about it */
+
+	}
+
+	return top/bot;
+}
+
+
+double stat_r2(double *obs, double *calc, unsigned int *c, int size,
+               double *scalep)
+{
+	double top = 0.0;
+	double bot = 0.0;
+	double scale;
+	int i;
+	scale = stat_scale_intensity(obs, calc, c, size);
+	*scalep = scale;
+
+	for ( i=0; i<size; i++ ) {
+
+		if ( c[i] > 0 ) {
+			double obsi;
+			obsi = obs[i] / (double)c[i];
+			top += fabs(obsi/scale - calc[i]);
+			bot += obsi/scale;
+		}
+
+	} /* else reflection not measured so don't worry about it */
+
+	return top/bot;
+}
diff --git a/src/statistics.h b/src/statistics.h
new file mode 100644
index 00000000..bc8184ca
--- /dev/null
+++ b/src/statistics.h
@@ -0,0 +1,22 @@
+/*
+ * statistics.h
+ *
+ * Structure-factor statistics
+ *
+ * (c) 2007-2009 Thomas White <thomas.white@desy.de>
+ *
+ * integr_sim - Test relrod integration
+ *
+ */
+
+#ifdef HAVE_CONFIG_H
+#include <config.h>
+#endif
+
+#ifndef STATISTICS_H
+#define STATISTICS_H
+
+double stat_r2(double *obs, double *calc, unsigned int *c, int size,
+               double *scalep);
+
+#endif	/* STATISTICS_H */
diff --git a/src/utils.h b/src/utils.h
index e81ea840..8e9fb6cc 100644
--- a/src/utils.h
+++ b/src/utils.h
@@ -22,6 +22,8 @@
 #include <stdlib.h>
 
 
+/* -------------------------- Fundamental constants  ------------------------ */
+
 /* Electron charge in C */
 #define ELECTRON_CHARGE (1.6021773e-19)
 
@@ -34,9 +36,16 @@
 /* Thomson scattering length (m) */
 #define THOMSON_LENGTH (2.81794e-15)
 
-/* Maxmimum index to go up to */
-#define INDMAX 40
-#define IDIM (INDMAX*2 +1)
+
+/* --------------------------- Useful datatypes ----------------------------- */
+
+struct quaternion
+{
+	double w;
+	double x;
+	double y;
+	double z;
+};
 
 
 extern unsigned int biggest(signed int a, signed int b);
@@ -59,6 +68,9 @@ extern void mapping_rotate(double x, double y, double z,
                            double omega, double tilt);
 extern void progress_bar(int val, int total);
 
+
+/* ----------------------------- Useful macros ------------------------------ */
+
 #define rad2deg(a) ((a)*180/M_PI)
 #define deg2rad(a) ((a)*M_PI/180)
 
@@ -77,51 +89,27 @@ extern void progress_bar(int val, int total);
 #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;
-
-	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;
+/* -------------------- Indexed lists for specified types ------------------- */
 
-	if ( (abs(h) > INDMAX) || (abs(k) > INDMAX) || (abs(l) > INDMAX) ) {
-		printf("\nReflection %i %i %i is out of range!\n", h, k, l);
-		printf("You need to re-configure INDMAX, delete the reflection"
-		       " cache file and re-run.\n");
-		exit(1);
-	}
-
-	if ( h < 0 ) h += IDIM;
-	if ( k < 0 ) k += IDIM;
-	if ( l < 0 ) l += IDIM;
+/* Maxmimum index to hold values up to (can be increased if necessary) */
+#define INDMAX 40
 
-	idx = h + (IDIM*k) + (IDIM*IDIM*l);
-	return ref[idx];
-}
+/* Array size */
+#define IDIM (INDMAX*2 +1)
 
+/* Create functions for storing reflection intensities indexed as h,k,l */
+#define LABEL(x) x##_intensity
+#define TYPE double
+#include "list_tmp.h"
 
-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;
-}
+/* As above, but for complex structure factors */
+#define LABEL(x) x##_sfac
+#define TYPE double complex
+#include "list_tmp.h"
 
+/* As above, but for (unsigned) integer counts */
+#define LABEL(x) x##_count
+#define TYPE unsigned int
+#include "list_tmp.h"
 
 #endif	/* UTILS_H */