From 20d05388833bd9f3b6819d82f983eea424663407 Mon Sep 17 00:00:00 2001 From: Thomas White Date: Mon, 22 Nov 2010 11:19:23 +0100 Subject: Fix post refinement maths and take clamping status into account --- src/facetron.c | 1 + src/geometry.c | 52 +++++++++++++++++++++++++++++--------------------- src/image.h | 3 +++ src/post-refinement.c | 53 ++++++++++++++++++++++++++++++++++++++++++--------- src/post-refinement.h | 5 ----- 5 files changed, 78 insertions(+), 36 deletions(-) (limited to 'src') diff --git a/src/facetron.c b/src/facetron.c index f561b554..2bac7511 100644 --- a/src/facetron.c +++ b/src/facetron.c @@ -467,6 +467,7 @@ int main(int argc, char *argv[]) images[i].det = det; images[i].beam = beam; images[i].osf = 1.0; + images[i].profile_radius = 0.005e9; /* Muppet proofing */ images[i].data = NULL; diff --git a/src/geometry.c b/src/geometry.c index 0ef1f1c3..664817e0 100644 --- a/src/geometry.c +++ b/src/geometry.c @@ -103,26 +103,17 @@ static double excitation_error(double xl, double yl, double zl, static double partiality(double r1, double r2, double r) { double q1, q2; - double p, p1, p2; + double p1, p2; /* Calculate degrees of penetration */ q1 = (r1 + r)/(2.0*r); q2 = (r2 + r)/(2.0*r); - /* Clamp */ - if ( q1 > 1.0 ) q1 = 1.0; - if ( q1 < 0.0 ) q1 = 0.0; - if ( q2 > 1.0 ) q2 = 1.0; - if ( q2 < 0.0 ) q2 = 0.0; - /* Convert to partiality */ p1 = 3.0*pow(q1,2.0) - 2.0*pow(q1,3.0); p2 = 3.0*pow(q2,2.0) - 2.0*pow(q2,3.0); - /* Input values may have been backwards */ - p = fabs(p2 - p1); - - return p; + return p2 - p1; } @@ -143,8 +134,6 @@ struct cpeak *find_intersections(struct image *image, UnitCell *cell, double klow, kcen, khigh; /* Wavenumber */ /* Bounding sphere for the shape transform approximation */ const double profile_cutoff = 0.02e9; /* 0.02 nm^-1 */ - /* Actual radius of the profile */ - const double profile_radius = 0.005e9; cpeaks = malloc(sizeof(struct cpeak)*MAX_CPEAKS); if ( cpeaks == NULL ) { @@ -178,6 +167,8 @@ struct cpeak *find_intersections(struct image *image, UnitCell *cell, double xda, yda; /* Position on detector */ int close, inside; double part; /* Partiality */ + int clamp_low = 0; + int clamp_high = 0; /* Ignore central beam */ if ( (h==0) && (k==0) && (l==0) ) continue; @@ -216,21 +207,36 @@ struct cpeak *find_intersections(struct image *image, UnitCell *cell, /* If the "lower" Ewald sphere is a long way away, use the * position at which the Ewald sphere would just touch the * reflection. */ - if ( rlow > profile_cutoff ) rlow = profile_cutoff; - if ( rlow < -profile_cutoff ) rlow = -profile_cutoff; + if ( rlow < -profile_cutoff ) { + rlow = -profile_cutoff; + clamp_low = -1; + } + if ( rlow > +profile_cutoff ) { + rlow = +profile_cutoff; + clamp_low = +1; + } + /* Likewise the "higher" Ewald sphere */ + if ( rhigh < -profile_cutoff ) { + rhigh = -profile_cutoff; + clamp_high = -1; + } + if ( rhigh > +profile_cutoff ) { + rhigh = +profile_cutoff; + clamp_high = +1; + } + /* The six possible combinations of clamp_{low,high} (including + * zero) correspond to the six situations in Table 3 of Rossmann + * et al. (1979). */ - /* As above, other side. */ - if ( rhigh > profile_cutoff ) rhigh = profile_cutoff; - if ( rhigh < -profile_cutoff ) rhigh = -profile_cutoff; + /* Calculate partiality and reject if too small */ + part = partiality(rlow, rhigh, image->profile_radius); + if ( part < 0.1 ) continue; /* Locate peak on detector. */ p = locate_peak(xl, yl, zl, kcen, image->det, &xda, &yda); if ( p == -1 ) continue; - part = partiality(rlow, rhigh, profile_radius); - - if ( part < 0.1 ) continue; - + /* Add peak to list */ cpeaks[np].h = h; cpeaks[np].k = k; cpeaks[np].l = l; @@ -239,6 +245,8 @@ struct cpeak *find_intersections(struct image *image, UnitCell *cell, cpeaks[np].r1 = rlow; cpeaks[np].r2 = rhigh; cpeaks[np].p = part; + cpeaks[np].clamp1 = clamp_low; + cpeaks[np].clamp2 = clamp_high; np++; if ( output ) { diff --git a/src/image.h b/src/image.h index ae2f3179..0d2f6c2e 100644 --- a/src/image.h +++ b/src/image.h @@ -66,6 +66,8 @@ struct cpeak double r1; /* First excitation error */ double r2; /* Second excitation error */ double p; /* Partiality */ + int clamp1; /* Clamp status for r1 */ + int clamp2; /* Clamp status for r2 */ /* Location in image */ int x; @@ -102,6 +104,7 @@ struct image { int f0_available; /* 0 if f0 wasn't available * from the input. */ double osf; /* Overall scaling factor */ + double profile_radius; /* Radius of reflection */ int width; int height; diff --git a/src/post-refinement.c b/src/post-refinement.c index 5ca0dc98..e2fa6d50 100644 --- a/src/post-refinement.c +++ b/src/post-refinement.c @@ -28,12 +28,33 @@ #include "geometry.h" +/* Returns dp/dr at "r" */ +static double partiality_gradient(double r, double profile_radius) +{ + double q, dpdq, dqdr; + + /* Calculate degree of penetration */ + q = (r + profile_radius)/(2.0*profile_radius); + + /* dp/dq */ + dpdq = 6.0*(q-pow(q, 2.0)); + + /* dq/dr */ + dqdr = 1.0 / (2.0*profile_radius); + + return dpdq * dqdr; +} + + /* Return the gradient of parameter 'k' given the current status of 'image'. */ double gradient(struct image *image, int k, - struct cpeak spot, double I_partial) + struct cpeak spot, double I_partial, double r) { double ds; double nom, den; + double r1g = 0.0; + double r2g = 0.0; + double g = 0.0; ds = 2.0 * resolution(image->indexed_cell, spot.h, spot.k, spot.l); @@ -43,9 +64,20 @@ double gradient(struct image *image, int k, return I_partial; case REF_DIV : - nom = sqrt(2.0) * ds * sin(image->div); - den = sqrt(1.0 - cos(image->div)); - return nom/den; + /* Calculate the gradient of r1 and r2 wrt divergence */ + if ( spot.clamp1 == 0 ) { + nom = sqrt(2.0) * ds * sin(image->div/2.0); + den = sqrt(1.0 - cos(image->div/2.0)); + r1g = nom/den; + g += r1g * partiality_gradient(spot.r1, r); + } + if ( spot.clamp2 == 0 ) { + nom = sqrt(2.0) * ds * sin(image->div/2.0); + den = sqrt(1.0 - cos(image->div/2.0)); + r2g = nom/den; + g += r2g * partiality_gradient(spot.r2, r); + } + return g; } @@ -171,7 +203,7 @@ double pr_iterate(struct image *image, double *i_full, const char *sym, for ( k=0; kprofile_radius) + * gradient(image, k, spots[h], I_partial, + image->profile_radius); M_c *= pow(I_full, 2.0); gsl_matrix_set(M, g, k, M_curr + M_c); } - v_c = delta_I * I_full * gradient(image, k, spots[h], - I_partial); + gr = gradient(image, k, spots[h], I_partial, + image->profile_radius); + v_c = delta_I * I_full * gr; gsl_vector_set(v, k, v_c); } diff --git a/src/post-refinement.h b/src/post-refinement.h index 1ef7a1fd..8dea914d 100644 --- a/src/post-refinement.h +++ b/src/post-refinement.h @@ -30,11 +30,6 @@ enum { NUM_PARAMS }; - -/* Return the gradient of parameter 'k' given the current status of 'image'. */ -double gradient(struct image *image, int k, struct cpeak spot, - double I_partial); - /* Apply the given shift to the 'k'th parameter of 'image'. */ void apply_shift(struct image *image, int k, double shift); -- cgit v1.2.3