aboutsummaryrefslogtreecommitdiff
path: root/libcrystfel/src/geometry.c
diff options
context:
space:
mode:
Diffstat (limited to 'libcrystfel/src/geometry.c')
-rw-r--r--libcrystfel/src/geometry.c110
1 files changed, 90 insertions, 20 deletions
diff --git a/libcrystfel/src/geometry.c b/libcrystfel/src/geometry.c
index 60fab488..704baa51 100644
--- a/libcrystfel/src/geometry.c
+++ b/libcrystfel/src/geometry.c
@@ -3,11 +3,11 @@
*
* Geometry of diffraction
*
- * Copyright © 2012 Deutsches Elektronen-Synchrotron DESY,
+ * Copyright © 2013 Deutsches Elektronen-Synchrotron DESY,
* a research centre of the Helmholtz Association.
*
* Authors:
- * 2010-2012 Thomas White <taw@physics.org>
+ * 2010-2013 Thomas White <taw@physics.org>
*
* This file is part of CrystFEL.
*
@@ -33,6 +33,7 @@
#include <stdlib.h>
#include <assert.h>
+#include <fenv.h>
#include "utils.h"
#include "cell.h"
@@ -43,6 +44,7 @@
#include "reflist.h"
#include "reflist-utils.h"
#include "symmetry.h"
+#include "geometry.h"
static signed int locate_peak(double x, double y, double z, double k,
@@ -117,7 +119,7 @@ static double partiality(double rlow, double rhigh, double r)
}
-static Reflection *check_reflection(struct image *image,
+static Reflection *check_reflection(struct image *image, Crystal *cryst,
signed int h, signed int k, signed int l,
double xl, double yl, double zl)
{
@@ -131,6 +133,9 @@ static Reflection *check_reflection(struct image *image,
double klow, khigh; /* Wavenumber */
Reflection *refl;
double cet, cez;
+ double pr;
+
+ pr = crystal_get_profile_radius(cryst);
/* "low" gives the largest Ewald sphere (wavelength short => k large)
* "high" gives the smallest Ewald sphere (wavelength long => k small)
@@ -152,8 +157,8 @@ static Reflection *check_reflection(struct image *image,
rlow = klow - distance(cet, cez, tl, zl); /* Loss of precision */
if ( (signbit(rlow) == signbit(rhigh))
- && (fabs(rlow) > image->profile_radius)
- && (fabs(rhigh) > image->profile_radius) ) return NULL;
+ && (fabs(rlow) > pr)
+ && (fabs(rhigh) > pr) ) return NULL;
/* If the "lower" Ewald sphere is a long way away, use the
* position at which the Ewald sphere would just touch the
@@ -164,26 +169,26 @@ static Reflection *check_reflection(struct image *image,
* et al. (1979).
*/
clamp_low = 0; clamp_high = 0;
- if ( rlow < -image->profile_radius ) {
- rlow = -image->profile_radius;
+ if ( rlow < -pr ) {
+ rlow = -pr;
clamp_low = -1;
}
- if ( rlow > +image->profile_radius ) {
- rlow = +image->profile_radius;
+ if ( rlow > +pr ) {
+ rlow = +pr;
clamp_low = +1;
}
- if ( rhigh < -image->profile_radius ) {
- rhigh = -image->profile_radius;
+ if ( rhigh < -pr ) {
+ rhigh = -pr;
clamp_high = -1;
}
- if ( rhigh > +image->profile_radius ) {
- rhigh = +image->profile_radius;
+ if ( rhigh > +pr ) {
+ rhigh = +pr;
clamp_high = +1;
}
assert(clamp_low >= clamp_high);
/* Calculate partiality */
- part = partiality(rlow, rhigh, image->profile_radius);
+ part = partiality(rlow, rhigh, pr);
/* Locate peak on detector. */
p = locate_peak(xl, yl, zl, 1.0/image->lambda, image->det, &xda, &yda);
@@ -205,7 +210,7 @@ static Reflection *check_reflection(struct image *image,
}
-RefList *find_intersections(struct image *image, UnitCell *cell)
+RefList *find_intersections(struct image *image, Crystal *cryst)
{
double ax, ay, az;
double bx, by, bz;
@@ -217,6 +222,10 @@ RefList *find_intersections(struct image *image, UnitCell *cell)
int hmax, kmax, lmax;
double mres;
signed int h, k, l;
+ UnitCell *cell;
+
+ cell = crystal_get_cell(cryst);
+ if ( cell == NULL ) return NULL;
reflections = reflist_new();
@@ -262,7 +271,7 @@ RefList *find_intersections(struct image *image, UnitCell *cell)
yl = h*asy + k*bsy + l*csy;
zl = h*asz + k*bsz + l*csz;
- refl = check_reflection(image, h, k, l, xl, yl, zl);
+ refl = check_reflection(image, cryst, h, k, l, xl, yl, zl);
if ( refl != NULL ) {
add_refl_to_list(refl, reflections);
@@ -276,8 +285,68 @@ RefList *find_intersections(struct image *image, UnitCell *cell)
}
+RefList *select_intersections(struct image *image, Crystal *cryst)
+{
+ double ax, ay, az;
+ double bx, by, bz;
+ double cx, cy, cz;
+ const double min_dist = 0.25;
+ RefList *list;
+ int i;
+
+ /* Round towards nearest */
+ fesetround(1);
+
+ /* Cell basis vectors for this image */
+ cell_get_cartesian(crystal_get_cell(cryst), &ax, &ay, &az,
+ &bx, &by, &bz, &cx, &cy, &cz);
+
+ list = reflist_new();
+ if ( list == NULL ) return NULL;
+
+ /* Loop over peaks, checking proximity to nearest reflection */
+ for ( i=0; i<image_feature_count(image->features); i++ ) {
+
+ struct imagefeature *f;
+ struct rvec q;
+ double h, k, l, hd, kd, ld;
+ double dsq;
+
+ f = image_get_feature(image->features, i);
+ if ( f == NULL ) continue;
+
+ /* Reciprocal space position of found peak */
+ q = get_q(image, f->fs, f->ss, NULL, 1.0/image->lambda);
+
+ /* Decimal and fractional Miller indices of nearest
+ * reciprocal lattice point */
+ hd = q.u * ax + q.v * ay + q.w * az;
+ kd = q.u * bx + q.v * by + q.w * bz;
+ ld = q.u * cx + q.v * cy + q.w * cz;
+ h = lrint(hd);
+ k = lrint(kd);
+ l = lrint(ld);
+
+ /* Check distance */
+ dsq = pow(h-hd, 2.0) + pow(k-kd, 2.0) + pow(l-ld, 2.0);
+
+ if ( sqrt(dsq) < min_dist ) {
+
+ Reflection *refl;
+
+ refl = add_refl(list, h, k, l);
+ set_detector_pos(refl, sqrt(dsq), f->fs, f->ss);
+
+ }
+
+ }
+
+ return list;
+}
+
+
/* Calculate partialities and apply them to the image's reflections */
-void update_partialities(struct image *image)
+void update_partialities(Crystal *cryst)
{
Reflection *refl;
RefListIterator *iter;
@@ -285,14 +354,15 @@ void update_partialities(struct image *image)
double asx, asy, asz;
double bsx, bsy, bsz;
double csx, csy, csz;
+ struct image *image = crystal_get_image(cryst);
- cell_get_reciprocal(image->indexed_cell, &asx, &asy, &asz,
+ cell_get_reciprocal(crystal_get_cell(cryst), &asx, &asy, &asz,
&bsx, &bsy, &bsz, &csx, &csy, &csz);
/* Scratch list to give check_reflection() something to add to */
predicted = reflist_new();
- for ( refl = first_refl(image->reflections, &iter);
+ for ( refl = first_refl(crystal_get_reflections(cryst), &iter);
refl != NULL;
refl = next_refl(refl, iter) )
{
@@ -309,7 +379,7 @@ void update_partialities(struct image *image)
yl = h*asy + k*bsy + l*csy;
zl = h*asz + k*bsz + l*csz;
- vals = check_reflection(image, h, k, l, xl, yl, zl);
+ vals = check_reflection(image, cryst, h, k, l, xl, yl, zl);
if ( vals == NULL ) {
set_redundancy(refl, 0);
generated by cgit v1.2.3 (git 2.25.1) at 2024-07-15 02:32:25 +0000