From 761da16096472d157b2dd9c871e81ab17cdef12d Mon Sep 17 00:00:00 2001 From: Thomas White Date: Wed, 6 Feb 2013 18:59:40 +0100 Subject: Stuff for partialator --- src/post-refinement.c | 64 +++++++++++++++++++++++++++++---------------------- 1 file changed, 36 insertions(+), 28 deletions(-) (limited to 'src/post-refinement.c') diff --git a/src/post-refinement.c b/src/post-refinement.c index 7410d931..31558ec7 100644 --- a/src/post-refinement.c +++ b/src/post-refinement.c @@ -3,11 +3,11 @@ * * Post refinement * - * Copyright © 2012 Deutsches Elektronen-Synchrotron DESY, - * a research centre of the Helmholtz Association. + * Copyright © 2012-2013 Deutsches Elektronen-Synchrotron DESY, + * a research centre of the Helmholtz Association. * * Authors: - * 2010-2012 Thomas White + * 2010-2013 Thomas White * * This file is part of CrystFEL. * @@ -88,7 +88,7 @@ static double partiality_rgradient(double r, double profile_radius) /* Return the gradient of parameter 'k' given the current status of 'image'. */ -double gradient(struct image *image, int k, Reflection *refl, double r) +double gradient(Crystal *cr, int k, Reflection *refl, double r) { double ds, azix, aziy; double ttlow, tthigh, tt; @@ -103,17 +103,18 @@ double gradient(struct image *image, int k, Reflection *refl, double r) int clamp_low, clamp_high; double klow, khigh; double gr; + struct image *image = crystal_get_image(cr); get_symmetric_indices(refl, &hs, &ks, &ls); - cell_get_reciprocal(image->indexed_cell, &asx, &asy, &asz, - &bsx, &bsy, &bsz, - &csx, &csy, &csz); + cell_get_reciprocal(crystal_get_cell(cr), &asx, &asy, &asz, + &bsx, &bsy, &bsz, + &csx, &csy, &csz); xl = hs*asx + ks*bsx + ls*csx; yl = hs*asy + ks*bsy + ls*csy; zl = hs*asz + ks*bsz + ls*csz; - ds = 2.0 * resolution(image->indexed_cell, hs, ks, ls); + ds = 2.0 * resolution(crystal_get_cell(cr), hs, ks, ls); get_partial(refl, &r1, &r2, &p, &clamp_low, &clamp_high); klow = 1.0/(image->lambda - image->lambda*image->bw/2.0); @@ -237,8 +238,11 @@ static void apply_cell_shift(UnitCell *cell, int k, double shift) /* Apply the given shift to the 'k'th parameter of 'image'. */ -static void apply_shift(struct image *image, int k, double shift) +static void apply_shift(Crystal *cr, int k, double shift) { + double t; + struct image *image = crystal_get_image(cr); + switch ( k ) { case REF_DIV : @@ -251,7 +255,9 @@ static void apply_shift(struct image *image, int k, double shift) break; case REF_R : - image->profile_radius += shift; + t = crystal_get_profile_radius(cr); + t += shift; + crystal_set_profile_radius(cr, t); break; case REF_ASX : @@ -263,7 +269,7 @@ static void apply_shift(struct image *image, int k, double shift) case REF_CSX : case REF_CSY : case REF_CSZ : - apply_cell_shift(image->indexed_cell, k, shift); + apply_cell_shift(crystal_get_cell(cr), k, shift); break; default : @@ -357,7 +363,7 @@ static gsl_vector *solve_svd(gsl_vector *v, gsl_matrix *M) /* Perform one cycle of post refinement on 'image' against 'full' */ -static double pr_iterate(struct image *image, const RefList *full) +static double pr_iterate(Crystal *cr, const RefList *full) { gsl_matrix *M; gsl_vector *v; @@ -369,7 +375,7 @@ static double pr_iterate(struct image *image, const RefList *full) double max_shift; int nref = 0; - reflections = image->reflections; + reflections = crystal_get_reflections(cr); M = gsl_matrix_calloc(NUM_PARAMS, NUM_PARAMS); v = gsl_vector_calloc(NUM_PARAMS); @@ -387,6 +393,7 @@ static double pr_iterate(struct image *image, const RefList *full) double p; Reflection *match; double gradients[NUM_PARAMS]; + const double osf = crystal_get_osf(cr); if ( !get_refinable(refl) ) continue; @@ -398,7 +405,7 @@ static double pr_iterate(struct image *image, const RefList *full) " is it marked as refinable?\n", ha, ka, la); continue; } - I_full = image->osf * get_intensity(match); + I_full = osf * get_intensity(match); /* Actual measurement of this reflection from this pattern? */ I_partial = get_intensity(refl); @@ -412,7 +419,8 @@ static double pr_iterate(struct image *image, const RefList *full) /* Calculate all gradients for this reflection */ for ( k=0; kprofile_radius); + gr = gradient(cr, k, refl, + crystal_get_profile_radius(cr)); gradients[k] = gr; } @@ -429,7 +437,7 @@ static double pr_iterate(struct image *image, const RefList *full) if ( g > k ) continue; M_c = gradients[g] * gradients[k]; - M_c *= w * pow(image->osf * I_full, 2.0); + M_c *= w * pow(osf * I_full, 2.0); M_curr = gsl_matrix_get(M, k, g); gsl_matrix_set(M, k, g, M_curr + M_c); @@ -450,7 +458,7 @@ static double pr_iterate(struct image *image, const RefList *full) //STATUS("%i reflections went into the equations.\n", nref); if ( nref == 0 ) { - image->pr_dud = 1; + crystal_set_user_flag(cr, 1); gsl_matrix_free(M); gsl_vector_free(v); return 0.0; @@ -462,7 +470,7 @@ static double pr_iterate(struct image *image, const RefList *full) for ( param=0; param max_shift ) max_shift = fabs(shift); } @@ -480,7 +488,7 @@ static double pr_iterate(struct image *image, const RefList *full) } -static double guide_dev(struct image *image, const RefList *full) +static double guide_dev(Crystal *cr, const RefList *full) { double dev = 0.0; @@ -488,7 +496,7 @@ static double guide_dev(struct image *image, const RefList *full) Reflection *refl; RefListIterator *iter; - for ( refl = first_refl(image->reflections, &iter); + for ( refl = first_refl(crystal_get_reflections(cr), &iter); refl != NULL; refl = next_refl(refl, iter) ) { @@ -508,7 +516,7 @@ static double guide_dev(struct image *image, const RefList *full) * scale_intensities() might not yet have been called, so the * full version may not have been calculated yet. */ - G = image->osf; + G = crystal_get_osf(cr); p = get_partiality(refl); I_partial = get_intensity(refl); I_full = get_intensity(full_version); @@ -524,21 +532,21 @@ static double guide_dev(struct image *image, const RefList *full) } -void pr_refine(struct image *image, const RefList *full) +void pr_refine(Crystal *cr, const RefList *full) { double max_shift, dev; int i; const int verbose = 0; if ( verbose ) { - dev = guide_dev(image, full); + dev = guide_dev(cr, full); STATUS("\n"); /* Deal with progress bar */ STATUS("Before iteration: dev = %10.5e\n", dev); } i = 0; - image->pr_dud = 0; + crystal_set_user_flag(cr, 0); do { double asx, asy, asz; @@ -546,15 +554,15 @@ void pr_refine(struct image *image, const RefList *full) double csx, csy, csz; double dev; - cell_get_reciprocal(image->indexed_cell, &asx, &asy, &asz, + cell_get_reciprocal(crystal_get_cell(cr), &asx, &asy, &asz, &bsx, &bsy, &bsz, &csx, &csy, &csz); - max_shift = pr_iterate(image, full); + max_shift = pr_iterate(cr, full); - update_partialities(image); + update_partialities(cr); if ( verbose ) { - dev = guide_dev(image, full); + dev = guide_dev(cr, full); STATUS("PR Iteration %2i: max shift = %10.2f" " dev = %10.5e\n", i+1, max_shift, dev); } -- cgit v1.2.3