From 015ffcd77acc4b7f3178753b21e6591fd8212e4c Mon Sep 17 00:00:00 2001 From: Thomas White Date: Tue, 30 Sep 2014 16:58:29 +0200 Subject: Cell vector gradients for SCSphere, plus general rationalisation --- libcrystfel/src/geometry.c | 122 +++++++++++++++++++++++---------------------- 1 file changed, 63 insertions(+), 59 deletions(-) (limited to 'libcrystfel/src/geometry.c') diff --git a/libcrystfel/src/geometry.c b/libcrystfel/src/geometry.c index 250f889c..e08d9120 100644 --- a/libcrystfel/src/geometry.c +++ b/libcrystfel/src/geometry.c @@ -102,17 +102,64 @@ static signed int locate_peak(double x, double y, double z, double k, } -static double partiality(PartialityModel pmodel, - double rlow, double rmid, double rhigh, - double r, double D) +double sphere_fraction(double rlow, double rhigh, double pr) { double qlow, qhigh; double plow, phigh; - const double ng = 2.6; + + /* If the "lower" Ewald sphere is a long way away, use the + * position at which the Ewald sphere would just touch the + * reflection. + * + * The six possible combinations of clamp_{low,high} (including + * zero) correspond to the six situations in Table 3 of Rossmann + * et al. (1979). + */ + if ( rlow < -pr ) rlow = -pr; + if ( rlow > +pr ) rlow = +pr; + if ( rhigh < -pr ) rhigh = -pr; + if ( rhigh > +pr ) rhigh = +pr; /* Calculate degrees of penetration */ - qlow = (rlow + r)/(2.0*r); - qhigh = (rhigh + r)/(2.0*r); + qlow = (rlow + pr)/(2.0*pr); + qhigh = (rhigh + pr)/(2.0*pr); + + plow = 3.0*qlow*qlow - 2.0*qlow*qlow*qlow; + phigh = 3.0*qhigh*qhigh - 2.0*qhigh*qhigh*qhigh; + + return plow - phigh; +} + + +double gaussian_fraction(double rlow, double rhigh, double pr) +{ + double plow, phigh; + const double ng = 2.6; + + /* If the "lower" Ewald sphere is a long way away, use the + * position at which the Ewald sphere would just touch the + * reflection. + * + * The six possible combinations of clamp_{low,high} (including + * zero) correspond to the six situations in Table 3 of Rossmann + * et al. (1979). + */ + if ( rlow < -pr ) rlow = -pr; + if ( rlow > +pr ) rlow = +pr; + if ( rhigh < -pr ) rhigh = -pr; + if ( rhigh > +pr ) rhigh = +pr; + + plow = 0.5 * gsl_sf_erf(ng * rlow / (sqrt(2.0)*pr)); + phigh = 0.5 * gsl_sf_erf(ng * rhigh / (sqrt(2.0)*pr)); + + return plow - phigh; +} + + +static double partiality(PartialityModel pmodel, + double rlow, double rhigh, double pr) +{ + double D = rlow - rhigh; /* Convert to partiality */ switch ( pmodel ) { @@ -121,15 +168,11 @@ static double partiality(PartialityModel pmodel, case PMODEL_UNITY: return 1.0; - case PMODEL_SCGAUSSIAN: - plow = 0.5 * gsl_sf_erf(ng * rlow / (sqrt(2.0)*r)); - phigh = 0.5 * gsl_sf_erf(ng * rhigh / (sqrt(2.0)*r)); - return 4.0*(plow-phigh)*r / (3.0*D); - case PMODEL_SCSPHERE: - plow = 3.0*qlow*qlow - 2.0*qlow*qlow*qlow; - phigh = 3.0*qhigh*qhigh - 2.0*qhigh*qhigh*qhigh; - return 4.0*(plow-phigh)*r / (3.0*D); + return 4.0*sphere_fraction(rlow, rhigh, pr)*pr / (3.0*D); + + case PMODEL_SCGAUSSIAN: + return 4.0*gaussian_fraction(rlow, rhigh, pr)*pr / (3.0*D); } } @@ -142,16 +185,13 @@ static Reflection *check_reflection(struct image *image, Crystal *cryst, { const int output = 0; double tl; - double rlow, rmid, rhigh; /* "Excitation error" */ + double rlow, rhigh; /* "Excitation error" */ double part; /* Partiality */ - int clamp_low, clamp_high; - double klow, kmid, khigh; /* Wavenumber */ + double klow, khigh; /* Wavenumber */ Reflection *refl; double cet, cez; /* Centre of Ewald sphere */ double pr; - double D; double del; - double rlowuc, rhighuc; /* Values before clamping */ /* Don't predict 000 */ if ( abs(h)+abs(k)+abs(l) == 0 ) return NULL; @@ -163,7 +203,6 @@ 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 */ @@ -175,10 +214,6 @@ 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 */ @@ -196,38 +231,8 @@ static Reflection *check_reflection(struct image *image, Crystal *cryst, && (fabs(rlow) > pr) && (fabs(rhigh) > pr) ) return NULL; - D = rlow - rhigh; - - rlowuc = rlow; rhighuc = rhigh; - - /* If the "lower" Ewald sphere is a long way away, use the - * position at which the Ewald sphere would just touch the - * reflection. - * - * The six possible combinations of clamp_{low,high} (including - * zero) correspond to the six situations in Table 3 of Rossmann - * et al. (1979). - */ - clamp_low = 0; clamp_high = 0; - if ( rlow < -pr ) { - rlow = -pr; - clamp_low = -1; - } - if ( rlow > +pr ) { - rlow = +pr; - clamp_low = +1; - } - if ( rhigh < -pr ) { - rhigh = -pr; - clamp_high = -1; - } - if ( rhigh > +pr ) { - rhigh = +pr; - clamp_high = +1; - } - /* Calculate partiality */ - part = partiality(pmodel, rlow, rmid, rhigh, pr, D); + part = partiality(pmodel, rlow, rhigh, pr); /* Add peak to list */ refl = reflection_new(h, k, l); @@ -246,7 +251,7 @@ static Reflection *check_reflection(struct image *image, Crystal *cryst, set_detector_pos(refl, 0.0, xda, yda); } - set_partial(refl, rlowuc, rhighuc, part, clamp_low, clamp_high); + set_partial(refl, rlow, rhigh, part); set_lorentz(refl, 1.0); set_symmetric_indices(refl, h, k, l); set_redundancy(refl, 1); @@ -442,7 +447,6 @@ void update_partialities_2(Crystal *cryst, PartialityModel pmodel, double r1, r2, L, p, x, y; double xl, yl, zl; signed int h, k, l; - int clamp1, clamp2; double old_p; get_symmetric_indices(refl, &h, &k, &l); @@ -472,8 +476,8 @@ void update_partialities_2(Crystal *cryst, PartialityModel pmodel, } /* Transfer partiality stuff */ - get_partial(vals, &r1, &r2, &p, &clamp1, &clamp2); - set_partial(refl, r1, r2, p, clamp1, clamp2); + get_partial(vals, &r1, &r2, &p); + set_partial(refl, r1, r2, p); L = get_lorentz(vals); set_lorentz(refl, L); -- cgit v1.2.3