aboutsummaryrefslogtreecommitdiff
diff options
context:
space:
mode:
Diffstat
-rw-r--r--src/facetron.c4
-rw-r--r--src/geometry.c224
-rw-r--r--src/geometry.h2
-rw-r--r--src/post-refinement.c6
-rw-r--r--src/templates.c2
5 files changed, 125 insertions, 113 deletions
diff --git a/src/facetron.c b/src/facetron.c
index 2bac7511..663a3236 100644
--- a/src/facetron.c
+++ b/src/facetron.c
@@ -101,7 +101,7 @@ static void refine_image(int mytask, void *tasks)
return;
}
- spots = find_intersections(image, image->indexed_cell, &n, 0);
+ spots = find_intersections(image, image->indexed_cell, &n, 0, NULL);
dev = +INFINITY;
i = 0;
do {
@@ -162,7 +162,7 @@ static void integrate_image(int mytask, void *tasks)
}
/* Figure out which spots should appear in this pattern */
- spots = find_intersections(image, image->indexed_cell, &n, 0);
+ spots = find_intersections(image, image->indexed_cell, &n, 0, NULL);
/* For each reflection, estimate the partiality */
for ( j=0; j<n; j++ ) {
diff --git a/src/geometry.c b/src/geometry.c
index 664817e0..09544120 100644
--- a/src/geometry.c
+++ b/src/geometry.c
@@ -117,8 +117,122 @@ static double partiality(double r1, double r2, double r)
}
+static int check_reflection(struct image *image, double mres, int output,
+ struct cpeak *cpeaks, int np,
+ signed int h, signed int k, signed int l,
+ double asx, double asy, double asz,
+ double bsx, double bsy, double bsz,
+ double csx, double csy, double csz)
+{
+ double xl, yl, zl;
+ double ds, ds_sq;
+ double rlow, rhigh; /* "Excitation error" */
+ signed int p; /* Panel number */
+ double xda, yda; /* Position on detector */
+ int close, inside;
+ double part; /* Partiality */
+ int clamp_low = 0;
+ int clamp_high = 0;
+ double bandwidth = image->bw;
+ double divergence = image->div;
+ double lambda = image->lambda;
+ double klow, kcen, khigh; /* Wavenumber */
+ /* Bounding sphere for the shape transform approximation */
+ const double profile_cutoff = 0.02e9; /* 0.02 nm^-1 */
+
+ /* "low" gives the largest Ewald sphere,
+ * "high" gives the smallest Ewald sphere. */
+ klow = 1.0/(lambda - lambda*bandwidth/2.0);
+ kcen = 1.0/lambda;
+ khigh = 1.0/(lambda + lambda*bandwidth/2.0);
+
+ /* Get the coordinates of the reciprocal lattice point */
+ zl = h*asz + k*bsz + l*csz;
+ /* Throw out if it's "in front" */
+ if ( zl > profile_cutoff ) return 0;
+ xl = h*asx + k*bsx + l*csx;
+ yl = h*asy + k*bsy + l*csy;
+
+ /* Calculate reciprocal lattice point modulus (and square) */
+ ds_sq = modulus_squared(xl, yl, zl); /* d*^2 */
+ ds = sqrt(ds_sq);
+ if ( ds > mres ) return 0; /* Outside resolution range */
+
+ /* Calculate excitation errors */
+ rlow = excitation_error(xl, yl, zl, ds, klow, -divergence);
+ rhigh = excitation_error(xl, yl, zl, ds, khigh, +divergence);
+
+ /* Is the reciprocal lattice point close to either extreme of
+ * the sphere, maybe just outside the "Ewald volume"? */
+ close = (fabs(rlow) < profile_cutoff)
+ || (fabs(rhigh) < profile_cutoff);
+
+ /* Is the reciprocal lattice point somewhere between the
+ * extremes of the sphere, i.e. inside the "Ewald volume"? */
+ inside = signbit(rlow) ^ signbit(rhigh);
+
+ /* Can't be both inside and close */
+ if ( inside ) close = 0;
+
+ /* Neither? Skip it. */
+ if ( !(close || inside) ) return 0;
+
+ /* If the "lower" Ewald sphere is a long way away, use the
+ * position at which the Ewald sphere would just touch the
+ * reflection. */
+ if ( rlow < -profile_cutoff ) {
+ rlow = -profile_cutoff;
+ clamp_low = -1;
+ }
+ if ( rlow > +profile_cutoff ) {
+ rlow = +profile_cutoff;
+ clamp_low = +1;
+ }
+ /* Likewise the "higher" Ewald sphere */
+ if ( rhigh < -profile_cutoff ) {
+ rhigh = -profile_cutoff;
+ clamp_high = -1;
+ }
+ if ( rhigh > +profile_cutoff ) {
+ rhigh = +profile_cutoff;
+ clamp_high = +1;
+ }
+ /* The six possible combinations of clamp_{low,high} (including
+ * zero) correspond to the six situations in Table 3 of Rossmann
+ * et al. (1979). */
+
+ /* Calculate partiality and reject if too small */
+ part = partiality(rlow, rhigh, image->profile_radius);
+ if ( part < 0.1 ) return 0;
+
+ /* Locate peak on detector. */
+ p = locate_peak(xl, yl, zl, kcen, image->det, &xda, &yda);
+ if ( p == -1 ) return 0;
+
+ /* Add peak to list */
+ cpeaks[np].h = h;
+ cpeaks[np].k = k;
+ cpeaks[np].l = l;
+ cpeaks[np].x = xda;
+ cpeaks[np].y = yda;
+ cpeaks[np].r1 = rlow;
+ cpeaks[np].r2 = rhigh;
+ cpeaks[np].p = part;
+ cpeaks[np].clamp1 = clamp_low;
+ cpeaks[np].clamp2 = clamp_high;
+ np++;
+
+ if ( output ) {
+ printf("%3i %3i %3i %6f (at %5.2f,%5.2f) %5.2f\n",
+ h, k, l, 0.0, xda, yda, part);
+ }
+
+ return 1;
+}
+
+
struct cpeak *find_intersections(struct image *image, UnitCell *cell,
- int *n, int output)
+ int *n, int output, struct cpeak *t)
{
double asx, asy, asz;
double bsx, bsy, bsz;
@@ -128,12 +242,6 @@ struct cpeak *find_intersections(struct image *image, UnitCell *cell,
int hmax, kmax, lmax;
double mres;
signed int h, k, l;
- double bandwidth = image->bw;
- double divergence = image->div;
- double lambda = image->lambda;
- double klow, kcen, khigh; /* Wavenumber */
- /* Bounding sphere for the shape transform approximation */
- const double profile_cutoff = 0.02e9; /* 0.02 nm^-1 */
cpeaks = malloc(sizeof(struct cpeak)*MAX_CPEAKS);
if ( cpeaks == NULL ) {
@@ -150,112 +258,14 @@ struct cpeak *find_intersections(struct image *image, UnitCell *cell,
kmax = mres / modulus(bsx, bsy, bsz);
lmax = mres / modulus(csx, csy, csz);
- /* "low" gives the largest Ewald sphere,
- * "high" gives the smallest Ewald sphere. */
- klow = 1.0/(lambda - lambda*bandwidth/2.0);
- kcen = 1.0/lambda;
- khigh = 1.0/(lambda + lambda*bandwidth/2.0);
-
for ( h=-hmax; h<hmax; h++ ) {
for ( k=-kmax; k<kmax; k++ ) {
for ( l=-lmax; l<lmax; l++ ) {
-
- double xl, yl, zl;
- double ds, ds_sq;
- double rlow, rhigh; /* "Excitation error" */
- signed int p; /* Panel number */
- double xda, yda; /* Position on detector */
- int close, inside;
- double part; /* Partiality */
- int clamp_low = 0;
- int clamp_high = 0;
-
/* Ignore central beam */
if ( (h==0) && (k==0) && (l==0) ) continue;
-
- /* Get the coordinates of the reciprocal lattice point */
- zl = h*asz + k*bsz + l*csz;
- /* Throw out if it's "in front" */
- if ( zl > profile_cutoff ) continue;
- xl = h*asx + k*bsx + l*csx;
- yl = h*asy + k*bsy + l*csy;
-
- /* Calculate reciprocal lattice point modulus (and square) */
- ds_sq = modulus_squared(xl, yl, zl); /* d*^2 */
- ds = sqrt(ds_sq);
- if ( ds > mres ) continue; /* Outside resolution range */
-
- /* Calculate excitation errors */
- rlow = excitation_error(xl, yl, zl, ds, klow, -divergence);
- rhigh = excitation_error(xl, yl, zl, ds, khigh, +divergence);
-
- /* Is the reciprocal lattice point close to either extreme of
- * the sphere, maybe just outside the "Ewald volume"? */
- close = (fabs(rlow) < profile_cutoff)
- || (fabs(rhigh) < profile_cutoff);
-
- /* Is the reciprocal lattice point somewhere between the
- * extremes of the sphere, i.e. inside the "Ewald volume"? */
- inside = signbit(rlow) ^ signbit(rhigh);
-
- /* Can't be both inside and close */
- if ( inside ) close = 0;
-
- /* Neither? Skip it. */
- if ( !(close || inside) ) continue;
-
- /* If the "lower" Ewald sphere is a long way away, use the
- * position at which the Ewald sphere would just touch the
- * reflection. */
- if ( rlow < -profile_cutoff ) {
- rlow = -profile_cutoff;
- clamp_low = -1;
- }
- if ( rlow > +profile_cutoff ) {
- rlow = +profile_cutoff;
- clamp_low = +1;
- }
- /* Likewise the "higher" Ewald sphere */
- if ( rhigh < -profile_cutoff ) {
- rhigh = -profile_cutoff;
- clamp_high = -1;
- }
- if ( rhigh > +profile_cutoff ) {
- rhigh = +profile_cutoff;
- clamp_high = +1;
- }
- /* The six possible combinations of clamp_{low,high} (including
- * zero) correspond to the six situations in Table 3 of Rossmann
- * et al. (1979). */
-
- /* Calculate partiality and reject if too small */
- part = partiality(rlow, rhigh, image->profile_radius);
- if ( part < 0.1 ) continue;
-
- /* Locate peak on detector. */
- p = locate_peak(xl, yl, zl, kcen, image->det, &xda, &yda);
- if ( p == -1 ) continue;
-
- /* Add peak to list */
- cpeaks[np].h = h;
- cpeaks[np].k = k;
- cpeaks[np].l = l;
- cpeaks[np].x = xda;
- cpeaks[np].y = yda;
- cpeaks[np].r1 = rlow;
- cpeaks[np].r2 = rhigh;
- cpeaks[np].p = part;
- cpeaks[np].clamp1 = clamp_low;
- cpeaks[np].clamp2 = clamp_high;
- np++;
-
- if ( output ) {
- printf("%3i %3i %3i %6f (at %5.2f,%5.2f) %5.2f\n",
- h, k, l, 0.0, xda, yda, part);
- }
-
+ np += check_reflection(image, mres, output, cpeaks, np, h, k, l,
+ asx,asy,asz,bsx,bsy,bsz,csx,csy,csz);
if ( np == MAX_CPEAKS ) goto out;
-
}
}
}
diff --git a/src/geometry.h b/src/geometry.h
index 0782a3e1..03fde52f 100644
--- a/src/geometry.h
+++ b/src/geometry.h
@@ -18,7 +18,7 @@
#endif
extern struct cpeak *find_intersections(struct image *image, UnitCell *cell,
- int *n, int output);
+ int *n, int output, struct cpeak *t);
extern double integrate_all(struct image *image, struct cpeak *cpeaks, int n);
diff --git a/src/post-refinement.c b/src/post-refinement.c
index e2fa6d50..c4ca03cc 100644
--- a/src/post-refinement.c
+++ b/src/post-refinement.c
@@ -165,6 +165,7 @@ double pr_iterate(struct image *image, double *i_full, const char *sym,
gsl_vector *shifts;
int h, param;
struct cpeak *spots = *pspots;
+ struct cpeak *spots_old;
M = gsl_matrix_calloc(NUM_PARAMS, NUM_PARAMS);
v = gsl_vector_calloc(NUM_PARAMS);
@@ -241,8 +242,9 @@ double pr_iterate(struct image *image, double *i_full, const char *sym,
gsl_vector_free(v);
gsl_vector_free(shifts);
- free(spots);
- spots = find_intersections(image, image->indexed_cell, n, 0);
+ spots_old = spots;
+ spots = find_intersections(image, image->indexed_cell, n, 0, spots_old);
+ free(spots_old);
*pspots = spots;
return mean_partial_dev(image, spots, *n, sym, i_full, NULL);
}
diff --git a/src/templates.c b/src/templates.c
index 4bd85a76..57434473 100644
--- a/src/templates.c
+++ b/src/templates.c
@@ -172,7 +172,7 @@ IndexingPrivate *generate_templates(UnitCell *cell, const char *filename,
cell_rot = rotate_cell(cell, omega, phi, 0.0);
- cpeaks = find_intersections(&image, cell_rot, &n, 0);
+ cpeaks = find_intersections(&image, cell_rot, &n, 0, NULL);
if ( cpeaks == NULL ) {
ERROR("Template calculation failed.\n");
return NULL;
generated by cgit v1.2.3 (git 2.25.1) at 2024-11-16 19:53:06 +0000