diff options
| -rw-r--r-- | src/compare_hkl.c | 95 |
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); } |