compare_hkl: Add Luzzati plot

1 files changed, 92 insertions, 3 deletions

diff --git a/src/compare_hkl.c b/src/compare_hkl.c
index b58e6d25..6f033400 100644
--- a/src/compare_hkl.c
+++ b/src/compare_hkl.c
@@ -28,6 +28,10 @@
 #include "symmetry.h"
 
 
+/* Number of bins for Luzzati plot */
+#define NBINS (10)
+
+
 static void show_help(const char *s)
 {
 	printf("Syntax: %s [options] <file1.hkl> <file2.hkl>\n\n", s);
@@ -41,6 +45,85 @@ static void show_help(const char *s)
 }
 
 
+static void plot_luzzati(const double *ref1, const double *ref2,
+                         ReflItemList *items, double scale, UnitCell *cell)
+{
+	double num[NBINS];
+	double den[NBINS];
+	double rmin, rmax, rstep;
+	int i;
+	FILE *fh;
+
+	if ( cell == NULL ) {
+		ERROR("Need the unit cell to plot the Luzzati plot.\n");
+		return;
+	}
+
+	fh = fopen("luzzati.dat", "w");
+	if ( fh == NULL ) {
+		ERROR("Couldn't open 'luzzati.dat'\n");
+		return;
+	}
+
+	for ( i=0; i<NBINS; i++ ) {
+		num[i] = 0.0;
+		den[i] = 0.0;
+	}
+
+	rmin = +INFINITY;
+	rmax = 0.0;
+	for ( i=0; i<num_items(items); i++ ) {
+
+		struct refl_item *it;
+		signed int h, k, l;
+		double res;
+
+		it = get_item(items, i);
+		h = it->h;  k = it->k;  l = it->l;
+
+		res = 2.0*resolution(cell, h, k, l);
+		if ( res > rmax ) rmax = res;
+		if ( res < rmin ) rmin = res;
+
+	}
+	rstep = (rmax-rmin) / NBINS;
+
+	for ( i=0; i<num_items(items); i++ ) {
+
+		struct refl_item *it;
+		signed int h, k, l;
+		double res;
+		int bin;
+		double i1, i2;
+
+		it = get_item(items, i);
+		h = it->h;  k = it->k;  l = it->l;
+
+		res = 2.0*resolution(cell, h, k, l);
+
+		bin = (res-rmin)/rstep;
+
+		i1 = lookup_intensity(ref1, h, k, l);
+		i2 = scale * lookup_intensity(ref2, h, k, l);
+
+		num[bin] += pow(i1 - i2, 2.0);
+		den[bin] += pow(i1, 2.0);
+
+	}
+
+	for ( i=0; i<NBINS; i++ ) {
+
+		double r2, cen;
+		cen = rmin + rstep*i + rstep/2.0;
+		r2 = sqrt(num[i]/den[i]);
+		fprintf(fh, "%f %f\n", cen/1.0e9, r2*100.0);
+
+	}
+
+	fclose(fh);
+}
+
+
 int main(int argc, char *argv[])
 {
 	int c;
@@ -53,15 +136,17 @@ int main(int argc, char *argv[])
 	char *afile = NULL;
 	char *bfile = NULL;
 	char *sym = NULL;
-	double scale, R1, R2, Rdiff, pearson;
+	double scale, scale_r2, R1, R2, Rdiff, pearson;
 	int i, ncom;
 	ReflItemList *i1, *i2, *icommon;
+	int config_luzzati = 0;
 
 	/* Long options */
 	const struct option longopts[] = {
 		{"help",               0, NULL,               'h'},
 		{"output",             1, NULL,               'o'},
 		{"symmetry",           1, NULL,               'y'},
+		{"luzzati",            0, &config_luzzati,     1},
 		{0, 0, NULL, 0}
 	};
 
@@ -149,13 +234,17 @@ int main(int argc, char *argv[])
 	       num_items(i1), num_items(i2), ncom);
 	R1 = stat_r1(ref1, ref2_transformed, icommon, &scale);
 	STATUS("R1 = %5.4f %% (scale=%5.2e)\n", R1*100.0, scale);
-	R2 = stat_r2(ref1, ref2_transformed, icommon, &scale);
-	STATUS("R2 = %5.4f %% (scale=%5.2e)\n", R2*100.0, scale);
+	R2 = stat_r2(ref1, ref2_transformed, icommon, &scale_r2);
+	STATUS("R2 = %5.4f %% (scale=%5.2e)\n", R2*100.0, scale_r2);
 	Rdiff = stat_rdiff(ref1, ref2_transformed, icommon, &scale);
 	STATUS("Rdiff = %5.4f %% (scale=%5.2e)\n", Rdiff*100.0, scale);
 	pearson = stat_pearson(ref1, ref2_transformed, icommon);
 	STATUS("Pearson r = %5.4f\n", pearson);
 
+	if ( config_luzzati ) {
+		plot_luzzati(ref1, ref2_transformed, icommon, scale_r2, cell);
+	}
+
 	if ( outfile != NULL ) {
 		write_reflections(outfile, icommon, out, NULL, NULL, cell);
 	}