Move gradients to predict-refine

Makes sense for the _dev and _gradient functions to be in the same
place.,

5 files changed, 234 insertions, 230 deletions

diff --git a/libcrystfel/src/crystfel-mille.h b/libcrystfel/src/crystfel-mille.h
index 00d0415d..bdd66f37 100644
--- a/libcrystfel/src/crystfel-mille.h
+++ b/libcrystfel/src/crystfel-mille.h
@@ -34,7 +34,6 @@ typedef void *Mille;
 #include "cell.h"
 #include "image.h"
 #include "predict-refine.h"
-#include "geometry.h"
 
 /**
  * \file crystfel-mille.h
diff --git a/libcrystfel/src/geometry.c b/libcrystfel/src/geometry.c
index 5289856b..9dc471ee 100644
--- a/libcrystfel/src/geometry.c
+++ b/libcrystfel/src/geometry.c
@@ -429,69 +429,6 @@ static Reflection *check_reflection(struct image *image, Crystal *cryst,
 }
 
 
-double r_gradient(UnitCell *cell, int k, Reflection *refl, struct image *image)
-{
-	double asx, asy, asz;
-	double bsx, bsy, bsz;
-	double csx, csy, csz;
-	double xl, yl, zl;
-	signed int hs, ks, ls;
-	double tl, phi, azi;
-
-	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;
-
-	tl = sqrt(xl*xl + yl*yl);
-	phi = angle_between_2d(tl, zl+1.0/image->lambda, 0.0, 1.0); /* 2theta */
-	azi = atan2(yl, xl); /* azimuth */
-
-	switch ( k ) {
-
-		case GPARAM_ASX :
-		return - hs * sin(phi) * cos(azi);
-
-		case GPARAM_BSX :
-		return - ks * sin(phi) * cos(azi);
-
-		case GPARAM_CSX :
-		return - ls * sin(phi) * cos(azi);
-
-		case GPARAM_ASY :
-		return - hs * sin(phi) * sin(azi);
-
-		case GPARAM_BSY :
-		return - ks * sin(phi) * sin(azi);
-
-		case GPARAM_CSY :
-		return - ls * sin(phi) * sin(azi);
-
-		case GPARAM_ASZ :
-		return - hs * cos(phi);
-
-		case GPARAM_BSZ :
-		return - ks * cos(phi);
-
-		case GPARAM_CSZ :
-		return - ls * cos(phi);
-
-		case GPARAM_DETX :
-		case GPARAM_DETY :
-		case GPARAM_CLEN :
-		return 0.0;
-
-	}
-
-	ERROR("No r gradient defined for parameter %i\n", k);
-	abort();
-}
-
-
 /**
  * \param cryst: A \ref Crystal
  * \param max_res: Maximum resolution to predict to (m^-1)
@@ -1084,132 +1021,3 @@ void polarisation_correction(RefList *list, UnitCell *cell,
 		set_esd_intensity(refl, sigma / pol);
 	}
 }
-
-
-/* Returns dfs_refl/dP, where P = any parameter *
- * fs_refl is the fast scan position of refl in panel 'p',
- * in metres (not pixels) */
-double fs_gradient(int param, Reflection *refl, UnitCell *cell,
-                   struct detgeom_panel *p)
-{
-	signed int h, k, l;
-	double xl, zl, kpred;
-	double asx, asy, asz, bsx, bsy, bsz, csx, csy, csz;
-
-	get_indices(refl, &h, &k, &l);
-	kpred = get_kpred(refl);
-	cell_get_reciprocal(cell, &asx, &asy, &asz,
-	                          &bsx, &bsy, &bsz,
-	                          &csx, &csy, &csz);
-	xl = h*asx + k*bsx + l*csx;
-	zl = h*asz + k*bsz + l*csz;
-
-	switch ( param ) {
-
-		case GPARAM_ASX :
-		return h * p->cnz * p->pixel_pitch / (kpred + zl);
-
-		case GPARAM_BSX :
-		return k * p->cnz * p->pixel_pitch / (kpred + zl);
-
-		case GPARAM_CSX :
-		return l * p->cnz * p->pixel_pitch / (kpred + zl);
-
-		case GPARAM_ASY :
-		return 0.0;
-
-		case GPARAM_BSY :
-		return 0.0;
-
-		case GPARAM_CSY :
-		return 0.0;
-
-		case GPARAM_ASZ :
-		return -h * xl * p->cnz * p->pixel_pitch / (kpred*kpred + 2.0*kpred*zl + zl*zl);
-
-		case GPARAM_BSZ :
-		return -k * xl * p->cnz * p->pixel_pitch / (kpred*kpred + 2.0*kpred*zl + zl*zl);
-
-		case GPARAM_CSZ :
-		return -l * xl * p->cnz * p->pixel_pitch / (kpred*kpred + 2.0*kpred*zl + zl*zl);
-
-		case GPARAM_DETX :
-		return -1;
-
-		case GPARAM_DETY :
-		return 0;
-
-		case GPARAM_CLEN :
-		return -xl / (kpred+zl);
-
-		default:
-		ERROR("Positional gradient requested for parameter %i?\n", param);
-		abort();
-
-	}
-
-}
-
-
-/* Returns dss_refl/dP, where P = any parameter
- * ss_refl is the slow scan position of refl in panel 'p',
- * in metres (not pixels) */
-double ss_gradient(int param, Reflection *refl, UnitCell *cell,
-                   struct detgeom_panel *p)
-{
-	signed int h, k, l;
-	double yl, zl, kpred;
-	double asx, asy, asz, bsx, bsy, bsz, csx, csy, csz;
-
-	get_indices(refl, &h, &k, &l);
-	kpred = get_kpred(refl);
-	cell_get_reciprocal(cell, &asx, &asy, &asz,
-	                          &bsx, &bsy, &bsz,
-	                          &csx, &csy, &csz);
-	yl = h*asy + k*bsy + l*csy;
-	zl = h*asz + k*bsz + l*csz;
-
-	switch ( param ) {
-
-		case GPARAM_ASX :
-		return 0.0;
-
-		case GPARAM_BSX :
-		return 0.0;
-
-		case GPARAM_CSX :
-		return 0.0;
-
-		case GPARAM_ASY :
-		return h * p->cnz * p->pixel_pitch / (kpred + zl);
-
-		case GPARAM_BSY :
-		return k * p->cnz * p->pixel_pitch / (kpred + zl);
-
-		case GPARAM_CSY :
-		return l * p->cnz * p->pixel_pitch / (kpred + zl);
-
-		case GPARAM_ASZ :
-		return -h * yl * p->cnz * p->pixel_pitch / (kpred*kpred + 2.0*kpred*zl + zl*zl);
-
-		case GPARAM_BSZ :
-		return -k * yl * p->cnz * p->pixel_pitch / (kpred*kpred + 2.0*kpred*zl + zl*zl);
-
-		case GPARAM_CSZ :
-		return -l * yl * p->cnz * p->pixel_pitch / (kpred*kpred + 2.0*kpred*zl + zl*zl);
-
-		case GPARAM_DETX :
-		return 0;
-
-		case GPARAM_DETY :
-		return -1;
-
-		case GPARAM_CLEN :
-		return -yl / (kpred+zl);
-
-		default :
-		ERROR("Positional gradient requested for parameter %i?\n", param);
-		abort();
-
-	}
-}
diff --git a/libcrystfel/src/geometry.h b/libcrystfel/src/geometry.h
index a7057f12..d05c8832 100644
--- a/libcrystfel/src/geometry.h
+++ b/libcrystfel/src/geometry.h
@@ -62,35 +62,6 @@ typedef enum {
 } PartialityModel;
 
 
-/** Enumeration of parameters which may want to be refined */
-enum gparam {
-	GPARAM_ASX,
-	GPARAM_ASY,
-	GPARAM_ASZ,
-	GPARAM_BSX,
-	GPARAM_BSY,
-	GPARAM_BSZ,
-	GPARAM_CSX,
-	GPARAM_CSY,
-	GPARAM_CSZ,
-	GPARAM_R,
-	GPARAM_DIV,
-	GPARAM_DETX,
-	GPARAM_DETY,
-	GPARAM_CLEN,
-	GPARAM_OSF,   /* Linear scale factor */
-	GPARAM_BFAC,  /* D-W scale factor */
-	GPARAM_ANG1,  /* Out of plane rotation angles of crystal */
-	GPARAM_ANG2,
-	GPARAM_WAVELENGTH,
-	GPARAM_ROTX,  /* Detector panel rotation about +x */
-	GPARAM_ROTY,  /* Detector panel rotation about +y */
-	GPARAM_ROTZ,  /* Detector panel rotation about +z */
-
-	GPARAM_EOL    /* End of list */
-};
-
-
 /**
  * This structure represents the polarisation of the incident radiation
  */
@@ -110,8 +81,6 @@ extern RefList *predict_to_res(Crystal *cryst, double max_res);
 
 extern void calculate_partialities(Crystal *cryst, PartialityModel pmodel);
 
-extern double r_gradient(UnitCell *cell, int k, Reflection *refl,
-                         struct image *image);
 extern void update_predictions(Crystal *cryst);
 extern struct polarisation parse_polarisation(const char *text);
 extern void polarisation_correction(RefList *list, UnitCell *cell,
@@ -120,11 +89,6 @@ extern void polarisation_correction(RefList *list, UnitCell *cell,
 extern double sphere_fraction(double rlow, double rhigh, double pr);
 extern double gaussian_fraction(double rlow, double rhigh, double pr);
 
-extern double fs_gradient(int param, Reflection *refl, UnitCell *cell,
-                          struct detgeom_panel *p);
-extern double ss_gradient(int param, Reflection *refl, UnitCell *cell,
-                          struct detgeom_panel *p);
-
 #ifdef __cplusplus
 }
 #endif
diff --git a/libcrystfel/src/predict-refine.c b/libcrystfel/src/predict-refine.c
index ad7f5584..27f97700 100644
--- a/libcrystfel/src/predict-refine.c
+++ b/libcrystfel/src/predict-refine.c
@@ -68,7 +68,7 @@ static const enum gparam rv[] =
 };
 
 
-static double r_dev(struct reflpeak *rp)
+double r_dev(struct reflpeak *rp)
 {
 	/* Excitation error term */
 	return get_exerr(rp->refl);
@@ -91,6 +91,198 @@ double ss_dev(struct reflpeak *rp, struct detgeom *det)
 }
 
 
+double r_gradient(UnitCell *cell, int k, Reflection *refl, struct image *image)
+{
+	double asx, asy, asz;
+	double bsx, bsy, bsz;
+	double csx, csy, csz;
+	double xl, yl, zl;
+	signed int hs, ks, ls;
+	double tl, phi, azi;
+
+	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;
+
+	tl = sqrt(xl*xl + yl*yl);
+	phi = angle_between_2d(tl, zl+1.0/image->lambda, 0.0, 1.0); /* 2theta */
+	azi = atan2(yl, xl); /* azimuth */
+
+	switch ( k ) {
+
+		case GPARAM_ASX :
+		return - hs * sin(phi) * cos(azi);
+
+		case GPARAM_BSX :
+		return - ks * sin(phi) * cos(azi);
+
+		case GPARAM_CSX :
+		return - ls * sin(phi) * cos(azi);
+
+		case GPARAM_ASY :
+		return - hs * sin(phi) * sin(azi);
+
+		case GPARAM_BSY :
+		return - ks * sin(phi) * sin(azi);
+
+		case GPARAM_CSY :
+		return - ls * sin(phi) * sin(azi);
+
+		case GPARAM_ASZ :
+		return - hs * cos(phi);
+
+		case GPARAM_BSZ :
+		return - ks * cos(phi);
+
+		case GPARAM_CSZ :
+		return - ls * cos(phi);
+
+		case GPARAM_DETX :
+		case GPARAM_DETY :
+		case GPARAM_CLEN :
+		return 0.0;
+
+	}
+
+	ERROR("No r gradient defined for parameter %i\n", k);
+	abort();
+}
+
+
+/* Returns dfs_refl/dP, where P = any parameter *
+ * fs_refl is the fast scan position of refl in panel 'p',
+ * in metres (not pixels) */
+double fs_gradient(int param, Reflection *refl, UnitCell *cell,
+                   struct detgeom_panel *p)
+{
+	signed int h, k, l;
+	double xl, zl, kpred;
+	double asx, asy, asz, bsx, bsy, bsz, csx, csy, csz;
+
+	get_indices(refl, &h, &k, &l);
+	kpred = get_kpred(refl);
+	cell_get_reciprocal(cell, &asx, &asy, &asz,
+	                          &bsx, &bsy, &bsz,
+	                          &csx, &csy, &csz);
+	xl = h*asx + k*bsx + l*csx;
+	zl = h*asz + k*bsz + l*csz;
+
+	switch ( param ) {
+
+		case GPARAM_ASX :
+		return h * p->cnz * p->pixel_pitch / (kpred + zl);
+
+		case GPARAM_BSX :
+		return k * p->cnz * p->pixel_pitch / (kpred + zl);
+
+		case GPARAM_CSX :
+		return l * p->cnz * p->pixel_pitch / (kpred + zl);
+
+		case GPARAM_ASY :
+		return 0.0;
+
+		case GPARAM_BSY :
+		return 0.0;
+
+		case GPARAM_CSY :
+		return 0.0;
+
+		case GPARAM_ASZ :
+		return -h * xl * p->cnz * p->pixel_pitch / (kpred*kpred + 2.0*kpred*zl + zl*zl);
+
+		case GPARAM_BSZ :
+		return -k * xl * p->cnz * p->pixel_pitch / (kpred*kpred + 2.0*kpred*zl + zl*zl);
+
+		case GPARAM_CSZ :
+		return -l * xl * p->cnz * p->pixel_pitch / (kpred*kpred + 2.0*kpred*zl + zl*zl);
+
+		case GPARAM_DETX :
+		return -1;
+
+		case GPARAM_DETY :
+		return 0;
+
+		case GPARAM_CLEN :
+		return -xl / (kpred+zl);
+
+		default:
+		ERROR("Positional gradient requested for parameter %i?\n", param);
+		abort();
+
+	}
+
+}
+
+
+/* Returns dss_refl/dP, where P = any parameter
+ * ss_refl is the slow scan position of refl in panel 'p',
+ * in metres (not pixels) */
+double ss_gradient(int param, Reflection *refl, UnitCell *cell,
+                   struct detgeom_panel *p)
+{
+	signed int h, k, l;
+	double yl, zl, kpred;
+	double asx, asy, asz, bsx, bsy, bsz, csx, csy, csz;
+
+	get_indices(refl, &h, &k, &l);
+	kpred = get_kpred(refl);
+	cell_get_reciprocal(cell, &asx, &asy, &asz,
+	                          &bsx, &bsy, &bsz,
+	                          &csx, &csy, &csz);
+	yl = h*asy + k*bsy + l*csy;
+	zl = h*asz + k*bsz + l*csz;
+
+	switch ( param ) {
+
+		case GPARAM_ASX :
+		return 0.0;
+
+		case GPARAM_BSX :
+		return 0.0;
+
+		case GPARAM_CSX :
+		return 0.0;
+
+		case GPARAM_ASY :
+		return h * p->cnz * p->pixel_pitch / (kpred + zl);
+
+		case GPARAM_BSY :
+		return k * p->cnz * p->pixel_pitch / (kpred + zl);
+
+		case GPARAM_CSY :
+		return l * p->cnz * p->pixel_pitch / (kpred + zl);
+
+		case GPARAM_ASZ :
+		return -h * yl * p->cnz * p->pixel_pitch / (kpred*kpred + 2.0*kpred*zl + zl*zl);
+
+		case GPARAM_BSZ :
+		return -k * yl * p->cnz * p->pixel_pitch / (kpred*kpred + 2.0*kpred*zl + zl*zl);
+
+		case GPARAM_CSZ :
+		return -l * yl * p->cnz * p->pixel_pitch / (kpred*kpred + 2.0*kpred*zl + zl*zl);
+
+		case GPARAM_DETX :
+		return 0;
+
+		case GPARAM_DETY :
+		return -1;
+
+		case GPARAM_CLEN :
+		return -yl / (kpred+zl);
+
+		default :
+		ERROR("Positional gradient requested for parameter %i?\n", param);
+		abort();
+
+	}
+}
+
+
 static int cmpd2(const void *av, const void *bv)
 {
 	struct reflpeak *a, *b;
diff --git a/libcrystfel/src/predict-refine.h b/libcrystfel/src/predict-refine.h
index e74e00a8..c667dc6d 100644
--- a/libcrystfel/src/predict-refine.h
+++ b/libcrystfel/src/predict-refine.h
@@ -30,6 +30,36 @@
 #define PREDICT_REFINE_H
 
 struct reflpeak;
+
+/** Enumeration of parameters which may want to be refined */
+enum gparam {
+	GPARAM_ASX,
+	GPARAM_ASY,
+	GPARAM_ASZ,
+	GPARAM_BSX,
+	GPARAM_BSY,
+	GPARAM_BSZ,
+	GPARAM_CSX,
+	GPARAM_CSY,
+	GPARAM_CSZ,
+	GPARAM_R,
+	GPARAM_DIV,
+	GPARAM_DETX,
+	GPARAM_DETY,
+	GPARAM_CLEN,
+	GPARAM_OSF,   /* Linear scale factor */
+	GPARAM_BFAC,  /* D-W scale factor */
+	GPARAM_ANG1,  /* Out of plane rotation angles of crystal */
+	GPARAM_ANG2,
+	GPARAM_WAVELENGTH,
+	GPARAM_ROTX,  /* Detector panel rotation about +x */
+	GPARAM_ROTY,  /* Detector panel rotation about +y */
+	GPARAM_ROTZ,  /* Detector panel rotation about +z */
+
+	GPARAM_EOL    /* End of list */
+};
+
+
 #include "crystal.h"
 #include "crystfel-mille.h"
 
@@ -50,8 +80,19 @@ extern int refine_prediction(struct image *image, Crystal *cr, Mille *mille);
 
 extern int refine_radius(Crystal *cr, struct image *image);
 
+extern double r_dev(struct reflpeak *rp);
+
 extern double fs_dev(struct reflpeak *rp, struct detgeom *det);
 
 extern double ss_dev(struct reflpeak *rp, struct detgeom *det);
 
+extern double r_gradient(UnitCell *cell, int k, Reflection *refl,
+                         struct image *image);
+
+extern double fs_gradient(int param, Reflection *refl, UnitCell *cell,
+                          struct detgeom_panel *p);
+
+extern double ss_gradient(int param, Reflection *refl, UnitCell *cell,
+                          struct detgeom_panel *p);
+
 #endif	/* PREDICT_REFINE_H */