From 5f053ff4194f65bc0b760e409dd8ac18abe2aee0 Mon Sep 17 00:00:00 2001 From: Thomas White Date: Mon, 15 Nov 2010 15:43:12 +0100 Subject: compare_hkl: Take resolution limits on command line --- src/compare_hkl.c | 87 ++++++++++++++++++++++++++----------------------------- 1 file 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= The symmetry of both the input files.\n" " -p, --pdb= PDB file to use (default: molecule.pdb).\n" " --shells Plot the figures of merit by resolution.\n" +" --rmin= Fix lower resolution limit for --shells (m^-1).\n" +" --rmax= 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; ih; 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; ih; 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 ) { -- cgit v1.2.3