From d3b0933dc3293c724bb6013dad8e909bc61e80b0 Mon Sep 17 00:00:00 2001 From: Thomas White Date: Mon, 27 Jun 2011 11:40:20 +0200 Subject: Gradient must take wavelength into account --- src/post-refinement.c | 24 ++++++++++++++++++++---- 1 file changed, 20 insertions(+), 4 deletions(-) (limited to 'src/post-refinement.c') 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; -- cgit v1.2.3