aboutsummaryrefslogtreecommitdiff
path: root/src
diff options
context:
space:
mode:
authorThomas White <taw@physics.org>2010-11-18 13:38:46 +0100
committerThomas White <taw@physics.org>2012-02-22 15:27:06 +0100
commit763a64fc497458e1920da555d8f5d4c831b21fc4 (patch)
treec536a7c906c5ace0cfe834f629f460e1aed4cd96 /src
parentd6ae8065f7c422300dba3b6256f4ace53b39a092 (diff)
Add divergence
Diffstat (limited to 'src')
-rw-r--r--src/geometry.c17
1 files changed, 14 insertions, 3 deletions
diff --git a/src/geometry.c b/src/geometry.c
index 9db1f103..e391715c 100644
--- a/src/geometry.c
+++ b/src/geometry.c
@@ -76,16 +76,24 @@ static signed int locate_peak(double x, double y, double z, double k,
static double excitation_error(double xl, double yl, double zl,
- double ds, double k)
+ double ds, double k, double divergence)
{
double tt, al;
double r;
+ double delta;
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) ) - k;
+ delta = sqrt(2.0 * pow(ds, 2.0) * (1-cos(divergence)));
+ if ( divergence > 0.0 ) {
+ r += delta;
+ } else {
+ r -= delta;
+ }
+
return r;
}
@@ -99,14 +107,17 @@ static double partiality(double r1, double r2, double r)
q1 = (r1 + r)/(2.0*r);
q2 = (r2 + r)/(2.0*r);
+ /* Clamp */
if ( q1 > 1.0 ) q1 = 1.0;
if ( q1 < 0.0 ) q1 = 0.0;
if ( q2 > 1.0 ) q2 = 1.0;
if ( q2 < 0.0 ) q2 = 0.0;
+ /* Convert to partiality */
p1 = 3.0*pow(q1,2.0) - 2.0*pow(q1,3.0);
p2 = 3.0*pow(q2,2.0) - 2.0*pow(q2,3.0);
+ /* Input values may have been backwards */
p = fabs(p2 - p1);
return p;
@@ -182,8 +193,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, klow);
- rhigh = excitation_error(xl, yl, zl, ds, khigh);
+ 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"? */