compare_hkl: Show Rmerge as well

1 files changed, 63 insertions, 0 deletions

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;
+}