From 8bcc2f1fd10f5c1efdf89c00d04af1f0de948741 Mon Sep 17 00:00:00 2001 From: Thomas White Date: Thu, 26 Jun 2014 16:34:59 +0200 Subject: Add Thin Ewald Sphere model --- libcrystfel/src/geometry.c | 83 +++++++++++++++++++------- libcrystfel/src/geometry.h | 2 + src/partialator.c | 2 + src/post-refinement.c | 135 +++++++++++++++++++++++++------------------ tests/pr_l_gradient_check.c | 1 + tests/pr_p_gradient_check.c | 9 ++- tests/pr_pl_gradient_check.c | 3 +- 7 files changed, 155 insertions(+), 80 deletions(-) diff --git a/libcrystfel/src/geometry.c b/libcrystfel/src/geometry.c index 5cc46b2e..3586ed0b 100644 --- a/libcrystfel/src/geometry.c +++ b/libcrystfel/src/geometry.c @@ -103,7 +103,8 @@ static signed int locate_peak(double x, double y, double z, double k, } -static double partiality(PartialityModel pmodel, double rlow, double rhigh, +static double partiality(PartialityModel pmodel, + double rlow, double rmid, double rhigh, double r) { double qlow, qhigh; @@ -119,23 +120,22 @@ static double partiality(PartialityModel pmodel, double rlow, double rhigh, default: case PMODEL_UNITY: - plow = 1.0; - phigh = 0.0; - break; + return 1.0; case PMODEL_SPHERE: - plow = 3.0*pow(qlow,2.0) - 2.0*pow(qlow,3.0); - phigh = 3.0*pow(qhigh,2.0) - 2.0*pow(qhigh,3.0); - break; + plow = 3.0*qlow*qlow - 2.0*qlow*qlow*qlow; + phigh = 3.0*qhigh*qhigh - 2.0*qhigh*qhigh*qhigh; + return plow - phigh; case PMODEL_GAUSSIAN: plow = 0.5 * gsl_sf_erf(ng * rlow / (sqrt(2.0)*r)); phigh = 0.5 * gsl_sf_erf(ng * rhigh / (sqrt(2.0)*r)); - break; + return plow - phigh; - } + case PMODEL_THIN: + return 1.0 - (rmid*rmid)/(r*r); - return plow - phigh; + } } @@ -146,12 +146,12 @@ static Reflection *check_reflection(struct image *image, Crystal *cryst, { const int output = 0; double tl; - double rlow, rhigh; /* "Excitation error" */ + double rlow, rmid, rhigh; /* "Excitation error" */ double part; /* Partiality */ int clamp_low, clamp_high; - double klow, khigh; /* Wavenumber */ + double klow, kmid, khigh; /* Wavenumber */ Reflection *refl; - double cet, cez; + double cet, cez; /* Centre of Ewald sphere */ double pr; double L; double del; @@ -166,6 +166,7 @@ static Reflection *check_reflection(struct image *image, Crystal *cryst, * "high" gives the smallest Ewald sphere (wavelength long => k small) */ klow = 1.0/(image->lambda - image->lambda*image->bw/2.0); + kmid = 1.0/image->lambda; khigh = 1.0/(image->lambda + image->lambda*image->bw/2.0); /* If the point is looking "backscattery", reject it straight away */ @@ -177,14 +178,14 @@ static Reflection *check_reflection(struct image *image, Crystal *cryst, cez = -cos(del/2.0) * khigh; rhigh = khigh - distance(cet, cez, tl, zl); /* Loss of precision */ + cet = 0.0; + cez = -kmid; + rmid = kmid - distance(cet, cez, tl, zl); /* Loss of precision */ + cet = sin(del/2.0) * klow; cez = -cos(del/2.0) * klow; rlow = klow - distance(cet, cez, tl, zl); /* Loss of precision */ - if ( (signbit(rlow) == signbit(rhigh)) - && (fabs(rlow) > pr) - && (fabs(rhigh) > pr) ) return NULL; - if ( unlikely(rlow < rhigh) ) { ERROR("Reflection with rlow < rhigh!\n"); ERROR("%3i %3i %3i rlow = %e, rhigh = %e\n", @@ -193,10 +194,41 @@ static Reflection *check_reflection(struct image *image, Crystal *cryst, return NULL; } + /* Conditions for reflection to be excited at all */ + switch ( pmodel ) { + + default: + case PMODEL_UNITY: /* PMODEL_UNITY shouldn't end up here */ + case PMODEL_SPHERE: + case PMODEL_GAUSSIAN: + if ( (signbit(rlow) == signbit(rhigh)) + && (fabs(rlow) > pr) + && (fabs(rhigh) > pr) ) return NULL; + break; + + case PMODEL_THIN: + if ( fabs(rmid) > pr ) return NULL; + break; + + } + /* Lorentz factor is determined direction from the r values, before * clamping. The multiplication by 0.01e9 to make the * correction factor vaguely near 1. */ - L = LORENTZ_SCALE / (rlow - rhigh); + switch ( pmodel ) { + + default: + case PMODEL_SPHERE: + case PMODEL_GAUSSIAN: + L = LORENTZ_SCALE / (rlow - rhigh); + break; + + case PMODEL_UNITY: /* PMODEL_UNITY shouldn't end up here */ + case PMODEL_THIN: + L = 1.0; + break; + + } /* If the "lower" Ewald sphere is a long way away, use the * position at which the Ewald sphere would just touch the @@ -225,7 +257,7 @@ static Reflection *check_reflection(struct image *image, Crystal *cryst, } /* Calculate partiality */ - part = partiality(pmodel, rlow, rhigh, pr); + part = partiality(pmodel, rlow, rmid, rhigh, pr); /* Add peak to list */ refl = reflection_new(h, k, l); @@ -244,7 +276,18 @@ static Reflection *check_reflection(struct image *image, Crystal *cryst, set_detector_pos(refl, 0.0, xda, yda); } - set_partial(refl, rlow, rhigh, part, clamp_low, clamp_high); + if ( pmodel != PMODEL_THIN ) { + set_partial(refl, rlow, rhigh, part, clamp_low, clamp_high); + } else { + /* If we are using the TES (Thin Ewald Sphere) model, we abuse + * the fields as follows: + * rlow = the r value for the middle (only) Ewald sphere + * rhigh = 0.0 + * clamp_low = 0 + * clamp_high = +1 + */ + set_partial(refl, rmid, 0.0, part, 0, +1); + } set_lorentz(refl, L); set_symmetric_indices(refl, h, k, l); set_redundancy(refl, 1); diff --git a/libcrystfel/src/geometry.h b/libcrystfel/src/geometry.h index 1f465167..8041936a 100644 --- a/libcrystfel/src/geometry.h +++ b/libcrystfel/src/geometry.h @@ -50,6 +50,7 @@ extern "C" { * space. * @PMODEL_UNITY : Set all all partialities and Lorentz factors to 1. * @PMODEL_GAUSSIAN : Gaussian profiles in 3D + * @PMODEL_THIN : Thin Ewald sphere intersecting sphere * * A %PartialityModel describes a geometrical model which can be used to * calculate spot partialities and Lorentz correction factors. @@ -59,6 +60,7 @@ typedef enum { PMODEL_SPHERE, PMODEL_UNITY, PMODEL_GAUSSIAN, + PMODEL_THIN, } PartialityModel; diff --git a/src/partialator.c b/src/partialator.c index 2ef01dba..e5083611 100644 --- a/src/partialator.c +++ b/src/partialator.c @@ -348,6 +348,8 @@ int main(int argc, char *argv[]) pmodel = PMODEL_UNITY; } else if ( strcmp(pmodel_str, "gaussian") == 0 ) { pmodel = PMODEL_GAUSSIAN; + } else if ( strcmp(pmodel_str, "thin") == 0 ) { + pmodel = PMODEL_THIN; } else { ERROR("Unknown partiality model '%s'.\n", pmodel_str); return 1; diff --git a/src/post-refinement.c b/src/post-refinement.c index 5de8a246..99ddf6b1 100644 --- a/src/post-refinement.c +++ b/src/post-refinement.c @@ -55,15 +55,23 @@ #define MAX_CYCLES (10) -static double dpdq(double r, double profile_radius, PartialityModel pmodel) +static double dpdq(double r, double profile_radius) { double q; - double ng = 3.0; /* Calculate degree of penetration */ q = (r + profile_radius)/(2.0*profile_radius); - /* dp/dq */ + return 6.0*(q-q*q); +} + + +/* Returns dp/dr at "r" */ +static double partiality_gradient(double r, double profile_radius, + PartialityModel pmodel) +{ + double dqdr; /* dq/dr */ + switch ( pmodel ) { default: @@ -71,39 +79,46 @@ static double dpdq(double r, double profile_radius, PartialityModel pmodel) return 0.0; case PMODEL_SPHERE: - return 6.0*(q-pow(q, 2.0)); + dqdr = 1.0 / (2.0*profile_radius); + return dpdq(r, profile_radius) * dqdr; case PMODEL_GAUSSIAN: - /* The flat sphere model is close enough */ - return 6.0*(q-pow(q, 2.0)); + /* FIXME: Get a proper gradient */ + dqdr = 1.0 / (2.0*profile_radius); + return dpdq(r, profile_radius) * dqdr; + + case PMODEL_THIN: + return -(2.0*r)/(profile_radius*profile_radius); } } -/* Returns dp/dr at "r" */ -static double partiality_gradient(double r, double profile_radius, - PartialityModel pmodel) +/* Returns dp/drad at "r" */ +static double partiality_rgradient(double r, double profile_radius, + PartialityModel pmodel) { - double dqdr; + double dqdrad; /* dq/drad */ - /* dq/dr */ - dqdr = 1.0 / (2.0*profile_radius); + switch ( pmodel ) { - return dpdq(r, profile_radius, pmodel) * dqdr; -} + default: + case PMODEL_UNITY: + return 0.0; + case PMODEL_SPHERE: + dqdrad = -0.5 * r / (profile_radius * profile_radius); + return dpdq(r, profile_radius) * dqdrad; -/* Returns dp/drad at "r" */ -static double partiality_rgradient(double r, double profile_radius, - PartialityModel pmodel) -{ - double dqdrad; + case PMODEL_GAUSSIAN: + /* FIXME: Get a proper gradient */ + dqdrad = -0.5 * r / (profile_radius * profile_radius); + return dpdq(r, profile_radius) * dqdrad; - /* dq/drad */ - dqdrad = -0.5 * r * pow(profile_radius, -2.0); + case PMODEL_THIN: + return 2.0*r*r*pow(profile_radius, -3.0); - return dpdq(r, profile_radius, pmodel) * dqdrad; + } } @@ -173,19 +188,48 @@ double p_gradient(Crystal *cr, int k, Reflection *refl, PartialityModel pmodel) ghigh = 0.0; } + if ( k == REF_R ) { + switch ( pmodel ) { + + default: + case PMODEL_UNITY: + return 0.0; + + case PMODEL_GAUSSIAN: + gr = partiality_rgradient(rlow, r, pmodel); + gr -= partiality_rgradient(rhigh, r, pmodel); + return gr; + + case PMODEL_SPHERE: + gr = partiality_rgradient(rlow, r, pmodel); + gr -= partiality_rgradient(rhigh, r, pmodel); + return gr; + + case PMODEL_THIN: + return 2.0*rlow*rlow/(r*r*r); + } + } + + if ( k == REF_DIV ) { + switch ( pmodel ) { + + default: + case PMODEL_UNITY: + return 0.0; + + case PMODEL_GAUSSIAN: + case PMODEL_SPHERE: + return (ds*glow + ds*ghigh) / 2.0; + + case PMODEL_THIN: + return 0.0; + } + } + /* For many gradients, just multiply the above number by the gradient * of excitation error wrt whatever. */ switch ( k ) { - case REF_DIV : - /* Small angle approximation */ - return (ds*glow + ds*ghigh) / 2.0; - - case REF_R : - gr = partiality_rgradient(rlow, r, pmodel); - gr -= partiality_rgradient(rhigh, r, pmodel); - return gr; - /* Cell parameters and orientation */ case REF_ASX : return hs * sin(phi) * cos(azi) * (ghigh-glow); @@ -227,37 +271,16 @@ double l_gradient(Crystal *cr, int k, Reflection *refl, PartialityModel pmodel) { double ds; signed int hs, ks, ls; + double L; - switch ( k ) { - - /* Cell parameters do not affect Lorentz factor */ - case REF_ASX : - case REF_BSX : - case REF_CSX : - case REF_ASY : - case REF_BSY : - case REF_CSY : - case REF_ASZ : - case REF_BSZ : - case REF_CSZ : - return 0.0; - - /* Nor does change of radius */ - case REF_R : - return 0.0; - - default: - break; - - } - - assert(k == REF_DIV); + if ( k != REF_DIV ) return 0.0; get_symmetric_indices(refl, &hs, &ks, &ls); ds = 2.0 * resolution(crystal_get_cell(cr), hs, ks, ls); - return -ds*pow(get_lorentz(refl), 2.0) / LORENTZ_SCALE; + L = get_lorentz(refl); + return -ds*L*L / LORENTZ_SCALE; } diff --git a/tests/pr_l_gradient_check.c b/tests/pr_l_gradient_check.c index b707629f..7342d5f8 100644 --- a/tests/pr_l_gradient_check.c +++ b/tests/pr_l_gradient_check.c @@ -350,6 +350,7 @@ int main(int argc, char *argv[]) pmodel = PMODEL_GAUSSIAN; STATUS("Testing Gaussian model:\n"); } + /* No point testing TES model, because it has no Lorentz factor */ orientation = random_quaternion(rng); rot = cell_rotate(cell, orientation); diff --git a/tests/pr_p_gradient_check.c b/tests/pr_p_gradient_check.c index 8f7a66e4..14a2b579 100644 --- a/tests/pr_p_gradient_check.c +++ b/tests/pr_p_gradient_check.c @@ -174,6 +174,7 @@ static Crystal *new_shifted_crystal(Crystal *cr, int refine, double incr_val) return cr_new; } + static void calc_either_side(Crystal *cr, double incr_val, int *valid, long double *vals[3], int refine, PartialityModel pmodel) @@ -225,7 +226,6 @@ static void calc_either_side(Crystal *cr, double incr_val, } - static double test_gradients(Crystal *cr, double incr_val, int refine, const char *str, const char *file, PartialityModel pmodel, int quiet, int plot) @@ -450,7 +450,7 @@ int main(int argc, char *argv[]) rng = gsl_rng_alloc(gsl_rng_mt19937); - for ( i=0; i<2; i++ ) { + for ( i=0; i<3; i++ ) { UnitCell *rot; double val; @@ -459,9 +459,12 @@ int main(int argc, char *argv[]) if ( i == 0 ) { pmodel = PMODEL_SPHERE; STATUS("Testing flat sphere model:\n"); - } else { + } else if ( i == 1 ) { pmodel = PMODEL_GAUSSIAN; STATUS("Testing Gaussian model:\n"); + } else { + pmodel = PMODEL_THIN; + STATUS("Testing thin Ewald sphere model:\n"); } orientation = random_quaternion(rng); diff --git a/tests/pr_pl_gradient_check.c b/tests/pr_pl_gradient_check.c index f83af335..69087a4b 100644 --- a/tests/pr_pl_gradient_check.c +++ b/tests/pr_pl_gradient_check.c @@ -174,6 +174,7 @@ static Crystal *new_shifted_crystal(Crystal *cr, int refine, double incr_val) return cr_new; } + static void calc_either_side(Crystal *cr, double incr_val, int *valid, long double *vals[3], int refine, PartialityModel pmodel) @@ -225,7 +226,6 @@ static void calc_either_side(Crystal *cr, double incr_val, } - static double test_gradients(Crystal *cr, double incr_val, int refine, const char *str, const char *file, PartialityModel pmodel, int quiet, int plot) @@ -465,6 +465,7 @@ int main(int argc, char *argv[]) pmodel = PMODEL_GAUSSIAN; STATUS("Testing Gaussian model:\n"); } + /* No point testing TES model */ orientation = random_quaternion(rng); rot = cell_rotate(cell, orientation); -- cgit v1.2.3