aboutsummaryrefslogtreecommitdiff
diff options
context:
space:
mode:
-rw-r--r--src/compare_hkl.c87
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 ) {