compare_hkl: Display R-factor

3 files changed, 32 insertions, 5 deletions

diff --git a/src/Makefile.am b/src/Makefile.am
index 4041e1e2..fcb1f2bd 100644
--- a/src/Makefile.am
+++ b/src/Makefile.am
@@ -38,7 +38,8 @@ 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_SOURCES = compare_hkl.c sfac.c cell.c utils.c reflections.c \
+                      statistics.c
 compare_hkl_LDADD = @LIBS@
 
 powder_plot_SOURCES = powder_plot.c cell.c utils.c image.c hdf5-file.c \
diff --git a/src/Makefile.in b/src/Makefile.in
index 1d7547a0..0ccef8a1 100644
--- a/src/Makefile.in
+++ b/src/Makefile.in
@@ -52,7 +52,8 @@ CONFIG_CLEAN_VPATH_FILES =
 am__installdirs = "$(DESTDIR)$(bindir)"
 PROGRAMS = $(bin_PROGRAMS)
 am_compare_hkl_OBJECTS = compare_hkl.$(OBJEXT) sfac.$(OBJEXT) \
-	cell.$(OBJEXT) utils.$(OBJEXT) reflections.$(OBJEXT)
+	cell.$(OBJEXT) utils.$(OBJEXT) reflections.$(OBJEXT) \
+	statistics.$(OBJEXT)
 compare_hkl_OBJECTS = $(am_compare_hkl_OBJECTS)
 compare_hkl_DEPENDENCIES =
 am_get_hkl_OBJECTS = get_hkl.$(OBJEXT) sfac.$(OBJEXT) cell.$(OBJEXT) \
@@ -239,7 +240,9 @@ indexamajig_LDADD = @LIBS@
 @HAVE_GTK_TRUE@hdfsee_LDADD = @LIBS@
 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_SOURCES = compare_hkl.c sfac.c cell.c utils.c reflections.c \
+                      statistics.c
+
 compare_hkl_LDADD = @LIBS@
 powder_plot_SOURCES = powder_plot.c cell.c utils.c image.c hdf5-file.c \
                       detector.c
diff --git a/src/compare_hkl.c b/src/compare_hkl.c
index 63daf318..433459ed 100644
--- a/src/compare_hkl.c
+++ b/src/compare_hkl.c
@@ -24,6 +24,7 @@
 #include "utils.h"
 #include "sfac.h"
 #include "reflections.h"
+#include "statistics.h"
 
 
 static void show_help(const char *s)
@@ -49,6 +50,12 @@ int main(int argc, char *argv[])
 	char *afile = NULL;
 	char *bfile = NULL;
 	signed int h, k, l;
+	double scale, R;
+	unsigned int *c1;
+	unsigned int *c2;
+	unsigned int *cjoint;
+	int i;
+	int nc1, nc2, ncom;
 
 	/* Long options */
 	const struct option longopts[] = {
@@ -93,12 +100,14 @@ int main(int argc, char *argv[])
 	}
 
 	cell = load_cell_from_pdb("molecule.pdb");
-	ref1 = read_reflections(afile, NULL);
+	c1 = new_list_count();
+	ref1 = read_reflections(afile, c1);
 	if ( ref1 == NULL ) {
 		ERROR("Couldn't open file '%s'\n", afile);
 		return 1;
 	}
-	ref2 = read_reflections(bfile, NULL);
+	c2 = new_list_count();
+	ref2 = read_reflections(bfile, c2);
 	if ( ref2 == NULL ) {
 		ERROR("Couldn't open file '%s'\n", bfile);
 		return 1;
@@ -122,6 +131,20 @@ int main(int argc, char *argv[])
 	}
 	}
 
+	cjoint = new_list_count();
+	nc1 = 0;
+	nc2 = 0;
+	ncom = 0;
+	for ( i=0; i<IDIM*IDIM*IDIM; i++ ) {
+		cjoint[i] = c1[i] && c2[i];
+		nc1 += c1[i];
+		nc2 += c2[i];
+		ncom += cjoint[i];
+	}
+	STATUS("%i,%i reflections: %i in common\n", nc1, nc2, ncom);
+	R = stat_r2(ref1, ref2, cjoint, IDIM*IDIM*IDIM, &scale);
+	STATUS("R2 = %5.4f %% (scale=%5.2f)\n", R*100.0, scale);
+
 	if ( outfile != NULL ) {
 		write_reflections(outfile, NULL, out, 1, cell);
 	}