compare_hkl: Add Rinternal(I)

3 files changed, 50 insertions, 2 deletions

diff --git a/src/compare_hkl.c b/src/compare_hkl.c
index 959f1fed..ed8636fb 100644
--- a/src/compare_hkl.c
+++ b/src/compare_hkl.c
@@ -293,6 +293,7 @@ int main(int argc, char *argv[])
 	char *bfile = NULL;
 	char *sym = NULL;
 	double scale, scale_r2, scale_rdig, R1, R2, R1i, Rdiff, pearson;
+	double scale_rintint;
 	int i, ncom;
 	ReflItemList *i1, *i2, *icommon;
 	int config_shells = 0;
@@ -451,11 +452,16 @@ int main(int argc, char *argv[])
 	STATUS("R1(I) = %5.4f %% (scale=%5.2e)\n", R1i*100.0, scale);
 
 	Rdiff = stat_rdiff_ignore(ref1, ref2_transformed, icommon, &scale_rdig);
-	STATUS("Rdiff(F) = %5.4f %% (scale=%5.2e) (ignoring negative intensities)\n",
+	STATUS("Rint(F) = %5.4f %% (scale=%5.2e) (ignoring negative intensities)\n",
 	       Rdiff*100.0, scale);
 
 	Rdiff = stat_rdiff_zero(ref1, ref2_transformed, icommon, &scale);
-	STATUS("Rdiff(F) = %5.4f %% (scale=%5.2e) (zeroing negative intensities)\n",
+	STATUS("Rint(F) = %5.4f %% (scale=%5.2e) (zeroing negative intensities)\n",
+	       Rdiff*100.0, scale);
+
+	Rdiff = stat_rdiff_intensity(ref1, ref2_transformed, icommon,
+	                             &scale_rintint);
+	STATUS("Rint(I) = %5.4f %% (scale=%5.2e)\n",
 	       Rdiff*100.0, scale);
 
 	pearson = stat_pearson_i(ref1, ref2_transformed, icommon);
diff --git a/src/statistics.c b/src/statistics.c
index 78b3f13e..d0d3ae85 100644
--- a/src/statistics.c
+++ b/src/statistics.c
@@ -38,6 +38,7 @@ enum {
 	R_1_I,
 	R_DIFF_ZERO,
 	R_DIFF_IGNORE,
+	R_DIFF_INTENSITY,
 };
 
 
@@ -232,6 +233,35 @@ static double internal_r_i(const double *ref1, const double *ref2,
 }
 
 
+static double internal_rdiff_intensity(const double *ref1, const double *ref2,
+                                       ReflItemList *items, double scale)
+{
+	double top = 0.0;
+	double bot = 0.0;
+	int i;
+
+	for ( i=0; i<num_items(items); i++ ) {
+
+		double i1, i2;
+		struct refl_item *it;
+		signed int h, k, l;
+
+		it = get_item(items, i);
+		h = it->h;  k = it->k;  l = it->l;
+
+		i1 = lookup_intensity(ref1, h, k, l);
+		i2 = lookup_intensity(ref2, h, k, l);
+		i2 *= scale;
+
+		top += fabs(i1 - i2);
+		bot += i1 + i2;
+
+	}
+
+	return 2.0*top/bot;
+}
+
+
 static double internal_rdiff_negstozero(const double *ref1, const double *ref2,
                                         ReflItemList *items, double scale)
 {
@@ -319,6 +349,9 @@ static double calc_r(double scale, void *params)
 	case R_DIFF_IGNORE :
 		return internal_rdiff_ignorenegs(rp->ref1, rp->ref2,
 		                                 rp->items, scale);
+	case R_DIFF_INTENSITY :
+		return internal_rdiff_intensity(rp->ref1, rp->ref2,
+		                                rp->items, scale);
 	}
 
 	ERROR("No such FoM!\n");
@@ -356,6 +389,7 @@ static double r_minimised(const double *ref1, const double *ref2,
 		break;
 	case R_2 :
 	case R_1_I :
+	case R_DIFF_INTENSITY :
 		scale = stat_scale_intensity(ref1, ref2, items);
 		break;
 	}
@@ -437,6 +471,12 @@ double stat_rdiff_ignore(const double *ref1, const double *ref2,
 }
 
 
+double stat_rdiff_intensity(const double *ref1, const double *ref2,
+                            ReflItemList *items, double *scalep)
+{
+	return r_minimised(ref1, ref2, items, scalep, R_DIFF_INTENSITY);
+}
+
 double stat_pearson_i(const double *ref1, const double *ref2,
                       ReflItemList *items)
 {
diff --git a/src/statistics.h b/src/statistics.h
index 6a8efb50..18603357 100644
--- a/src/statistics.h
+++ b/src/statistics.h
@@ -37,6 +37,8 @@ extern double stat_rdiff_zero(const double *ref1, const double *ref2,
                          ReflItemList *items, double *scalep);
 extern double stat_rdiff_ignore(const double *ref1, const double *ref2,
                          ReflItemList *items, double *scalep);
+extern double stat_rdiff_intensity(const double *ref1, const double *ref2,
+                         ReflItemList *items, double *scalep);
 
 extern double stat_pearson_i(const double *ref1, const double *ref2,
                              ReflItemList *items);