Add compare_hkl utility for working out Lorentz factors

4 files changed, 173 insertions, 1 deletions

diff --git a/src/Makefile.am b/src/Makefile.am
index 824211f8..19e8ca41 100644
--- a/src/Makefile.am
+++ b/src/Makefile.am
@@ -1,4 +1,4 @@
-bin_PROGRAMS = pattern_sim process_hkl get_hkl indexamajig
+bin_PROGRAMS = pattern_sim process_hkl get_hkl indexamajig compare_hkl
 
 if HAVE_GTK
 bin_PROGRAMS += hdfsee
@@ -33,3 +33,6 @@ endif
 
 get_hkl_SOURCES = get_hkl.c sfac.c cell.c utils.c reflections.c
 get_hkl_LDADD = @LIBS@
+
+compare_hkl_SOURCES = compare_hkl.c sfac.c cell.c utils.c reflections.c
+compare_hkl_LDADD = @LIBS@
diff --git a/src/compare_hkl.c b/src/compare_hkl.c
new file mode 100644
index 00000000..00b8bfc1
--- /dev/null
+++ b/src/compare_hkl.c
@@ -0,0 +1,133 @@
+/*
+ * compare_hkl.c
+ *
+ * Compare reflection lists
+ *
+ * (c) 2006-2010 Thomas White <taw@physics.org>
+ *
+ * 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 <getopt.h>
+
+#include "utils.h"
+#include "sfac.h"
+#include "reflections.h"
+
+
+static void show_help(const char *s)
+{
+	printf("Syntax: %s [options] <file1.hkl> <file2.hkl>\n\n", s);
+	printf(
+"Compare intensity lists.\n"
+"\n"
+"  -h, --help                 Display this help message.\n"
+"  -o, --output=<filename>    Specify output filename for correction factor.\n"
+"\n");
+}
+
+
+int main(int argc, char *argv[])
+{
+	int c;
+	double *ref1;
+	double *ref2;
+	double *out;
+	struct molecule *mol;
+	char *outfile = NULL;
+	char *afile = NULL;
+	char *bfile = NULL;
+	signed int h, k, l;
+
+	/* Long options */
+	const struct option longopts[] = {
+		{"help",               0, NULL,               'h'},
+		{"output",             1, NULL,               'o'},
+		{0, 0, NULL, 0}
+	};
+
+	/* Short options */
+	while ((c = getopt_long(argc, argv, "ho:a:b:", longopts, NULL)) != -1) {
+
+		switch (c) {
+		case 'h' : {
+			show_help(argv[0]);
+			return 0;
+		}
+
+		case 'o' : {
+			outfile = strdup(optarg);
+			break;
+		}
+
+		case 'a' : {
+			afile = strdup(optarg);
+			break;
+		}
+
+		case 'b' : {
+			bfile = strdup(optarg);
+			break;
+		}
+
+		case 0 : {
+			break;
+		}
+
+		default : {
+			return 1;
+		}
+		}
+
+	}
+
+	if ( outfile == NULL ) {
+		ERROR("You must specify the output filename with -o\n");
+		return 1;
+	}
+
+	mol = load_molecule();
+	ref1 = read_reflections(afile);
+	if ( ref1 == NULL ) {
+		ERROR("Couldn't open file '%s'\n", afile);
+		return 1;
+	}
+	ref2 = read_reflections(bfile);
+	if ( ref2 == NULL ) {
+		ERROR("Couldn't open file '%s'\n", bfile);
+		return 1;
+	}
+	out = new_list_intensity();
+
+	for ( h=-INDMAX; h<INDMAX; h++ ) {
+	for ( k=-INDMAX; k<INDMAX; k++ ) {
+	for ( l=-INDMAX; l<INDMAX; l++ ) {
+
+		double i1, i2;
+
+		i1 = lookup_intensity(ref1, h, k, l);
+		i2 = lookup_intensity(ref2, h, k, l);
+
+		if ( (i1 != 0.0) && (i2 != 0.0) ) {
+			set_intensity(out, h, k, l, i1/i2);
+		}
+
+	}
+	}
+	}
+
+	write_reflections(outfile, NULL, out, 1, mol->cell);
+
+	return 0;
+}
diff --git a/src/reflections.c b/src/reflections.c
index 88b2ab59..ab7738e4 100644
--- a/src/reflections.c
+++ b/src/reflections.c
@@ -77,6 +77,40 @@ void write_reflections(const char *filename, unsigned int *counts,
 }
 
 
+double *read_reflections(const char *filename)
+{
+	double *ref;
+	FILE *fh;
+	char *rval;
+
+	fh = fopen(filename, "r");
+	if ( fh == NULL ) {
+		ERROR("Failed to open input file\n");
+		return NULL;
+	}
+
+	ref = new_list_intensity();
+
+	do {
+
+		char line[1024];
+		signed int h, k, l, intensity;
+		int r;
+
+		rval = fgets(line, 1023, fh);
+		r = sscanf(line, "%i %i %i %i", &h, &k, &l, &intensity);
+		if ( r != 4 ) continue;
+
+		set_intensity(ref, h, k, l, intensity);
+
+	} while ( rval != NULL );
+
+	fclose(fh);
+
+	return ref;
+}
+
+
 double *ideal_intensities(double complex *sfac)
 {
 	double *ref;
diff --git a/src/reflections.h b/src/reflections.h
index 7bc426ca..14209c34 100644
--- a/src/reflections.h
+++ b/src/reflections.h
@@ -23,6 +23,8 @@
 extern void write_reflections(const char *filename, unsigned int *counts,
                               double *ref, int zone_axis, UnitCell *cell);
 
+extern double *read_reflections(const char *filename);
+
 extern double *ideal_intensities(double complex *sfac);