diff options
| -rw-r--r-- | src/geometry.c | 17 |
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"? */ |