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 +++++++++++++++++++++++---------------------- libcrystfel/src/geometry.h | 3 +- libcrystfel/src/reflist.c | 37 +++++--------- libcrystfel/src/reflist.h | 7 ++- 4 files changed, 81 insertions(+), 88 deletions(-) (limited to 'libcrystfel') 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); diff --git a/libcrystfel/src/geometry.h b/libcrystfel/src/geometry.h index ccdea1c2..93aaf619 100644 --- a/libcrystfel/src/geometry.h +++ b/libcrystfel/src/geometry.h @@ -74,7 +74,8 @@ extern void update_partialities_2(Crystal *cryst, PartialityModel pmodel, extern void polarisation_correction(RefList *list, UnitCell *cell, struct image *image); -#define LORENTZ_SCALE (0.01e9) +extern double sphere_fraction(double rlow, double rhigh, double pr); +extern double gaussian_fraction(double rlow, double rhigh, double pr); #ifdef __cplusplus } diff --git a/libcrystfel/src/reflist.c b/libcrystfel/src/reflist.c index afc76caa..33c1f948 100644 --- a/libcrystfel/src/reflist.c +++ b/libcrystfel/src/reflist.c @@ -67,12 +67,10 @@ struct _refldata { signed int ls; /* Partiality and related geometrical stuff */ - double r1; /* First excitation error */ - double r2; /* Second excitation error */ + double rlow; /* Low excitation error */ + double rhigh; /* High excitation error */ double p; /* Partiality */ double L; /* Lorentz factor */ - int clamp1; /* Clamp status for r1 */ - int clamp2; /* Clamp status for r2 */ /* Location in image */ double fs; @@ -424,24 +422,20 @@ double get_intensity(const Reflection *refl) /** * get_partial: * @refl: A %Reflection - * @r1: Location at which to store the first excitation error - * @r2: Location at which to store the second excitation error + * @rlow: Location at which to store the "low" excitation error + * @rhigh: Location at which to store the "high" excitation error * @p: Location at which to store the partiality - * @clamp_low: Location at which to store the first clamp status - * @clamp_high: Location at which to store the second clamp status * * This function is used during post refinement (in conjunction with * set_partial()) to get access to the details of the partiality calculation. * **/ -void get_partial(const Reflection *refl, double *r1, double *r2, double *p, - int *clamp_low, int *clamp_high) +void get_partial(const Reflection *refl, double *rlow, double *rhigh, + double *p) { - *r1 = refl->data.r1; - *r2 = refl->data.r2; + *rlow = refl->data.rlow; + *rhigh = refl->data.rhigh; *p = get_partiality(refl); - *clamp_low = refl->data.clamp1; - *clamp_high = refl->data.clamp2; } @@ -592,24 +586,19 @@ void set_detector_pos(Reflection *refl, double exerr, double fs, double ss) /** * set_partial: * @refl: A %Reflection - * @r1: The first excitation error - * @r2: The second excitation error + * @rlow: The "low" excitation error + * @rhigh: The "high" excitation error * @p: The partiality - * @clamp_low: The first clamp status - * @clamp_high: The second clamp status * * This function is used during post refinement (in conjunction with * get_partial()) to get access to the details of the partiality calculation. * **/ -void set_partial(Reflection *refl, double r1, double r2, double p, - double clamp_low, double clamp_high) +void set_partial(Reflection *refl, double rlow, double rhigh, double p) { - refl->data.r1 = r1; - refl->data.r2 = r2; + refl->data.rlow = rlow; + refl->data.rhigh = rhigh; refl->data.p = p; - refl->data.clamp1 = clamp_low; - refl->data.clamp2 = clamp_high; } diff --git a/libcrystfel/src/reflist.h b/libcrystfel/src/reflist.h index 6892d7f2..85a87c54 100644 --- a/libcrystfel/src/reflist.h +++ b/libcrystfel/src/reflist.h @@ -94,8 +94,8 @@ extern void get_symmetric_indices(const Reflection *refl, signed int *hs, signed int *ks, signed int *ls); extern double get_intensity(const Reflection *refl); -extern void get_partial(const Reflection *refl, double *r1, double *r2, - double *p, int *clamp_low, int *clamp_high); +extern void get_partial(const Reflection *refl, double *rlow, double *rhigh, + double *p); extern int get_redundancy(const Reflection *refl); extern double get_temp1(const Reflection *refl); extern double get_temp2(const Reflection *refl); @@ -108,8 +108,7 @@ extern double get_mean_bg(const Reflection *refl); extern void copy_data(Reflection *to, const Reflection *from); extern void set_detector_pos(Reflection *refl, double exerr, double fs, double ss); -extern void set_partial(Reflection *refl, double r1, double r2, double p, - double clamp_low, double clamp_high); +extern void set_partial(Reflection *refl, double rlow, double rhigh, double p); extern void set_partiality(Reflection *refl, double p); extern void set_lorentz(Reflection *refl, double L); extern void set_intensity(Reflection *refl, double intensity); -- cgit v1.2.3