aboutsummaryrefslogtreecommitdiff
path: root/src
diff options
context:
space:
mode:
authorThomas White <taw@physics.org>2010-10-28 15:07:43 +0200
committerThomas White <taw@physics.org>2012-02-22 15:27:04 +0100
commitea3422faf57e288c00bf971550601cb54773d07f (patch)
tree1d087e870d336159e0d1b23eedbcefe25f71c9b8 /src
parentf9c54e91c2a8e0c2e0a17a99bc49f05b45bc9ba7 (diff)
More geometry work (fixes a previous commit)
Diffstat (limited to 'src')
-rw-r--r--src/geometry.c73
-rw-r--r--src/templates.c3
2 files changed, 33 insertions, 43 deletions
diff --git a/src/geometry.c b/src/geometry.c
index 2b94e7ba..1f21c96f 100644
--- a/src/geometry.c
+++ b/src/geometry.c
@@ -76,28 +76,19 @@ static signed int locate_peak(double x, double y, double z, double lambda,
static double excitation_error(double xl, double yl, double zl,
- double ds_sq, double lambda)
+ double ds, double lambda)
{
- double tt;
- double a, b, c;
- double r1, r2;
- int n;
+ double tt, al;
+ double r;
tt = angle_between(0.0, 0.0, 1.0, xl, yl, zl+1.0/lambda);
+ al = M_PI_2 - asin(-zl/ds);
- a = 1.0;
- b = - cos(tt) * 2.0/lambda;
- c = pow(lambda, -2.0) - ds_sq;
+ r = ds * sin(al) / sin(tt);
- /* FIXME: I don't trust GSL's quadratic solver */
- n = gsl_poly_solve_quadratic(a, b, c, &r1, &r2);
- assert(n == 2);
+ r -= 1.0/lambda;
- r1 -= 1.0/lambda;
- r2 -= 1.0/lambda;
-
- if ( r1 > r2 ) return r2;
- return r1;
+ return r;
}
@@ -115,7 +106,7 @@ struct cpeak *find_intersections(struct image *image, UnitCell *cell,
double bandwidth = image->bw;
double divergence = image->div;
double lambda = image->lambda;
- const double profile_cutoff = 0.05e9; /* 0.1 nm^-1 */
+ const double profile_cutoff = 0.005e9; /* 0.1 nm^-1 */
double llow, lhigh; /* Wavelength */
cpeaks = malloc(sizeof(struct cpeak)*MAX_CPEAKS);
@@ -148,10 +139,10 @@ struct cpeak *find_intersections(struct image *image, UnitCell *cell,
double xl, yl, zl;
double ds, ds_sq;
double rlow, rhigh; /* "Excitation error" */
- signed int plow, phigh;
- double xdalow, xdahigh;
- double ydalow, ydahigh;
+ signed int p;
+ double lcen;
double xda, yda;
+ int close, straddled;
/* Ignore central beam */
if ( (h==0) && (k==0) && (l==0) ) continue;
@@ -168,32 +159,32 @@ struct cpeak *find_intersections(struct image *image, UnitCell *cell,
if ( ds > mres ) continue; /* Outside resolution range */
/* Calculate excitation errors */
- rlow = excitation_error(xl, yl, zl, ds_sq, llow);
- rhigh = excitation_error(xl, yl, zl, ds_sq, lhigh);
-
- if ( (fabs(rlow) > profile_cutoff)
- && (fabs(rhigh) > profile_cutoff) ) {
- /* Reflection is not close to Bragg at either of
- * the two extremes of wavelength, so skip it. */
- continue;
- }
+ rlow = excitation_error(xl, yl, zl, ds, llow);
+ rhigh = excitation_error(xl, yl, zl, ds, lhigh);
+
+ /* 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"? */
+ straddled = signbit(rlow) ^ signbit(rhigh);
+
+ /* Neither? Skip it. */
+ if ( !(close || straddled) ) continue;
+
+ lcen = -2.0*zl / ds_sq;
/* Locate peak on detector, and check it doesn't span panels */
- plow = locate_peak(xl, yl, zl, llow, image->det,
- &xdalow, &ydalow);
- if ( plow == -1 ) continue;
- phigh = locate_peak(xl, yl, zl, lhigh, image->det,
- &xdahigh, &ydahigh);
- if ( phigh == -1 ) continue;
- if ( plow != phigh ) continue;
-
- xda = xdahigh;
- yda = ydahigh;
+ p = locate_peak(xl, yl, zl, lcen, image->det, &xda, &yda);
+ if ( p == -1 ) continue;
+
cpeaks[np].h = h;
cpeaks[np].k = k;
cpeaks[np].l = l;
- cpeaks[np].x = xdahigh;
- cpeaks[np].y = ydahigh;
+ cpeaks[np].x = xda;
+ cpeaks[np].y = yda;
np++;
if ( output ) {
diff --git a/src/templates.c b/src/templates.c
index e9bd7b72..4f47c572 100644
--- a/src/templates.c
+++ b/src/templates.c
@@ -173,8 +173,7 @@ IndexingPrivate *generate_templates(UnitCell *cell, const char *filename,
cell_rot = rotate_cell(cell, omega, phi, 0.0);
- cpeaks = find_intersections(&image, cell_rot, 5.0e-3,
- 3.0/100.0, &n, 0);
+ cpeaks = find_intersections(&image, cell_rot, &n, 0);
if ( cpeaks == NULL ) {
ERROR("Template calculation failed.\n");
return NULL;