compare_hkl: Take resolution limits on command line

1 files changed, 41 insertions, 46 deletions

diff --git a/src/compare_hkl.c b/src/compare_hkl.c
index 2bdbd5a4..a0481559 100644
--- a/src/compare_hkl.c
+++ b/src/compare_hkl.c
@@ -43,6 +43,8 @@ static void show_help(const char *s)
 "  -y, --symmetry=<sym>       The symmetry of both the input files.\n"
 "  -p, --pdb=<filename>       PDB file to use (default: molecule.pdb).\n"
 "      --shells               Plot the figures of merit by resolution.\n"
+"      --rmin=<res>           Fix lower resolution limit for --shells (m^-1).\n"
+"      --rmax=<res>           Fix upper resolution limit for --shells (m^-1).\n"
 "\n");
 }
 
@@ -50,7 +52,8 @@ static void show_help(const char *s)
 static void plot_shells(const double *ref1, const double *ref2,
                         ReflItemList *items, double scale, UnitCell *cell,
                         const char *sym, ReflItemList *characterise,
-                        unsigned int *char_counts, const double *sigma)
+                        unsigned int *char_counts, const double *sigma,
+                        double rmin_fix, double rmax_fix)
 {
 	double num[NBINS];
 	int cts[NBINS];
@@ -71,6 +74,7 @@ static void plot_shells(const double *ref1, const double *ref2,
 	double snr_total = 0;
 	int nmeas = 0;
 	int nmeastot = 0;
+	int nout = 0;
 
 	if ( cell == NULL ) {
 		ERROR("Need the unit cell to plot resolution shells.\n");
@@ -109,14 +113,15 @@ static void plot_shells(const double *ref1, const double *ref2,
 
 	}
 
-	STATUS("%f -> %f\n", rmin/1e9, rmax/1e9);
+	STATUS("1/d goes from %f to %f nm^-1\n", rmin/1e9, rmax/1e9);
 
-	/* Increase the max just a little bit */
+	/* Widen the range just a little bit */
+	rmin -= 0.001e9;
 	rmax += 0.001e9;
 
-	/* FIXME: Fixed resolution shells */
-	rmin = 0.120e9;
-	rmax = 1.172e9;
+	/* Fixed resolution shells if needed */
+	if ( rmin_fix > 0.0 ) rmin = rmin_fix;
+	if ( rmax_fix > 0.0 ) rmax = rmax_fix;
 
 	total_vol = pow(rmax, 3.0) - pow(rmin, 3.0);
 	vol_per_shell = total_vol / NBINS;
@@ -144,24 +149,6 @@ static void plot_shells(const double *ref1, const double *ref2,
 	STATUS("Shell %i: %f to %f\n", NBINS-1,
 	       rmins[NBINS-1]/1e9, rmaxs[NBINS-1]/1e9);
 
-#if 0
-	/* FIXME: Fixed resolution shells */
-	rmins[0] = 0.121065;  rmaxs[0] = 0.552486;
-	rmins[1] = 0.552486;  rmaxs[1] = 0.690186;
-	rmins[2] = 0.690186;  rmaxs[2] = 0.787787;
-	rmins[3] = 0.787787;  rmaxs[3] = 0.865813;
-	rmins[4] = 0.865813;  rmaxs[4] = 0.931853;
-	rmins[5] = 0.931853;  rmaxs[5] = 0.989663;
-	rmins[6] = 0.989663;  rmaxs[6] = 1.041409;
-	rmins[7] = 1.041409;  rmaxs[7] = 1.088467;
-	rmins[8] = 1.088467;  rmaxs[8] = 1.131775;
-	rmins[9] = 1.131775;  rmaxs[9] = 1.172000;
-	for ( i=0; i<NBINS; i++ ) {
-		rmins[i] *= 1e9;
-		rmaxs[i] *= 1e9;
-	}
-#endif
-
 	/* Count the number of reflections possible in each shell */
 	counted = new_list_count();
 	for ( h=-50; h<=+50; h++ ) {
@@ -172,9 +159,6 @@ static void plot_shells(const double *ref1, const double *ref2,
 		signed int hs, ks, ls;
 		int bin;
 
-		/* FIXME: Reflection condition */
-		if ( (h==0) && (k==0) && (l%2) ) continue;
-
 		get_asymm(h, k, l, &hs, &ks, &ls, sym);
 		if ( lookup_count(counted, hs, ks, ls) ) continue;
 		set_count(counted, hs, ks, ls, 1);
@@ -209,9 +193,6 @@ static void plot_shells(const double *ref1, const double *ref2,
 		it = get_item(characterise, i);
 		h = it->h;  k = it->k;  l = it->l;
 
-		/* FIXME: Reflection condition */
-		if ( (h==0) && (k==0) && (l%2) ) continue;
-
 		d = resolution(cell, h, k, l) * 2.0;
 
 		bin = -1;
@@ -222,7 +203,7 @@ static void plot_shells(const double *ref1, const double *ref2,
 			}
 		}
 		if ( bin == -1 ) {
-			ERROR("Warnung! %i %i %i %f\n", h, k, l, d/1e9);
+			nout++;
 			continue;
 		}
 
@@ -240,6 +221,11 @@ static void plot_shells(const double *ref1, const double *ref2,
 	STATUS("%i measurements in total.\n", nmeastot);
 	STATUS("%i reflections in total.\n", nmeas);
 
+	if ( nout ) {
+		STATUS("Warning; %i reflections outside resolution range.\n",
+		       nout);
+	}
+
 	den = 0.0;
 	ctot = 0;
 	for ( i=0; i<num_items(items); i++ ) {
@@ -248,15 +234,12 @@ static void plot_shells(const double *ref1, const double *ref2,
 		signed int h, k, l;
 		double d;
 		int bin;
-		double i1, i2, f1, f2;
+		double i1, i2;
 		int j;
 
 		it = get_item(items, i);
 		h = it->h;  k = it->k;  l = it->l;
 
-		/* FIXME: Reflection condition */
-		if ( (h==0) && (k==0) && (l%2) ) continue;
-
 		d = resolution(cell, h, k, l) * 2.0;
 
 		bin = -1;
@@ -266,22 +249,16 @@ static void plot_shells(const double *ref1, const double *ref2,
 				break;
 			}
 		}
-		if ( bin == -1 ) {
-			ERROR("Warnung! %i %i %i %f\n", h, k, l, d/1e9);
-			abort();
-			continue;
-		}
+
+		/* Outside resolution range? */
+		if ( bin == -1 ) continue;
 
 		i1 = lookup_intensity(ref1, h, k, l);
-		//if ( i1 < 0.0 ) continue;
-		//f1 = sqrt(i1);
 		i2 = lookup_intensity(ref2, h, k, l);
-		//if ( i2 < 0.0 ) continue;
-		//f2 = sqrt(i2);
 		i2 *= scale;
 
 		num[bin] += fabs(i1 - i2);
-		den += i1;// + i2) / 2.0;
+		den += i1;
 		ctot++;
 		cts[bin]++;
 
@@ -327,6 +304,8 @@ int main(int argc, char *argv[])
 	int rej1 = 0;
 	int rej2 = 0;
 	unsigned int *cts1;
+	float rmin_fix = -1.0;
+	float rmax_fix = -1.0;
 
 	/* Long options */
 	const struct option longopts[] = {
@@ -335,6 +314,8 @@ int main(int argc, char *argv[])
 		{"symmetry",           1, NULL,               'y'},
 		{"shells",             0, &config_shells,     1},
 		{"pdb",                1, NULL,               'p'},
+		{"rmin",               1, NULL,               2},
+		{"rmax",               1, NULL,               3},
 		{0, 0, NULL, 0}
 	};
 
@@ -361,6 +342,20 @@ int main(int argc, char *argv[])
 		case 0 :
 			break;
 
+		case 2 :
+			if ( sscanf(optarg, "%e", &rmin_fix) != 1 ) {
+				ERROR("Invalid value for --rmin\n");
+				return 1;
+			}
+			break;
+
+		case 3 :
+			if ( sscanf(optarg, "%e", &rmax_fix) != 1 ) {
+				ERROR("Invalid value for --rmax\n");
+				return 1;
+			}
+			break;
+
 		default :
 			return 1;
 		}
@@ -501,7 +496,7 @@ int main(int argc, char *argv[])
 
 	if ( config_shells ) {
 		plot_shells(ref1, ref2_transformed, icommon, scale_r1fi,
-		            cell, sym, i1, cts1, esd1);
+		            cell, sym, i1, cts1, esd1, rmin_fix, rmax_fix);
 	}
 
 	if ( outfile != NULL ) {