Fix post refinement maths and take clamping status into account

5 files changed, 78 insertions, 36 deletions

diff --git a/src/facetron.c b/src/facetron.c
index f561b554..2bac7511 100644
--- a/src/facetron.c
+++ b/src/facetron.c
@@ -467,6 +467,7 @@ int main(int argc, char *argv[])
 		images[i].det = det;
 		images[i].beam = beam;
 		images[i].osf = 1.0;
+		images[i].profile_radius = 0.005e9;
 
 		/* Muppet proofing */
 		images[i].data = NULL;
diff --git a/src/geometry.c b/src/geometry.c
index 0ef1f1c3..664817e0 100644
--- a/src/geometry.c
+++ b/src/geometry.c
@@ -103,26 +103,17 @@ static double excitation_error(double xl, double yl, double zl,
 static double partiality(double r1, double r2, double r)
 {
 	double q1, q2;
-	double p, p1, p2;
+	double p1, p2;
 
 	/* Calculate degrees of penetration */
 	q1 = (r1 + r)/(2.0*r);
 	q2 = (r2 + r)/(2.0*r);
 
-	/* Clamp */
-	if ( q1 > 1.0 ) q1 = 1.0;
-	if ( q1 < 0.0 ) q1 = 0.0;
-	if ( q2 > 1.0 ) q2 = 1.0;
-	if ( q2 < 0.0 ) q2 = 0.0;
-
 	/* Convert to partiality */
 	p1 = 3.0*pow(q1,2.0) - 2.0*pow(q1,3.0);
 	p2 = 3.0*pow(q2,2.0) - 2.0*pow(q2,3.0);
 
-	/* Input values may have been backwards */
-	p = fabs(p2 - p1);
-
-	return p;
+	return p2 - p1;
 }
 
 
@@ -143,8 +134,6 @@ struct cpeak *find_intersections(struct image *image, UnitCell *cell,
 	double klow, kcen, khigh;    /* Wavenumber */
 	/* Bounding sphere for the shape transform approximation */
 	const double profile_cutoff = 0.02e9;  /* 0.02 nm^-1 */
-	/* Actual radius of the profile */
-	const double profile_radius = 0.005e9;
 
 	cpeaks = malloc(sizeof(struct cpeak)*MAX_CPEAKS);
 	if ( cpeaks == NULL ) {
@@ -178,6 +167,8 @@ struct cpeak *find_intersections(struct image *image, UnitCell *cell,
 		double xda, yda;        /* Position on detector */
 		int close, inside;
 		double part;            /* Partiality */
+		int clamp_low = 0;
+		int clamp_high = 0;
 
 		/* Ignore central beam */
 		if ( (h==0) && (k==0) && (l==0) ) continue;
@@ -216,21 +207,36 @@ struct cpeak *find_intersections(struct image *image, UnitCell *cell,
 		/* If the "lower" Ewald sphere is a long way away, use the
 		 * position at which the Ewald sphere would just touch the
 		 * reflection. */
-		if ( rlow > profile_cutoff ) rlow = profile_cutoff;
-		if ( rlow < -profile_cutoff ) rlow = -profile_cutoff;
+		if ( rlow < -profile_cutoff ) {
+			rlow = -profile_cutoff;
+			clamp_low = -1;
+		}
+		if ( rlow > +profile_cutoff ) {
+			rlow = +profile_cutoff;
+			clamp_low = +1;
+		}
+		/* Likewise the "higher" Ewald sphere */
+		if ( rhigh < -profile_cutoff ) {
+			rhigh = -profile_cutoff;
+			clamp_high = -1;
+		}
+		if ( rhigh > +profile_cutoff ) {
+			rhigh = +profile_cutoff;
+			clamp_high = +1;
+		}
+		/* The six possible combinations of clamp_{low,high} (including
+		 * zero) correspond to the six situations in Table 3 of Rossmann
+		 * et al. (1979). */
 
-		/* As above, other side. */
-		if ( rhigh > profile_cutoff ) rhigh = profile_cutoff;
-		if ( rhigh < -profile_cutoff ) rhigh = -profile_cutoff;
+		/* Calculate partiality and reject if too small */
+		part = partiality(rlow, rhigh, image->profile_radius);
+		if ( part < 0.1 ) continue;
 
 		/* Locate peak on detector. */
 		p = locate_peak(xl, yl, zl, kcen, image->det, &xda, &yda);
 		if ( p == -1 ) continue;
 
-		part = partiality(rlow, rhigh, profile_radius);
-
-		if ( part < 0.1 ) continue;
-
+		/* Add peak to list */
 		cpeaks[np].h = h;
 		cpeaks[np].k = k;
 		cpeaks[np].l = l;
@@ -239,6 +245,8 @@ struct cpeak *find_intersections(struct image *image, UnitCell *cell,
 		cpeaks[np].r1 = rlow;
 		cpeaks[np].r2 = rhigh;
 		cpeaks[np].p = part;
+		cpeaks[np].clamp1 = clamp_low;
+		cpeaks[np].clamp2 = clamp_high;
 		np++;
 
 		if ( output ) {
diff --git a/src/image.h b/src/image.h
index ae2f3179..0d2f6c2e 100644
--- a/src/image.h
+++ b/src/image.h
@@ -66,6 +66,8 @@ struct cpeak
 	double r1;  /* First excitation error */
 	double r2;  /* Second excitation error */
 	double p;   /* Partiality */
+	int clamp1; /* Clamp status for r1 */
+	int clamp2; /* Clamp status for r2 */
 
 	/* Location in image */
 	int x;
@@ -102,6 +104,7 @@ struct image {
 	int                     f0_available;  /* 0 if f0 wasn't available
 	                                        * from the input. */
 	double                  osf;           /* Overall scaling factor */
+	double                  profile_radius; /* Radius of reflection */
 
 	int                     width;
 	int                     height;
diff --git a/src/post-refinement.c b/src/post-refinement.c
index 5ca0dc98..e2fa6d50 100644
--- a/src/post-refinement.c
+++ b/src/post-refinement.c
@@ -28,12 +28,33 @@
 #include "geometry.h"
 
 
+/* Returns dp/dr at "r" */
+static double partiality_gradient(double r, double profile_radius)
+{
+	double q, dpdq, dqdr;
+
+	/* Calculate degree of penetration */
+	q = (r + profile_radius)/(2.0*profile_radius);
+
+	/* dp/dq */
+	dpdq = 6.0*(q-pow(q, 2.0));
+
+	/* dq/dr */
+	dqdr = 1.0 / (2.0*profile_radius);
+
+	return dpdq * dqdr;
+}
+
+
 /* Return the gradient of parameter 'k' given the current status of 'image'. */
 double gradient(struct image *image, int k,
-                       struct cpeak spot, double I_partial)
+                struct cpeak spot, double I_partial, double r)
 {
 	double ds;
 	double nom, den;
+	double r1g = 0.0;
+	double r2g = 0.0;
+	double g = 0.0;
 
 	ds = 2.0 * resolution(image->indexed_cell, spot.h, spot.k, spot.l);
 
@@ -43,9 +64,20 @@ double gradient(struct image *image, int k,
 		return I_partial;
 
 	case REF_DIV :
-		nom = sqrt(2.0) * ds * sin(image->div);
-		den = sqrt(1.0 - cos(image->div));
-		return nom/den;
+		/* Calculate the gradient of r1 and r2 wrt divergence */
+		if ( spot.clamp1 == 0 ) {
+			nom = sqrt(2.0) * ds * sin(image->div/2.0);
+			den = sqrt(1.0 - cos(image->div/2.0));
+			r1g = nom/den;
+			g += r1g * partiality_gradient(spot.r1, r);
+		}
+		if ( spot.clamp2 == 0 ) {
+			nom = sqrt(2.0) * ds * sin(image->div/2.0);
+			den = sqrt(1.0 - cos(image->div/2.0));
+			r2g = nom/den;
+			g += r2g * partiality_gradient(spot.r2, r);
+		}
+		return g;
 
 	}
 
@@ -171,7 +203,7 @@ double pr_iterate(struct image *image, double *i_full, const char *sym,
 		for ( k=0; k<NUM_PARAMS; k++ ) {
 
 			int g;
-			double v_c;
+			double v_c, gr;
 
 			for ( g=0; g<NUM_PARAMS; g++ ) {
 
@@ -179,16 +211,19 @@ double pr_iterate(struct image *image, double *i_full, const char *sym,
 
 				M_curr = gsl_matrix_get(M, g, k);
 
-				M_c = gradient(image, g, spots[h], I_partial)
-				    * gradient(image, k, spots[h], I_partial);
+				M_c = gradient(image, g, spots[h], I_partial,
+				               image->profile_radius)
+				    * gradient(image, k, spots[h], I_partial,
+				               image->profile_radius);
 				M_c *= pow(I_full, 2.0);
 
 				gsl_matrix_set(M, g, k, M_curr + M_c);
 
 			}
 
-			v_c = delta_I * I_full * gradient(image, k, spots[h],
-			                                  I_partial);
+			gr = gradient(image, k, spots[h], I_partial,
+			              image->profile_radius);
+			v_c = delta_I * I_full * gr;
 			gsl_vector_set(v, k, v_c);
 
 		}
diff --git a/src/post-refinement.h b/src/post-refinement.h
index 1ef7a1fd..8dea914d 100644
--- a/src/post-refinement.h
+++ b/src/post-refinement.h
@@ -30,11 +30,6 @@ enum {
 	NUM_PARAMS
 };
 
-
-/* Return the gradient of parameter 'k' given the current status of 'image'. */
-double gradient(struct image *image, int k, struct cpeak spot,
-                double I_partial);
-
 /* Apply the given shift to the 'k'th parameter of 'image'. */
 void apply_shift(struct image *image, int k, double shift);