diff options
Diffstat (limited to 'src/geometry.c')
| -rw-r--r-- | src/geometry.c | 40 |
1 files changed, 37 insertions, 3 deletions
diff --git a/src/geometry.c b/src/geometry.c index 11020909..9db1f103 100644 --- a/src/geometry.c +++ b/src/geometry.c @@ -90,6 +90,29 @@ static double excitation_error(double xl, double yl, double zl, } +static double partiality(double r1, double r2, double r) +{ + double q1, q2; + double p, p1, p2; + + /* Calculate degrees of penetration */ + q1 = (r1 + r)/(2.0*r); + q2 = (r2 + r)/(2.0*r); + + 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; + + 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); + + p = fabs(p2 - p1); + + return p; +} + + struct cpeak *find_intersections(struct image *image, UnitCell *cell, int *n, int output) { @@ -104,8 +127,11 @@ 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.005e9; /* 0.1 nm^-1 */ double klow, kcen, khigh; /* Wavenumber */ + /* Bounding sphere for the shape transform approximation */ + const double profile_cutoff = 0.02e9; /* 0.02 nm^-1 */ + /* Actual radius of the profile */ + const double profile_radius = 0.005e9; cpeaks = malloc(sizeof(struct cpeak)*MAX_CPEAKS); if ( cpeaks == NULL ) { @@ -138,6 +164,7 @@ struct cpeak *find_intersections(struct image *image, UnitCell *cell, signed int p; /* Panel number */ double xda, yda; /* Position on detector */ int close, inside; + double part; /* Partiality */ /* Ignore central beam */ if ( (h==0) && (k==0) && (l==0) ) continue; @@ -187,16 +214,23 @@ struct cpeak *find_intersections(struct image *image, UnitCell *cell, p = locate_peak(xl, yl, zl, kcen, image->det, &xda, &yda); if ( p == -1 ) continue; + part = partiality(rlow, rhigh, profile_radius); + + if ( part < 0.1 ) continue; + 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; np++; if ( output ) { - printf("%3i %3i %3i %6f (at %5.2f,%5.2f) %9e %9e\n", - h, k, l, 0.0, xda, yda, rlow, rhigh); + printf("%3i %3i %3i %6f (at %5.2f,%5.2f) %5.2f\n", + h, k, l, 0.0, xda, yda, part); } if ( np == MAX_CPEAKS ) goto out; |