aboutsummaryrefslogtreecommitdiff
path: root/src
diff options
context:
space:
mode:
authorThomas White <taw@physics.org>2013-03-08 10:49:49 +0100
committerThomas White <taw@physics.org>2013-04-17 17:33:48 +0200
commit3a348d3f7e2586440590747c1b921b8ce6dbc0e1 (patch)
treeb7f357e7c68b5f549491769bd194826eeb1dc64b /src
parentf8de3ad3f3410b95180f569d99b9b1f9bd9523ab (diff)
Work on post refinement
Brought across from "tom/pr" Conflicts: src/indexamajig.c src/post-refinement.c tests/pr_gradient_check.c
Diffstat (limited to 'src')
-rw-r--r--src/post-refinement.c69
1 files changed, 37 insertions, 32 deletions
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;
}