diff options
| -rw-r--r-- | libcrystfel/src/utils.h | 20 | ||||
| -rw-r--r-- | src/post-refinement.c | 69 | ||||
| -rw-r--r-- | tests/pr_gradient_check.c | 3 |
3 files changed, 60 insertions, 32 deletions
diff --git a/libcrystfel/src/utils.h b/libcrystfel/src/utils.h index b75693db..1adb69e6 100644 --- a/libcrystfel/src/utils.h +++ b/libcrystfel/src/utils.h @@ -139,6 +139,11 @@ static inline double modulus(double x, double y, double z) return sqrt(x*x + y*y + z*z); } +static inline double modulus2d(double x, double y) +{ + return sqrt(x*x + y*y); +} + static inline double modulus_squared(double x, double y, double z) { return x*x + y*y + z*z; } @@ -165,6 +170,21 @@ static inline double angle_between(double x1, double y1, double z1, return acos(cosine); } +/* Answer in radians */ +static inline double angle_between_2d(double x1, double y1, + double x2, double y2) +{ + double mod1 = modulus2d(x1, y1); + double mod2 = modulus2d(x2, y2); + double cosine = (x1*x2 + y1*y2) / (mod1*mod2); + + /* Fix domain if outside due to rounding */ + if ( cosine > 1.0 ) cosine = 1.0; + if ( cosine < -1.0 ) cosine = -1.0; + + return acos(cosine); +} + static inline int within_tolerance(double a, double b, double percent) { double tol = fabs(a) * (percent/100.0); diff --git a/src/post-refinement.c b/src/post-refinement.c index f6f4caa5..9e63736f 100644 --- a/src/post-refinement.c +++ b/src/post-refinement.c @@ -90,8 +90,7 @@ static double partiality_rgradient(double r, double profile_radius) /* Return the gradient of parameter 'k' given the current status of 'image'. */ double gradient(Crystal *cr, int k, Reflection *refl, PartialityModel pmodel) { - double ds, azix, aziy; - double ttlow, tthigh, tt; + double ds, azi; double nom, den; double g; double asx, asy, asz; @@ -99,9 +98,11 @@ double gradient(Crystal *cr, int k, Reflection *refl, PartialityModel pmodel) double csx, csy, csz; double xl, yl, zl; signed int hs, ks, ls; - double r1, r2, p; + double rlow, rhigh, p; int clamp_low, clamp_high; - double klow, khigh; + double philow, phihigh, phi; + double khigh, klow; + double tl, cet, cez; double gr; struct image *image = crystal_get_image(cr); double r = crystal_get_profile_radius(cr); @@ -116,33 +117,37 @@ double gradient(Crystal *cr, int k, Reflection *refl, PartialityModel pmodel) zl = hs*asz + ks*bsz + ls*csz; ds = 2.0 * resolution(crystal_get_cell(cr), hs, ks, ls); - get_partial(refl, &r1, &r2, &p, &clamp_low, &clamp_high); + get_partial(refl, &rlow, &rhigh, &p, &clamp_low, &clamp_high); + /* "low" gives the largest Ewald sphere (wavelength short => k large) + * "high" gives the smallest Ewald sphere (wavelength long => k small) + */ 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 + image->div; - } else if ( clamp_low == 0 ) { - tt = ttlow - image->div; - } 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); + tl = sqrt(xl*xl + yl*yl); + ds = modulus(xl, yl, zl); + + cet = -sin(image->div/2.0) * klow; + cez = -cos(image->div/2.0) * klow; + philow = M_PI_2 - angle_between_2d(tl-cet, zl-cez, 1.0, 0.0); + + cet = -sin(image->div/2.0) * khigh; + cez = -cos(image->div/2.0) * khigh; + phihigh = M_PI_2 - angle_between_2d(tl-cet, zl-cez, 1.0, 0.0); + + /* Approximation: philow and phihigh are very similar */ + phi = (philow + phihigh) / 2.0; + + azi = atan2(yl, xl); /* Calculate the gradient of partiality wrt excitation error. */ g = 0.0; if ( clamp_low == 0 ) { - g -= partiality_gradient(r1, r); + g -= partiality_gradient(rlow, r); } if ( clamp_high == 0 ) { - g += partiality_gradient(r2, r); + g += partiality_gradient(rhigh, r); } /* For many gradients, just multiply the above number by the gradient @@ -167,40 +172,40 @@ double gradient(Crystal *cr, int k, Reflection *refl, PartialityModel pmodel) case REF_R : g = 0.0; if ( clamp_low == 0 ) { - g += partiality_rgradient(r1, r); + g += partiality_rgradient(rlow, r); } if ( clamp_high == 0 ) { - g += partiality_rgradient(r2, r); + g += partiality_rgradient(rhigh, r); } return g; /* Cell parameters and orientation */ case REF_ASX : - return hs * sin(tt) * cos(azix) * g; + return hs * sin(phi) * cos(azi) * g; case REF_BSX : - return ks * sin(tt) * cos(azix) * g; + return ks * sin(phi) * cos(azi) * g; case REF_CSX : - return ls * sin(tt) * cos(azix) * g; + return ls * sin(phi) * cos(azi) * g; case REF_ASY : - return hs * sin(tt) * cos(aziy) * g; + return hs * sin(phi) * sin(azi) * g; case REF_BSY : - return ks * sin(tt) * cos(aziy) * g; + return ks * sin(phi) * sin(azi) * g; case REF_CSY : - return ls * sin(tt) * cos(aziy) * g; + return ls * sin(phi) * sin(azi) * g; case REF_ASZ : - return hs * cos(tt) * g; + return hs * cos(phi) * g; case REF_BSZ : - return ks * cos(tt) * g; + return ks * cos(phi) * g; case REF_CSZ : - return ls * cos(tt) * g; + return ls * cos(phi) * g; } diff --git a/tests/pr_gradient_check.c b/tests/pr_gradient_check.c index cf11d18c..d5344801 100644 --- a/tests/pr_gradient_check.c +++ b/tests/pr_gradient_check.c @@ -425,6 +425,9 @@ int main(int argc, char *argv[]) val = test_gradients(cr, incr_val, REF_R, "R", pmodel); if ( val > 0.1 ) fail = 1; + //incr_val = incr_frac * image.profile_radius; + //val += test_gradients(&image, incr_val, REF_R, "rad"); + incr_val = incr_frac * ax; val = test_gradients(cr, incr_val, REF_ASX, "ax*", pmodel); if ( val > 0.1 ) fail = 1; |