Factorise dr/da part of gradient calculation for PR and prediction

3 files changed, 50 insertions, 213 deletions

diff --git a/src/post-refinement.c b/src/post-refinement.c
index 745e5bd7..a532cd07 100644
--- a/src/post-refinement.c
+++ b/src/post-refinement.c
@@ -54,6 +54,8 @@
 /* Maximum number of iterations of NLSq to do for each image per macrocycle. */
 #define MAX_CYCLES (10)
 
+/* Number of parameters to refine, in the order they appear in the enum */
+#define NUM_PARAMS (9)
 
 /* Returns dp(gauss)/dr at "r" */
 static double gaussian_fraction_gradient(double r, double R)
@@ -179,27 +181,17 @@ static double volume_fraction(double rlow, double rhigh, double pr,
  * of 'image'. */
 double p_gradient(Crystal *cr, int k, Reflection *refl, PartialityModel pmodel)
 {
-	double azi;
 	double glow, ghigh;
-	double asx, asy, asz;
-	double bsx, bsy, bsz;
-	double csx, csy, csz;
-	double xl, yl, zl;
-	double ds;
-	signed int hs, ks, ls;
 	double rlow, rhigh, p;
-	double philow, phihigh, phi;
-	double khigh, klow;
-	double tl, cet, cez;
 	struct image *image = crystal_get_image(cr);
 	double R = crystal_get_profile_radius(cr);
-	double D, psph;
 
 	get_partial(refl, &rlow, &rhigh, &p);
 
-	if ( k == REF_R ) {
+	if ( k == GPARAM_R ) {
 
 		double Rglow, Rghigh;
+		double D, psph;
 
 		D = rlow - rhigh;
 		psph = volume_fraction(rlow, rhigh, R, pmodel);
@@ -215,78 +207,21 @@ double p_gradient(Crystal *cr, int k, Reflection *refl, PartialityModel pmodel)
 	glow = partiality_gradient(rlow, R, pmodel, rlow, rhigh);
 	ghigh = partiality_gradient(rhigh, R, pmodel, rlow, rhigh);
 
-	get_symmetric_indices(refl, &hs, &ks, &ls);
-	ds = 2.0 * resolution(crystal_get_cell(cr), hs, ks, ls);
-
-	cell_get_reciprocal(crystal_get_cell(cr), &asx, &asy, &asz,
-	                                          &bsx, &bsy, &bsz,
-	                                          &csx, &csy, &csz);
-	xl = hs*asx + ks*bsx + ls*csx;
-	yl = hs*asy + ks*bsy + ls*csy;
-	zl = hs*asz + ks*bsz + ls*csz;
-
-	/* "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);
-
-	tl = sqrt(xl*xl + yl*yl);
+	if ( k == GPARAM_DIV ) {
 
-	cet = -sin(image->div/2.0) * klow;
-	cez = -cos(image->div/2.0) * klow;
-	philow = angle_between_2d(tl-cet, zl-cez, 0.0, 1.0);
-
-	cet = -sin(image->div/2.0) * khigh;
-	cez = -cos(image->div/2.0) * khigh;
-	phihigh = angle_between_2d(tl-cet, zl-cez, 0.0, 1.0);
-
-	/* Approximation: philow and phihigh are very similar */
-	phi = (philow + phihigh) / 2.0;
-
-	azi = atan2(yl, xl);
-
-	/* For many gradients, just multiply the above number by the gradient
-	 * of excitation error wrt whatever. */
-	switch ( k ) {
+		double D, psph, ds;
+		signed int hs, ks, ls;
 
-		/* Cell parameters and orientation */
-		case REF_ASX :
-		return - hs * sin(phi) * cos(azi) * (glow-ghigh);
-
-		case REF_BSX :
-		return - ks * sin(phi) * cos(azi) * (glow-ghigh);
-
-		case REF_CSX :
-		return - ls * sin(phi) * cos(azi) * (glow-ghigh);
-
-		case REF_ASY :
-		return - hs * sin(phi) * sin(azi) * (glow-ghigh);
-
-		case REF_BSY :
-		return - ks * sin(phi) * sin(azi) * (glow-ghigh);
-
-		case REF_CSY :
-		return - ls * sin(phi) * sin(azi) * (glow-ghigh);
-
-		case REF_ASZ :
-		return - hs * cos(phi) * (glow-ghigh);
-
-		case REF_BSZ :
-		return - ks * cos(phi) * (glow-ghigh);
-
-		case REF_CSZ :
-		return - ls * cos(phi) * (glow-ghigh);
-
-		case REF_DIV :
 		D = rlow - rhigh;
 		psph = volume_fraction(rlow, rhigh, R, pmodel);
+		get_symmetric_indices(refl, &hs, &ks, &ls);
+		ds = 2.0 * resolution(crystal_get_cell(cr), hs, ks, ls);
+
 		return (ds/2.0)*(glow+ghigh) - 4.0*R*psph*ds/(3.0*D*D);
 
 	}
 
-	ERROR("No gradient defined for parameter %i\n", k);
-	abort();
+	return r_gradient(crystal_get_cell(cr), k, refl, image) * (glow-ghigh);
 }
 
 
@@ -302,15 +237,15 @@ static void apply_cell_shift(UnitCell *cell, int k, double shift)
 
 	switch ( k )
 	{
-		case REF_ASX :  asx += shift;  break;
-		case REF_ASY :  asy += shift;  break;
-		case REF_ASZ :  asz += shift;  break;
-		case REF_BSX :  bsx += shift;  break;
-		case REF_BSY :  bsy += shift;  break;
-		case REF_BSZ :  bsz += shift;  break;
-		case REF_CSX :  csx += shift;  break;
-		case REF_CSY :  csy += shift;  break;
-		case REF_CSZ :  csz += shift;  break;
+		case GPARAM_ASX :  asx += shift;  break;
+		case GPARAM_ASY :  asy += shift;  break;
+		case GPARAM_ASZ :  asz += shift;  break;
+		case GPARAM_BSX :  bsx += shift;  break;
+		case GPARAM_BSY :  bsy += shift;  break;
+		case GPARAM_BSZ :  bsz += shift;  break;
+		case GPARAM_CSX :  csx += shift;  break;
+		case GPARAM_CSY :  csy += shift;  break;
+		case GPARAM_CSZ :  csz += shift;  break;
 	}
 
 	cell_set_reciprocal(cell, asx, asy, asz,
@@ -327,7 +262,7 @@ static void apply_shift(Crystal *cr, int k, double shift)
 
 	switch ( k ) {
 
-		case REF_DIV :
+		case GPARAM_DIV :
 		if ( isnan(shift) ) {
 			ERROR("NaN divergence shift\n");
 		} else {
@@ -336,21 +271,21 @@ static void apply_shift(Crystal *cr, int k, double shift)
 		}
 		break;
 
-		case REF_R :
+		case GPARAM_R :
 		t = crystal_get_profile_radius(cr);
 		t += shift;
 		crystal_set_profile_radius(cr, t);
 		break;
 
-		case REF_ASX :
-		case REF_ASY :
-		case REF_ASZ :
-		case REF_BSX :
-		case REF_BSY :
-		case REF_BSZ :
-		case REF_CSX :
-		case REF_CSY :
-		case REF_CSZ :
+		case GPARAM_ASX :
+		case GPARAM_ASY :
+		case GPARAM_ASZ :
+		case GPARAM_BSX :
+		case GPARAM_BSY :
+		case GPARAM_BSZ :
+		case GPARAM_CSX :
+		case GPARAM_CSY :
+		case GPARAM_CSZ :
 		apply_cell_shift(crystal_get_cell(cr), k, shift);
 		break;
 
diff --git a/src/post-refinement.h b/src/post-refinement.h
index a9c79ed6..caf5f9af 100644
--- a/src/post-refinement.h
+++ b/src/post-refinement.h
@@ -43,24 +43,6 @@
 #include "geometry.h"
 
 
-/* Refineable parameters.
- * Don't forget to add new things to backup_crystal() et al.! */
-enum {
-	REF_ASX,
-	REF_ASY,
-	REF_ASZ,
-	REF_BSX,
-	REF_BSY,
-	REF_BSZ,
-	REF_CSX,
-	REF_CSY,
-	REF_CSZ,
-	NUM_PARAMS,
-	REF_R,
-	REF_DIV,
-};
-
-
 struct prdata
 {
 	int refined;
diff --git a/src/predict-refine.c b/src/predict-refine.c
index 6af9211f..a2a4ce53 100644
--- a/src/predict-refine.c
+++ b/src/predict-refine.c
@@ -37,7 +37,7 @@
 #include <gsl/gsl_vector.h>
 
 #include "image.h"
-#include "post-refinement.h"
+#include "geometry.h"
 #include "cell-utils.h"
 
 
@@ -127,86 +127,6 @@ static void twod_mapping(double fs, double ss, double *px, double *py,
 }
 
 
-static double r_gradient(UnitCell *cell, int k, Reflection *refl,
-                         struct image *image)
-{
-	double azi;
-	double asx, asy, asz;
-	double bsx, bsy, bsz;
-	double csx, csy, csz;
-	double xl, yl, zl;
-	signed int hs, ks, ls;
-	double rlow, rhigh, p;
-	double philow, phihigh, phi;
-	double khigh, klow;
-	double tl, cet, cez;
-
-	get_partial(refl, &rlow, &rhigh, &p);
-
-	get_symmetric_indices(refl, &hs, &ks, &ls);
-
-	cell_get_reciprocal(cell, &asx, &asy, &asz,
-	                          &bsx, &bsy, &bsz,
-	                          &csx, &csy, &csz);
-	xl = hs*asx + ks*bsx + ls*csx;
-	yl = hs*asy + ks*bsy + ls*csy;
-	zl = hs*asz + ks*bsz + ls*csz;
-
-	/* "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);
-
-	tl = sqrt(xl*xl + yl*yl);
-
-	cet = -sin(image->div/2.0) * klow;
-	cez = -cos(image->div/2.0) * klow;
-	philow = angle_between_2d(tl-cet, zl-cez, 0.0, 1.0);
-
-	cet = -sin(image->div/2.0) * khigh;
-	cez = -cos(image->div/2.0) * khigh;
-	phihigh = angle_between_2d(tl-cet, zl-cez, 0.0, 1.0);
-
-	/* Approximation: philow and phihigh are very similar */
-	phi = (philow + phihigh) / 2.0;
-
-	azi = atan2(yl, xl);
-
-	switch ( k ) {
-
-		case REF_ASX :
-		return - hs * sin(phi) * cos(azi);
-
-		case REF_BSX :
-		return - ks * sin(phi) * cos(azi);
-
-		case REF_CSX :
-		return - ls * sin(phi) * cos(azi);
-
-		case REF_ASY :
-		return - hs * sin(phi) * sin(azi);
-
-		case REF_BSY :
-		return - ks * sin(phi) * sin(azi);
-
-		case REF_CSY :
-		return - ls * sin(phi) * sin(azi);
-
-		case REF_ASZ :
-		return - hs * cos(phi);
-
-		case REF_BSZ :
-		return - ks * cos(phi);
-
-		case REF_CSZ :
-		return - ls * cos(phi);
-
-	}
-
-	ERROR("No gradient defined for parameter %i\n", k);
-	abort();
-}
 
 
 /* Returns d(xh-xpk)/dP + d(yh-ypk)/dP, where P = any parameter */
@@ -230,31 +150,31 @@ static double x_gradient(int param, struct reflpeak *rp, struct detector *det,
 
 	switch ( param ) {
 
-		case REF_ASX :
+		case GPARAM_ASX :
 		return h * lambda * clen / cos(tt);
 
-		case REF_BSX :
+		case GPARAM_BSX :
 		return k * lambda * clen / cos(tt);
 
-		case REF_CSX :
+		case GPARAM_CSX :
 		return l * lambda * clen / cos(tt);
 
-		case REF_ASY :
+		case GPARAM_ASY :
 		return 0.0;
 
-		case REF_BSY :
+		case GPARAM_BSY :
 		return 0.0;
 
-		case REF_CSY :
+		case GPARAM_CSY :
 		return 0.0;
 
-		case REF_ASZ :
+		case GPARAM_ASZ :
 		return -h * lambda * clen * azf * sin(tt) / (cos(tt)*cos(tt));
 
-		case REF_BSZ :
+		case GPARAM_BSZ :
 		return -k * lambda * clen * azf * sin(tt) / (cos(tt)*cos(tt));
 
-		case REF_CSZ :
+		case GPARAM_CSZ :
 		return -l * lambda * clen * azf * sin(tt) / (cos(tt)*cos(tt));
 
 	}
@@ -285,31 +205,31 @@ static double y_gradient(int param, struct reflpeak *rp, struct detector *det,
 
 	switch ( param ) {
 
-		case REF_ASX :
+		case GPARAM_ASX :
 		return 0.0;
 
-		case REF_BSX :
+		case GPARAM_BSX :
 		return 0.0;
 
-		case REF_CSX :
+		case GPARAM_CSX :
 		return 0.0;
 
-		case REF_ASY :
+		case GPARAM_ASY :
 		return h * lambda * clen / cos(tt);
 
-		case REF_BSY :
+		case GPARAM_BSY :
 		return k * lambda * clen / cos(tt);
 
-		case REF_CSY :
+		case GPARAM_CSY :
 		return l * lambda * clen / cos(tt);
 
-		case REF_ASZ :
+		case GPARAM_ASZ :
 		return -h * lambda * clen * azf * sin(tt) / (cos(tt)*cos(tt));
 
-		case REF_BSZ :
+		case GPARAM_BSZ :
 		return -k * lambda * clen * azf * sin(tt) / (cos(tt)*cos(tt));
 
-		case REF_CSZ :
+		case GPARAM_CSZ :
 		return -l * lambda * clen * azf * sin(tt) / (cos(tt)*cos(tt));
 
 	}