aboutsummaryrefslogtreecommitdiff
path: root/src
diff options
context:
space:
mode:
authorThomas White <taw@physics.org>2011-06-27 11:40:20 +0200
committerThomas White <taw@physics.org>2012-02-22 15:27:30 +0100
commitd3b0933dc3293c724bb6013dad8e909bc61e80b0 (patch)
treef6f020ecf843331090fc4eff3780f742715b5238 /src
parent2cfcbecc2effde939f10789949ab4c8f4a3e2cd9 (diff)
Gradient must take wavelength into account
Diffstat (limited to 'src')
-rw-r--r--src/post-refinement.c24
1 files changed, 20 insertions, 4 deletions
diff --git a/src/post-refinement.c b/src/post-refinement.c
index 10c24e0c..61ea4e20 100644
--- a/src/post-refinement.c
+++ b/src/post-refinement.c
@@ -73,7 +73,8 @@ static double partiality_rgradient(double r, double profile_radius)
/* Return the gradient of parameter 'k' given the current status of 'image'. */
double gradient(struct image *image, int k, Reflection *refl, double r)
{
- double ds, tt, azix, aziy;
+ double ds, azix, aziy;
+ double ttlow, tthigh, tt;
double nom, den;
double g;
double asx, asy, asz;
@@ -83,6 +84,7 @@ double gradient(struct image *image, int k, Reflection *refl, double r)
signed int hs, ks, ls;
double r1, r2, p;
int clamp_low, clamp_high;
+ double klow, khigh;
get_symmetric_indices(refl, &hs, &ks, &ls);
@@ -94,11 +96,25 @@ double gradient(struct image *image, int k, Reflection *refl, double r)
zl = hs*asz + ks*bsz + ls*csz;
ds = 2.0 * resolution(image->indexed_cell, hs, ks, ls);
- /* FIXME: Should use k value corresponding to clamp status here */
- tt = angle_between(0.0, 0.0, 1.0, xl, yl, zl+1.0/image->lambda);
+ get_partial(refl, &r1, &r2, &p, &clamp_low, &clamp_high);
+
+ klow = 1.0/(image->lambda - image->lambda*image->bw/2.0);
+ khigh = 1.0/(image->lambda + image->lambda*image->bw/2.0);
+ ttlow = angle_between(0.0, 0.0, 1.0, xl, yl, zl+klow);
+ tthigh = angle_between(0.0, 0.0, 1.0, xl, yl, zl+khigh);
+ if ( (clamp_low == 0) && (clamp_high == 0) ) {
+ tt = (ttlow+tthigh)/2.0;
+ } else if ( clamp_high == 0 ) {
+ tt = tthigh;
+ } else if ( clamp_low == 0 ) {
+ tt = ttlow;
+ } else {
+ tt = 0.0;
+ /* Gradient should come out as zero in this case */
+ }
+
azix = angle_between(1.0, 0.0, 0.0, xl, yl, 0.0);
aziy = angle_between(0.0, 1.0, 0.0, xl, yl, 0.0);
- get_partial(refl, &r1, &r2, &p, &clamp_low, &clamp_high);
/* Calculate the gradient of partiality wrt excitation error. */
g = 0.0;