aboutsummaryrefslogtreecommitdiff
path: root/src
diff options
context:
space:
mode:
authorThomas White <taw@physics.org>2010-10-29 16:44:05 +0200
committerThomas White <taw@physics.org>2012-02-22 15:27:04 +0100
commit82284b116dab90f697970f1a2b1e798f352d7bde (patch)
tree3bd9a244fcdc2af1c739d80130af53e9696724bf /src
parentcc23b6eb31f8c4f59be3acf9f1bc8c51daf25e55 (diff)
Satisfactory peak location
Diffstat (limited to 'src')
-rw-r--r--src/geometry.c59
1 files changed, 33 insertions, 26 deletions
diff --git a/src/geometry.c b/src/geometry.c
index 1f21c96f..11020909 100644
--- a/src/geometry.c
+++ b/src/geometry.c
@@ -28,12 +28,12 @@
#define MAX_CPEAKS (256 * 256)
-static signed int locate_peak(double x, double y, double z, double lambda,
+static signed int locate_peak(double x, double y, double z, double k,
struct detector *det, double *xdap, double *ydap)
{
int p;
signed int found = -1;
- const double den = 1.0/lambda + z;
+ const double den = k + z;
for ( p=0; p<det->n_panels; p++ ) {
@@ -76,17 +76,15 @@ 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, double lambda)
+ double ds, double k)
{
double tt, al;
double r;
- tt = angle_between(0.0, 0.0, 1.0, xl, yl, zl+1.0/lambda);
+ tt = angle_between(0.0, 0.0, 1.0, xl, yl, zl+k);
al = M_PI_2 - asin(-zl/ds);
- r = ds * sin(al) / sin(tt);
-
- r -= 1.0/lambda;
+ r = ( ds * sin(al) / sin(tt) ) - k;
return r;
}
@@ -107,7 +105,7 @@ struct cpeak *find_intersections(struct image *image, UnitCell *cell,
double divergence = image->div;
double lambda = image->lambda;
const double profile_cutoff = 0.005e9; /* 0.1 nm^-1 */
- double llow, lhigh; /* Wavelength */
+ double klow, kcen, khigh; /* Wavenumber */
cpeaks = malloc(sizeof(struct cpeak)*MAX_CPEAKS);
if ( cpeaks == NULL ) {
@@ -119,9 +117,6 @@ struct cpeak *find_intersections(struct image *image, UnitCell *cell,
&bsx, &bsy, &bsz,
&csx, &csy, &csz);
- /* FIXME: Account for left-handed indexing */
- asz = -asz; bsz = -bsz; csz = -csz;
-
mres = 1.0 / 8.0e-10; /* 8 Angstroms */
hmax = mres / modulus(asx, asy, asz);
kmax = mres / modulus(bsx, bsy, bsz);
@@ -129,8 +124,9 @@ struct cpeak *find_intersections(struct image *image, UnitCell *cell,
/* "low" gives the largest Ewald sphere,
* "high" gives the smallest Ewald sphere. */
- llow = lambda - lambda*bandwidth/2.0;
- lhigh = lambda + lambda*bandwidth/2.0;
+ 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++ ) {
@@ -138,18 +134,18 @@ struct cpeak *find_intersections(struct image *image, UnitCell *cell,
double xl, yl, zl;
double ds, ds_sq;
- double rlow, rhigh; /* "Excitation error" */
- signed int p;
- double lcen;
- double xda, yda;
- int close, straddled;
+ double rlow, rhigh; /* "Excitation error" */
+ signed int p; /* Panel number */
+ double xda, yda; /* Position on detector */
+ int close, inside;
/* 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;
- if ( zl > 0.0 ) continue; /* Throw out if it's "in front" */
+ /* 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;
@@ -159,8 +155,8 @@ 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, llow);
- rhigh = excitation_error(xl, yl, zl, ds, lhigh);
+ rlow = excitation_error(xl, yl, zl, ds, klow);
+ rhigh = excitation_error(xl, yl, zl, ds, khigh);
/* Is the reciprocal lattice point close to either extreme of
* the sphere, maybe just outside the "Ewald volume"? */
@@ -169,15 +165,26 @@ struct cpeak *find_intersections(struct image *image, UnitCell *cell,
/* Is the reciprocal lattice point somewhere between the
* extremes of the sphere, i.e. inside the "Ewald volume"? */
- straddled = signbit(rlow) ^ signbit(rhigh);
+ inside = signbit(rlow) ^ signbit(rhigh);
+
+ /* Can't be both inside and close */
+ if ( inside ) close = 0;
/* Neither? Skip it. */
- if ( !(close || straddled) ) continue;
+ 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;
+ if ( rlow < -profile_cutoff ) rlow = -profile_cutoff;
- lcen = -2.0*zl / ds_sq;
+ /* As above, other side. */
+ if ( rhigh > profile_cutoff ) rhigh = profile_cutoff;
+ if ( rhigh < -profile_cutoff ) rhigh = -profile_cutoff;
- /* Locate peak on detector, and check it doesn't span panels */
- p = locate_peak(xl, yl, zl, lcen, image->det, &xda, &yda);
+ /* Locate peak on detector. */
+ p = locate_peak(xl, yl, zl, kcen, image->det, &xda, &yda);
if ( p == -1 ) continue;
cpeaks[np].h = h;