compare_hkl: Show Rmerge as well

3 files changed, 72 insertions, 3 deletions

diff --git a/src/compare_hkl.c b/src/compare_hkl.c
index b14a8dab..b2b40da1 100644
--- a/src/compare_hkl.c
+++ b/src/compare_hkl.c
@@ -50,7 +50,7 @@ int main(int argc, char *argv[])
 	char *afile = NULL;
 	char *bfile = NULL;
 	signed int h, k, l;
-	double scale, R;
+	double scale, R2, Rmerge;
 	unsigned int *c1;
 	unsigned int *c2;
 	int i;
@@ -141,8 +141,10 @@ int main(int argc, char *argv[])
 		ncom += c1[i] && c2[i];
 	}
 	STATUS("%i,%i reflections: %i in common\n", nc1, nc2, ncom);
-	R = stat_r2(ref1, c1, ref2, c2, &scale);
-	STATUS("R2 = %5.4f %% (scale=%5.2f)\n", R*100.0, scale);
+	R2 = stat_r2(ref1, c1, ref2, c2, &scale);
+	STATUS("R2 = %5.4f %% (scale=%5.2f)\n", R2*100.0, scale);
+	Rmerge = stat_rmerge(ref1, c1, ref2, c2, &scale);
+	STATUS("Rmerge = %5.4f %% (scale=%5.2f)\n", Rmerge*100.0, scale);
 
 	if ( outfile != NULL ) {
 		write_reflections(outfile, NULL, out, NULL, 1, cell, 1);
diff --git a/src/statistics.c b/src/statistics.c
index b1589ebe..d78dde9e 100644
--- a/src/statistics.c
+++ b/src/statistics.c
@@ -73,3 +73,66 @@ double stat_r2(const double *ref1, const unsigned int *c1,
 
 	return sqrt(top/bot);
 }
+
+
+double stat_scale_sqrti(const double *ref1, const unsigned int *c1,
+                            const double *ref2, const unsigned int *c2)
+{
+	double top = 0.0;
+	double bot = 0.0;
+	int i;
+
+	/* Start from i=1 -> skip central beam */
+	for ( i=1; i<LIST_SIZE; i++ ) {
+
+		if ( c1[i] && c2[i] ) {
+
+			double f1, f2;
+
+			if ( (ref1[i]<0.0) || (ref2[i]<0.0) ) continue;
+
+			f1 = sqrt(ref1[i]) / (double)c1[i];
+			f2 = sqrt(ref2[i]) / (double)c2[i];
+
+			top += f1 * f2;
+			bot += f2 * f2;
+
+		} /* else reflection not common so don't worry about it */
+
+	}
+
+	return top/bot;
+}
+
+
+double stat_rmerge(const double *ref1, const unsigned int *c1,
+                   const double *ref2, const unsigned int *c2, double *scalep)
+{
+	double top = 0.0;
+	double bot = 0.0;
+	double scale;
+	int i;
+	scale = stat_scale_sqrti(ref1, c1, ref2, c2);
+	*scalep = scale;
+
+	/* Start from i=1 -> skip central beam */
+	for ( i=1; i<LIST_SIZE; i++ ) {
+
+		if ( c1[i] && c2[i] ) {
+
+			double f1, f2;
+
+			if ( (ref1[i]<0.0) || (ref2[i]<0.0) ) continue;
+
+			f1 = sqrt(ref1[i]) / (scale*(double)c1[i]);
+			f2 = sqrt(ref2[i]) / (double)c2[i];
+
+			top += fabs(f1 - f2);
+			bot += f1 + f2;
+
+		} /* else reflection not measured so don't worry about it */
+
+	}
+
+	return 2.0*top/bot;
+}
diff --git a/src/statistics.h b/src/statistics.h
index fe7f57e1..e19d4b29 100644
--- a/src/statistics.h
+++ b/src/statistics.h
@@ -23,4 +23,8 @@ extern double stat_r2(const double *ref1, const unsigned int *c1,
                       const double *ref2, const unsigned int *c2,
                       double *scalep);
 
+extern double stat_rmerge(const double *ref1, const unsigned int *c1,
+                          const double *ref2, const unsigned int *c2,
+                          double *scalep);
+
 #endif	/* STATISTICS_H */