Work on post refinement

Brought across from "tom/pr"

Conflicts:
	src/indexamajig.c
	src/post-refinement.c
	tests/pr_gradient_check.c

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;
 
 	}