Add process_hkl program (replaces integr_sim)

.gitignore

@@ -10,5 +10,5 @@ stamp-h1
 src/*.o
 src/.deps
 src/pattern_sim
-src/integr_sim
+src/process_hkl
 *~

src/Makefile.am

-bin_PROGRAMS = pattern_sim integr_sim
+bin_PROGRAMS = pattern_sim process_hkl
 
 AM_CFLAGS = -Wall
 
@@ -6,6 +6,5 @@ pattern_sim_SOURCES = pattern_sim.c diffraction.c utils.c image.c cell.c \
                       hdf5-file.c ewald.c detector.c sfac.c intensities.c
 pattern_sim_LDADD = @LIBS@
 
-integr_sim_SOURCES = integr_sim.c diffraction.c utils.c ewald.c sfac.c \
-                     statistics.c cell.c
-integr_sim_LDADD = @LIBS@
+process_hkl_SOURCES = process_hkl.c sfac.c statistics.c cell.c utils.c
+process_hkl_LDADD = @LIBS@

src/integr_sim.c → src/process_hkl.c

 /*
- * integr_sim.c
+ * process_hkl.c
  *
- * Test integration of intensities
+ * Assemble and process FEL Bragg intensities
  *
  * (c) 2006-2009 Thomas White <thomas.white@desy.de>
  *
@@ -19,19 +19,19 @@
 #include <stdio.h>
 #include <string.h>
 #include <unistd.h>
+#include <getopt.h>
 
-#include "cell.h"
-#include "image.h"
 #include "utils.h"
 #include "statistics.h"
 #include "sfac.h"
 
 
-static void main_show_help(const char *s)
+static void show_help(const char *s)
 {
-	printf("Syntax: %s\n\n", s);
-	printf("Test relrod integration\n\n");
-	printf("  -h              Display this help message\n");
+	printf("Syntax: %s [options]\n\n", s);
+	printf("Assemble and process FEL Bragg intensities.\n\n");
+	printf("  -h              Display this help message.\n");
+	printf("  -i <filename>   Specify input filename (\"-\" for stdin).\n");
 }
 
 
@@ -94,139 +94,175 @@ static void write_RvsQ(const char *name, double *ref, double *trueref,
 }
 
 
-int main(int argc, char *argv[])
+static void write_reflections(const char *filename, unsigned int *counts,
+                              double *ref)
 {
-	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;
+	signed int h, k, l;
 
-	while ((c = getopt(argc, argv, "h")) != -1) {
+	fh = fopen(filename, "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);
+}
 
-		switch ( c ) {
 
-			case 'h' : {
-				main_show_help(argv[0]);
-				return 0;
-			}
+static double *ideal_intensities(double complex *sfac)
+{
+	double *ref;
+	signed int h, k, l;
 
-		}
+	ref = new_list_intensity();
 
+	/* Generate ideal reflections from complex structure factors */
+	for ( h=-INDMAX; h<=INDMAX; h++ ) {
+	for ( k=-INDMAX; k<=INDMAX; k++ ) {
+	for ( l=-INDMAX; l<=INDMAX; l++ ) {
+		double complex F = lookup_sfac(sfac, h, k, l);
+		double intensity = pow(cabs(F), 2.0);
+		set_intensity(ref, h, k, l, intensity);
+	}
+	}
 	}
 
-	/* 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();
+	return ref;
+}
 
-	/* 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 */
+static void process_reflections(double *ref, double *trueref,
+                                unsigned int *counts, unsigned int n_patterns,
+                                UnitCell *cell)
+{
+	int j;
+	double mean_counts;
+	int ctot = 0;
+	int nmeas = 0;
+	double ff;
+	char name[64];
+	double R, scale;
 
-	for ( i=1; i<=10e3; i++ ) {
+	for ( j=0; j<LIST_SIZE; j++ ) {
+		ctot += counts[j];
+		if ( counts[j] > 0 ) nmeas++;
+	}
+	mean_counts = (double)ctot/nmeas;
 
-		#if 0
-		image.orientation = random_quaternion();
+	ff = lookup_intensity(ref, 2, 2, 2) / lookup_count(counts, 2, 2, 2);
 
-		/* Calculate reflections using large smax
-		 * (rather than the actual value) */
-		//get_reflections(&image, cell, 1.0e9);
+	R = stat_r2(ref, trueref, counts, LIST_SIZE, &scale);
+	STATUS("%8u: R=%5.2f%% (sf=%7.4e)"
+	       " mean meas/refl=%5.2f,"
+	       " %i reflections measured. %f\n",
+	       n_patterns, R*100.0, scale, mean_counts, nmeas, ff);
 
-		nrefl = image_feature_count(image.rflist);
-		for ( j=0; j<nrefl; j++ ) {
+	/* Record graph of R against q for this N */
+	snprintf(name, 63, "results/R_vs_q-%u.dat", n_patterns);
+	write_RvsQ(name, ref, trueref, counts, scale, cell);
+}
 
-			struct imagefeature *f;
-			double t, s, intensity, F;
 
-			f = image_get_feature(image.rflist, j);
+int main(int argc, char *argv[])
+{
+	int c;
+	char *filename = NULL;
+	FILE *fh;
+	unsigned int n_patterns;
+	double *ref, *trueref;
+	unsigned int *counts;
+	char *rval;
+	struct molecule *mol;
+
+	/* Long options */
+	const struct option longopts[] = {
+		{"help",               0, NULL,               'h'},
+		{"input",              1, NULL,               'i'},
+		{0, 0, NULL, 0}
+	};
+
+	/* Short options */
+	while ((c = getopt_long(argc, argv, "hi:", longopts, NULL)) != -1) {
+
+		switch (c) {
+		case 'h' : {
+			show_help(argv[0]);
+			return 0;
+		}
 
-			t = 100.0e-9;  /* Thickness 100 nm */
-			s = f->s;      /* Get excitation error */
-			F = structure_factor(f->h, f->k, f->l);
+		case 'i' : {
+			filename = strdup(optarg);
+			break;
+		}
 
-			/* Calculate intensity from this reflection */
-			intensity = pow( F * SINC(M_PI*t*s), 2.0);
+		case 0 : {
+			break;
+		}
 
-			if ( intensity < 0.1 ) continue;
+		default : {
+			return 1;
+		}
+		}
 
-			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);
-			}
+	if ( filename == NULL ) {
+		ERROR("Please specify filename using the -i option\n");
+		return 1;
+	}
 
-			integrate_reflection(ref, f->h, f->k, f->l, intensity);
-			add_count(counts, f->h, f->k, f->l, 1);
-		}
-		#endif
+	mol = load_molecule();
+	get_reflections_cached(mol, eV_to_J(2.0e3));
 
-		if ( i % 1000 == 0 ) {
+	ref = new_list_intensity();
+	counts = new_list_count();
+	trueref = ideal_intensities(mol->reflections);
 
-			int j;
-			double mean_counts;
-			int ctot = 0;
-			int nmeas = 0;
-			double ff;
-			char name[64];
+	if ( strcmp(filename, "-") == 0 ) {
+		fh = stdin;
+	} else {
+		fh = fopen(filename, "r");
+	}
+	free(filename);
+	if ( fh == NULL ) {
+		ERROR("Failed to open input file\n");
+		return 1;
+	}
 
-			for ( j=0; j<ref_size; j++ ) {
-				ctot += counts[j];
-				if ( counts[j] > 0 ) nmeas++;
-			}
-			mean_counts = (double)ctot/nmeas;
+	n_patterns = 0;
+	do {
 
-			ff = lookup_intensity(ref, 2, 2, 2)
-			                        / lookup_count(counts, 2, 2, 2);
+		char line[1024];
+		int h, k, l, intensity;
+		int r;
 
-			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);
+		rval = fgets(line, 1023, fh);
+		if ( strncmp(line, "New pattern", 11) == 0 ) {
+			n_patterns++;
+			if ( n_patterns % 1000 == 0 ) {
+				process_reflections(ref, trueref, counts,
+				                    n_patterns, mol->cell);
+			}
+		}
 
-			/* 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);
+		r = sscanf(line, "%i %i %i %i", &h, &k, &l, &intensity);
+		if ( r != 4 ) continue;
 
-		}
+		integrate_intensity(ref, h, k, l, intensity);
+		integrate_count(counts, h, k, l, 1);
 
-	}
+	} while ( rval != NULL );
 
-	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);
 
+	write_reflections("results/reflections.hkl", counts, ref);
+
 	return 0;
 }

src/utils.h

@@ -131,6 +131,8 @@ static inline double distance3d(double x1, double y1, double z1,
 /* Array size */
 #define IDIM (INDMAX*2 +1)
 
+#define LIST_SIZE (IDIM*IDIM*IDIM)
+
 /* Create functions for storing reflection intensities indexed as h,k,l */
 #define LABEL(x) x##_intensity
 #define TYPE double