diff options
Diffstat (limited to 'libcrystfel')
| -rw-r--r-- | libcrystfel/src/geometry.c | 48 | ||||
| -rw-r--r-- | libcrystfel/src/geometry.h | 18 | ||||
| -rw-r--r-- | libcrystfel/src/reflist.c | 51 | ||||
| -rw-r--r-- | libcrystfel/src/reflist.h | 9 | ||||
| -rw-r--r-- | libcrystfel/src/utils.h | 20 |
5 files changed, 133 insertions, 13 deletions
diff --git a/libcrystfel/src/geometry.c b/libcrystfel/src/geometry.c index 621c9484..cdc62ab5 100644 --- a/libcrystfel/src/geometry.c +++ b/libcrystfel/src/geometry.c @@ -134,6 +134,10 @@ static Reflection *check_reflection(struct image *image, Crystal *cryst, Reflection *refl; double cet, cez; double pr; + double L; + + /* Don't predict 000 */ + if ( abs(h)+abs(k)+abs(l) == 0 ) return NULL; pr = crystal_get_profile_radius(cryst); @@ -160,6 +164,19 @@ static Reflection *check_reflection(struct image *image, Crystal *cryst, && (fabs(rlow) > pr) && (fabs(rhigh) > pr) ) return NULL; + if ( rlow < rhigh ) { + ERROR("Reflection with rlow < rhigh!\n"); + ERROR("%3i %3i %3i rlow = %e, rhigh = %e\n", + h, k, l, rlow, rhigh); + ERROR("div = %e\n", image->div); + return NULL; + } + + /* Lorentz factor is determined direction from the r values, before + * clamping. The multiplication by the profile radius is to make the + * correction factor vaguely near 1. */ + L = pr / (rlow - rhigh); + /* If the "lower" Ewald sphere is a long way away, use the * position at which the Ewald sphere would just touch the * reflection. @@ -185,7 +202,6 @@ static Reflection *check_reflection(struct image *image, Crystal *cryst, rhigh = +pr; clamp_high = +1; } - assert(clamp_low >= clamp_high); /* Calculate partiality */ part = partiality(rlow, rhigh, pr); @@ -198,6 +214,7 @@ static Reflection *check_reflection(struct image *image, Crystal *cryst, refl = reflection_new(h, k, l); set_detector_pos(refl, 0.0, xda, yda); set_partial(refl, rlow, rhigh, part, clamp_low, clamp_high); + set_lorentz(refl, L); set_symmetric_indices(refl, h, k, l); set_redundancy(refl, 1); @@ -345,8 +362,23 @@ RefList *select_intersections(struct image *image, Crystal *cryst) } +static void set_unity_partialities(Crystal *cryst) +{ + Reflection *refl; + RefListIterator *iter; + + for ( refl = first_refl(crystal_get_reflections(cryst), &iter); + refl != NULL; + refl = next_refl(refl, iter) ) + { + set_partiality(refl, 1.0); + set_lorentz(refl, 1.0); + } +} + + /* Calculate partialities and apply them to the image's reflections */ -void update_partialities(Crystal *cryst) +void update_partialities(Crystal *cryst, PartialityModel pmodel) { Reflection *refl; RefListIterator *iter; @@ -356,6 +388,14 @@ void update_partialities(Crystal *cryst) double csx, csy, csz; struct image *image = crystal_get_image(cryst); + if ( pmodel == PMODEL_UNITY ) { + /* It isn't strictly necessary to set the partialities to 1, + * because the scaling stuff will just a correction factor of + * 1 anyway. This is just to help things make sense. */ + set_unity_partialities(cryst); + return; + } + cell_get_reciprocal(crystal_get_cell(cryst), &asx, &asy, &asz, &bsx, &bsy, &bsz, &csx, &csy, &csz); @@ -367,7 +407,7 @@ void update_partialities(Crystal *cryst) refl = next_refl(refl, iter) ) { Reflection *vals; - double r1, r2, p, x, y; + double r1, r2, L, p, x, y; double xl, yl, zl; signed int h, k, l; int clamp1, clamp2; @@ -390,6 +430,8 @@ void update_partialities(Crystal *cryst) /* Transfer partiality stuff */ get_partial(vals, &r1, &r2, &p, &clamp1, &clamp2); set_partial(refl, r1, r2, p, clamp1, clamp2); + L = get_lorentz(vals); + set_lorentz(refl, L); /* Transfer detector location */ get_detector_pos(vals, &x, &y); diff --git a/libcrystfel/src/geometry.h b/libcrystfel/src/geometry.h index aecdc28a..2ac5982b 100644 --- a/libcrystfel/src/geometry.h +++ b/libcrystfel/src/geometry.h @@ -43,10 +43,26 @@ extern "C" { #endif +/** + * PartialityModel: + * @PMODEL_SPHERE : Intersection of sphere with excited volume of reciprocal + * space. + * @PMODEL_UNITY : Set all all partialities and Lorentz factors to 1. + * + * A %PartialityModel describes a geometrical model which can be used to + * calculate spot partialities and Lorentz correction factors. + **/ +typedef enum { + + PMODEL_SPHERE, + PMODEL_UNITY, + +} PartialityModel; + extern RefList *find_intersections(struct image *image, Crystal *cryst); extern RefList *select_intersections(struct image *image, Crystal *cryst); -extern void update_partialities(Crystal *cryst); +extern void update_partialities(Crystal *cryst, PartialityModel pmodel); #ifdef __cplusplus } diff --git a/libcrystfel/src/reflist.c b/libcrystfel/src/reflist.c index 64154774..b3b9f85b 100644 --- a/libcrystfel/src/reflist.c +++ b/libcrystfel/src/reflist.c @@ -3,11 +3,11 @@ * * Fast reflection/peak list * - * 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: - * 2011-2012 Thomas White <taw@physics.org> + * 2011-2013 Thomas White <taw@physics.org> * * This file is part of CrystFEL. * @@ -70,6 +70,7 @@ struct _refldata { double r1; /* First excitation error */ double r2; /* Second excitation error */ double p; /* Partiality */ + double L; /* Lorentz factor */ int clamp1; /* Clamp status for r1 */ int clamp2; /* Clamp status for r2 */ @@ -388,7 +389,7 @@ void get_symmetric_indices(const Reflection *refl, * get_partiality: * @refl: A %Reflection * - * Returns: The partiality of the reflection. + * Returns: The partiality of the reflection. See get_lorentz(). **/ double get_partiality(const Reflection *refl) { @@ -397,6 +398,19 @@ double get_partiality(const Reflection *refl) /** + * get_lorentz: + * @refl: A %Reflection + * + * Returns: The Lorentz factor for the reflection. To "scale up" a partial + * reflection, divide by this multiplied by the partiality. + **/ +double get_lorentz(const Reflection *refl) +{ + return refl->data.L; +} + + +/** * get_intensity: * @refl: A %Reflection * @@ -602,10 +616,35 @@ void set_partial(Reflection *refl, double r1, double r2, double p, /** * set_intensity: * @refl: A %Reflection + * @p: The partiality for the reflection. + * + * Set the partiality for the reflection. See set_lorentz(). + **/ +void set_partiality(Reflection *refl, double p) +{ + refl->data.p = p; +} + +/** + * set_lorentz: + * @refl: A %Reflection + * @L: The Lorentz factor for the reflection. + * + * Set the Lorentz factor for the reflection. To "scale up" a partial + * reflection, divide by this multiplied by the partiality. + **/ +void set_lorentz(Reflection *refl, double L) +{ + refl->data.L = L; +} + + +/** + * set_intensity: + * @refl: A %Reflection * @intensity: The intensity for the reflection. * - * Set the intensity for the reflection. Note that retrieval is done with - * get_intensity(). + * Set the intensity for the reflection. **/ void set_intensity(Reflection *refl, double intensity) { diff --git a/libcrystfel/src/reflist.h b/libcrystfel/src/reflist.h index 246ef885..7969235c 100644 --- a/libcrystfel/src/reflist.h +++ b/libcrystfel/src/reflist.h @@ -3,11 +3,11 @@ * * Fast reflection/peak list * - * 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: - * 2011-2012 Thomas White <taw@physics.org> + * 2011-2013 Thomas White <taw@physics.org> * * This file is part of CrystFEL. * @@ -78,6 +78,7 @@ extern Reflection *next_found_refl(Reflection *refl); extern double get_excitation_error(const Reflection *refl); extern void get_detector_pos(const Reflection *refl, double *fs, double *ss); extern double get_partiality(const Reflection *refl); +extern double get_lorentz(const Reflection *refl); extern void get_indices(const Reflection *refl, signed int *h, signed int *k, signed int *l); extern void get_symmetric_indices(const Reflection *refl, @@ -100,6 +101,8 @@ 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_partiality(Reflection *refl, double p); +extern void set_lorentz(Reflection *refl, double L); extern void set_intensity(Reflection *refl, double intensity); extern void set_scalable(Reflection *refl, int scalable); extern void set_refinable(Reflection *refl, int refinable); diff --git a/libcrystfel/src/utils.h b/libcrystfel/src/utils.h index b75693db..1adb69e6 100644 --- a/libcrystfel/src/utils.h +++ b/libcrystfel/src/utils.h @@ -139,6 +139,11 @@ static inline double modulus(double x, double y, double z) return sqrt(x*x + y*y + z*z); } +static inline double modulus2d(double x, double y) +{ + return sqrt(x*x + y*y); +} + static inline double modulus_squared(double x, double y, double z) { return x*x + y*y + z*z; } @@ -165,6 +170,21 @@ static inline double angle_between(double x1, double y1, double z1, return acos(cosine); } +/* Answer in radians */ +static inline double angle_between_2d(double x1, double y1, + double x2, double y2) +{ + double mod1 = modulus2d(x1, y1); + double mod2 = modulus2d(x2, y2); + double cosine = (x1*x2 + y1*y2) / (mod1*mod2); + + /* Fix domain if outside due to rounding */ + if ( cosine > 1.0 ) cosine = 1.0; + if ( cosine < -1.0 ) cosine = -1.0; + + return acos(cosine); +} + static inline int within_tolerance(double a, double b, double percent) { double tol = fabs(a) * (percent/100.0); |