From b227f195cd8e3dfeae0212b1de7bdcb5a6fc837a Mon Sep 17 00:00:00 2001 From: Thomas White Date: Fri, 4 Jan 2013 11:41:34 +0100 Subject: compare_hkl: Add --intensity-shells --- src/compare_hkl.c | 274 +++++++++++++++++++++++++++++++++++++++++++----------- 1 file changed, 218 insertions(+), 56 deletions(-) diff --git a/src/compare_hkl.c b/src/compare_hkl.c index 99ab77f2..f7776598 100644 --- a/src/compare_hkl.c +++ b/src/compare_hkl.c @@ -315,99 +315,238 @@ static double fom_shell(struct fom_context *fctx, int i) } -static void do_fom(RefList *list1, RefList *list2, UnitCell *cell, - double rmin, double rmax, enum fom fom, - int config_unity, int nshells, const char *filename) +struct shells { - int *cts; + int config_intshells; + int nshells; double *rmins; double *rmaxs; - double total_vol, vol_per_shell; +}; + + +static struct shells *set_intensity_shells(double min_I, double max_I, + int nshells) +{ + struct shells *s; int i; - Reflection *refl1; - RefListIterator *iter; - FILE *fh; - double scale; - struct fom_context *fctx; - cts = malloc(nshells*sizeof(int)); - rmins = malloc(nshells*sizeof(double)); - rmaxs = malloc(nshells*sizeof(double)); - fctx = init_fom(fom, num_reflections(list1), nshells); + if ( min_I >= max_I ) { + ERROR("Invalid intensity range.\n"); + return NULL; + } - if ( (fctx==NULL) || (cts==NULL) || (rmins==NULL) || (rmaxs==NULL) ) - { - ERROR("Couldn't allocate memory for resolution shells.\n"); - return; + /* Adjust minimum and maximum intensities to get the most densely + * populated part of the reflections */ + max_I = min_I + (max_I-min_I)/5000.0; + + s = malloc(sizeof(struct shells)); + if ( s == NULL ) return NULL; + + s->rmins = malloc(nshells*sizeof(double)); + s->rmaxs = malloc(nshells*sizeof(double)); + + if ( (s->rmins==NULL) || (s->rmaxs==NULL) ) { + ERROR("Couldn't allocate memory for shells.\n"); + free(s); + return NULL; } + s->config_intshells = 1; + s->nshells = nshells; + for ( i=0; irmins[i] = min_I + i*(max_I - min_I)/nshells;; + s->rmaxs[i] = min_I + (i+1)*(max_I - min_I)/nshells;; + } - if ( config_unity ) { - scale = 1.0; - } else { - scale = stat_scale_intensity(list1, list2); + return s; +} + + +static struct shells *set_resolution_shells(double rmin, double rmax, + int nshells) +{ + struct shells *s; + double total_vol, vol_per_shell; + int i; + + s = malloc(sizeof(struct shells)); + if ( s == NULL ) return NULL; + + s->rmins = malloc(nshells*sizeof(double)); + s->rmaxs = malloc(nshells*sizeof(double)); + + if ( (s->rmins==NULL) || (s->rmaxs==NULL) ) { + ERROR("Couldn't allocate memory for resolution shells.\n"); + free(s); + return NULL; } - /* Calculate the resolution bins */ + s->config_intshells = 0; + s->nshells = nshells; + total_vol = pow(rmax, 3.0) - pow(rmin, 3.0); vol_per_shell = total_vol / nshells; - rmins[0] = rmin; + s->rmins[0] = rmin; for ( i=1; irmins[i-1], 3.0); r = pow(r, 1.0/3.0); /* Shells of constant volume */ - rmaxs[i-1] = r; - rmins[i] = r; + s->rmaxs[i-1] = r; + s->rmins[i] = r; /* Shells of constant thickness */ //rmins[i] = rmins[i-1] + (rmax-rmin)/NBINS; //rmaxs[i-1] = rmins[i-1] + (rmax-rmin)/NBINS; } - rmaxs[nshells-1] = rmax; + s->rmaxs[nshells-1] = rmax; + + return s; +} + + +static double shell_label(struct shells *s, int i) +{ + if ( s->config_intshells ) { + return (i+0.5) / s->nshells; + } else { + return s->rmins[i] + (s->rmaxs[i] - s->rmins[i])/2.0; + } +} + + +static int get_bin(struct shells *s, Reflection *refl, UnitCell *cell) +{ + if ( s->config_intshells ) { + + double intensity; + int bin, j; + + intensity = get_intensity(refl); + + bin = -1; + for ( j=0; jnshells; j++ ) { + if ( (intensity>s->rmins[j]) + && (intensity<=s->rmaxs[j]) ) + { + bin = j; + break; + } + } + + return bin; + + + } else { - for ( refl1 = first_refl(list1, &iter); - refl1 != NULL; - refl1 = next_refl(refl1, iter) ) - { - signed int h, k, l; double d; - int bin; - double i1, i2; - int j; - Reflection *refl2; + int bin, j; + signed int h, k, l; - get_indices(refl1, &h, &k, &l); + get_indices(refl, &h, &k, &l); d = 2.0 * resolution(cell, h, k, l); bin = -1; - for ( j=0; jrmins[j]) && (d<=rmaxs[j]) ) { + for ( j=0; jnshells; j++ ) { + if ( (d>s->rmins[j]) && (d<=s->rmaxs[j]) ) { bin = j; break; } } /* Allow for slight rounding errors */ - if ( (bin == -1) && (d <= rmins[0]) ) bin = 0; + if ( (bin == -1) && (d <= s->rmins[0]) ) bin = 0; assert(bin != -1); + return bin; + + } +} + + +static void do_fom(RefList *list1, RefList *list2, UnitCell *cell, + double rmin, double rmax, enum fom fom, + int config_unity, int nshells, const char *filename, + int config_intshells, double min_I, double max_I) +{ + int *cts; + int i; + Reflection *refl1; + RefListIterator *iter; + FILE *fh; + double scale; + struct fom_context *fctx; + struct shells *shells; + const char *t1, *t2; + int n_out; + + cts = malloc(nshells*sizeof(int)); + fctx = init_fom(fom, num_reflections(list1), nshells); + + if ( (fctx==NULL) || (cts==NULL) ) { + ERROR("Couldn't allocate memory for resolution shells.\n"); + return; + } + + for ( i=0; i max_I ) max_I = val1; + if ( val1 < min_I ) min_I = val1; + ncom++; } @@ -832,7 +994,7 @@ int main(int argc, char *argv[]) rmin/1e9, rmax/1e9, 1e10/rmin, 1e10/rmax); } do_fom(list1_acc, list2_acc, cell, rmin, rmax, fom, config_unity, - nshells, shell_file); + nshells, shell_file, config_intshells, min_I, max_I); free(shell_file); reflist_free(list1_acc); -- cgit v1.2.3