Cell vector gradients for SCSphere, plus general rationalisation

4 files changed, 81 insertions, 88 deletions

diff --git a/libcrystfel/src/geometry.c b/libcrystfel/src/geometry.c
index 250f889c..e08d9120 100644
--- a/libcrystfel/src/geometry.c
+++ b/libcrystfel/src/geometry.c
@@ -102,17 +102,64 @@ static signed int locate_peak(double x, double y, double z, double k,
 }
 
 
-static double partiality(PartialityModel pmodel,
-                         double rlow, double rmid, double rhigh,
-                         double r, double D)
+double sphere_fraction(double rlow, double rhigh, double pr)
 {
 	double qlow, qhigh;
 	double plow, phigh;
-	const double ng = 2.6;
+
+	/* If the "lower" Ewald sphere is a long way away, use the
+	 * position at which the Ewald sphere would just touch the
+	 * reflection.
+	 *
+	 * The six possible combinations of clamp_{low,high} (including
+	 * zero) correspond to the six situations in Table 3 of Rossmann
+	 * et al. (1979).
+	 */
+	if ( rlow < -pr ) rlow = -pr;
+	if ( rlow > +pr ) rlow = +pr;
+	if ( rhigh < -pr ) rhigh = -pr;
+	if ( rhigh > +pr ) rhigh = +pr;
 
 	/* Calculate degrees of penetration */
-	qlow  = (rlow + r)/(2.0*r);
-	qhigh = (rhigh + r)/(2.0*r);
+	qlow  = (rlow + pr)/(2.0*pr);
+	qhigh = (rhigh + pr)/(2.0*pr);
+
+	plow  = 3.0*qlow*qlow - 2.0*qlow*qlow*qlow;
+	phigh = 3.0*qhigh*qhigh - 2.0*qhigh*qhigh*qhigh;
+
+	return plow - phigh;
+}
+
+
+double gaussian_fraction(double rlow, double rhigh, double pr)
+{
+	double plow, phigh;
+	const double ng = 2.6;
+
+	/* If the "lower" Ewald sphere is a long way away, use the
+	 * position at which the Ewald sphere would just touch the
+	 * reflection.
+	 *
+	 * The six possible combinations of clamp_{low,high} (including
+	 * zero) correspond to the six situations in Table 3 of Rossmann
+	 * et al. (1979).
+	 */
+	if ( rlow < -pr ) rlow = -pr;
+	if ( rlow > +pr ) rlow = +pr;
+	if ( rhigh < -pr ) rhigh = -pr;
+	if ( rhigh > +pr ) rhigh = +pr;
+
+	plow = 0.5 * gsl_sf_erf(ng * rlow / (sqrt(2.0)*pr));
+	phigh = 0.5 * gsl_sf_erf(ng * rhigh / (sqrt(2.0)*pr));
+
+	return plow - phigh;
+}
+
+
+static double partiality(PartialityModel pmodel,
+                         double rlow, double rhigh, double pr)
+{
+	double D = rlow - rhigh;
 
 	/* Convert to partiality */
 	switch ( pmodel ) {
@@ -121,15 +168,11 @@ static double partiality(PartialityModel pmodel,
 		case PMODEL_UNITY:
 		return 1.0;
 
-		case PMODEL_SCGAUSSIAN:
-		plow = 0.5 * gsl_sf_erf(ng * rlow / (sqrt(2.0)*r));
-		phigh = 0.5 * gsl_sf_erf(ng * rhigh / (sqrt(2.0)*r));
-		return 4.0*(plow-phigh)*r / (3.0*D);
-
 		case PMODEL_SCSPHERE:
-		plow  = 3.0*qlow*qlow  - 2.0*qlow*qlow*qlow;
-		phigh = 3.0*qhigh*qhigh - 2.0*qhigh*qhigh*qhigh;
-		return 4.0*(plow-phigh)*r / (3.0*D);
+		return 4.0*sphere_fraction(rlow, rhigh, pr)*pr / (3.0*D);
+
+		case PMODEL_SCGAUSSIAN:
+		return 4.0*gaussian_fraction(rlow, rhigh, pr)*pr / (3.0*D);
 
 	}
 }
@@ -142,16 +185,13 @@ static Reflection *check_reflection(struct image *image, Crystal *cryst,
 {
 	const int output = 0;
 	double tl;
-	double rlow, rmid, rhigh;     /* "Excitation error" */
+	double rlow, rhigh;     /* "Excitation error" */
 	double part;            /* Partiality */
-	int clamp_low, clamp_high;
-	double klow, kmid, khigh;    /* Wavenumber */
+	double klow, khigh;    /* Wavenumber */
 	Reflection *refl;
 	double cet, cez;  /* Centre of Ewald sphere */
 	double pr;
-	double D;
 	double del;
-	double rlowuc, rhighuc;  /* Values before clamping */
 
 	/* Don't predict 000 */
 	if ( abs(h)+abs(k)+abs(l) == 0 ) return NULL;
@@ -163,7 +203,6 @@ static Reflection *check_reflection(struct image *image, Crystal *cryst,
 	 * "high" gives the smallest Ewald sphere (wavelength long => k small)
 	 */
 	klow = 1.0/(image->lambda - image->lambda*image->bw/2.0);
-	kmid = 1.0/image->lambda;
 	khigh = 1.0/(image->lambda + image->lambda*image->bw/2.0);
 
 	/* If the point is looking "backscattery", reject it straight away */
@@ -175,10 +214,6 @@ static Reflection *check_reflection(struct image *image, Crystal *cryst,
 	cez = -cos(del/2.0) * khigh;
 	rhigh = khigh - distance(cet, cez, tl, zl);  /* Loss of precision */
 
-	cet = 0.0;
-	cez = -kmid;
-	rmid = kmid - distance(cet, cez, tl, zl);  /* Loss of precision */
-
 	cet =  sin(del/2.0) * klow;
 	cez = -cos(del/2.0) * klow;
 	rlow = klow - distance(cet, cez, tl, zl);  /* Loss of precision */
@@ -196,38 +231,8 @@ static Reflection *check_reflection(struct image *image, Crystal *cryst,
 	     && (fabs(rlow) > pr)
 	     && (fabs(rhigh) > pr) ) return NULL;
 
-	D = rlow - rhigh;
-
-	rlowuc = rlow;  rhighuc = rhigh;
-
-	/* If the "lower" Ewald sphere is a long way away, use the
-	 * position at which the Ewald sphere would just touch the
-	 * reflection.
-	 *
-	 * The six possible combinations of clamp_{low,high} (including
-	 * zero) correspond to the six situations in Table 3 of Rossmann
-	 * et al. (1979).
-	 */
-	clamp_low = 0;  clamp_high = 0;
-	if ( rlow < -pr ) {
-		rlow = -pr;
-		clamp_low = -1;
-	}
-	if ( rlow > +pr ) {
-		rlow = +pr;
-		clamp_low = +1;
-	}
-	if ( rhigh < -pr ) {
-		rhigh = -pr;
-		clamp_high = -1;
-	}
-	if ( rhigh > +pr ) {
-		rhigh = +pr;
-		clamp_high = +1;
-	}
-
 	/* Calculate partiality */
-	part = partiality(pmodel, rlow, rmid, rhigh, pr, D);
+	part = partiality(pmodel, rlow, rhigh, pr);
 
 	/* Add peak to list */
 	refl = reflection_new(h, k, l);
@@ -246,7 +251,7 @@ static Reflection *check_reflection(struct image *image, Crystal *cryst,
 		set_detector_pos(refl, 0.0, xda, yda);
 	}
 
-	set_partial(refl, rlowuc, rhighuc, part, clamp_low, clamp_high);
+	set_partial(refl, rlow, rhigh, part);
 	set_lorentz(refl, 1.0);
 	set_symmetric_indices(refl, h, k, l);
 	set_redundancy(refl, 1);
@@ -442,7 +447,6 @@ void update_partialities_2(Crystal *cryst, PartialityModel pmodel,
 		double r1, r2, L, p, x, y;
 		double xl, yl, zl;
 		signed int h, k, l;
-		int clamp1, clamp2;
 		double old_p;
 
 		get_symmetric_indices(refl, &h, &k, &l);
@@ -472,8 +476,8 @@ void update_partialities_2(Crystal *cryst, PartialityModel pmodel,
 			}
 
 			/* Transfer partiality stuff */
-			get_partial(vals, &r1, &r2, &p, &clamp1, &clamp2);
-			set_partial(refl, r1, r2, p, clamp1, clamp2);
+			get_partial(vals, &r1, &r2, &p);
+			set_partial(refl, r1, r2, p);
 			L = get_lorentz(vals);
 			set_lorentz(refl, L);
 
diff --git a/libcrystfel/src/geometry.h b/libcrystfel/src/geometry.h
index ccdea1c2..93aaf619 100644
--- a/libcrystfel/src/geometry.h
+++ b/libcrystfel/src/geometry.h
@@ -74,7 +74,8 @@ extern void update_partialities_2(Crystal *cryst, PartialityModel pmodel,
 extern void polarisation_correction(RefList *list, UnitCell *cell,
                                     struct image *image);
 
-#define LORENTZ_SCALE (0.01e9)
+extern double sphere_fraction(double rlow, double rhigh, double pr);
+extern double gaussian_fraction(double rlow, double rhigh, double pr);
 
 #ifdef __cplusplus
 }
diff --git a/libcrystfel/src/reflist.c b/libcrystfel/src/reflist.c
index afc76caa..33c1f948 100644
--- a/libcrystfel/src/reflist.c
+++ b/libcrystfel/src/reflist.c
@@ -67,12 +67,10 @@ struct _refldata {
 	signed int ls;
 
 	/* Partiality and related geometrical stuff */
-	double r1;  /* First excitation error */
-	double r2;  /* Second excitation error */
+	double rlow;  /* Low excitation error */
+	double rhigh;  /* High excitation error */
 	double p;   /* Partiality */
 	double L;   /* Lorentz factor */
-	int clamp1; /* Clamp status for r1 */
-	int clamp2; /* Clamp status for r2 */
 
 	/* Location in image */
 	double fs;
@@ -424,24 +422,20 @@ double get_intensity(const Reflection *refl)
 /**
  * get_partial:
  * @refl: A %Reflection
- * @r1: Location at which to store the first excitation error
- * @r2: Location at which to store the second excitation error
+ * @rlow: Location at which to store the "low" excitation error
+ * @rhigh: Location at which to store the "high" excitation error
  * @p: Location at which to store the partiality
- * @clamp_low: Location at which to store the first clamp status
- * @clamp_high: Location at which to store the second clamp status
  *
  * This function is used during post refinement (in conjunction with
  * set_partial()) to get access to the details of the partiality calculation.
  *
  **/
-void get_partial(const Reflection *refl, double *r1, double *r2, double *p,
-                 int *clamp_low, int *clamp_high)
+void get_partial(const Reflection *refl, double *rlow, double *rhigh,
+                 double *p)
 {
-	*r1 = refl->data.r1;
-	*r2 = refl->data.r2;
+	*rlow = refl->data.rlow;
+	*rhigh = refl->data.rhigh;
 	*p = get_partiality(refl);
-	*clamp_low = refl->data.clamp1;
-	*clamp_high = refl->data.clamp2;
 }
 
 
@@ -592,24 +586,19 @@ void set_detector_pos(Reflection *refl, double exerr, double fs, double ss)
 /**
  * set_partial:
  * @refl: A %Reflection
- * @r1: The first excitation error
- * @r2: The second excitation error
+ * @rlow: The "low" excitation error
+ * @rhigh: The "high" excitation error
  * @p: The partiality
- * @clamp_low: The first clamp status
- * @clamp_high: The second clamp status
  *
  * This function is used during post refinement (in conjunction with
  * get_partial()) to get access to the details of the partiality calculation.
  *
  **/
-void set_partial(Reflection *refl, double r1, double r2, double p,
-                 double clamp_low, double clamp_high)
+void set_partial(Reflection *refl, double rlow, double rhigh, double p)
 {
-	refl->data.r1 = r1;
-	refl->data.r2 = r2;
+	refl->data.rlow = rlow;
+	refl->data.rhigh = rhigh;
 	refl->data.p = p;
-	refl->data.clamp1 = clamp_low;
-	refl->data.clamp2 = clamp_high;
 }
 
 
diff --git a/libcrystfel/src/reflist.h b/libcrystfel/src/reflist.h
index 6892d7f2..85a87c54 100644
--- a/libcrystfel/src/reflist.h
+++ b/libcrystfel/src/reflist.h
@@ -94,8 +94,8 @@ extern void get_symmetric_indices(const Reflection *refl,
                                   signed int *hs, signed int *ks,
                                   signed int *ls);
 extern double get_intensity(const Reflection *refl);
-extern void get_partial(const Reflection *refl, double *r1, double *r2,
-                        double *p, int *clamp_low, int *clamp_high);
+extern void get_partial(const Reflection *refl, double *rlow, double *rhigh,
+                        double *p);
 extern int get_redundancy(const Reflection *refl);
 extern double get_temp1(const Reflection *refl);
 extern double get_temp2(const Reflection *refl);
@@ -108,8 +108,7 @@ extern double get_mean_bg(const Reflection *refl);
 extern void copy_data(Reflection *to, const Reflection *from);
 extern void set_detector_pos(Reflection *refl, double exerr,
                              double fs, double ss);
-extern void set_partial(Reflection *refl, double r1, double r2, double p,
-                        double clamp_low, double clamp_high);
+extern void set_partial(Reflection *refl, double rlow, double rhigh, double p);
 extern void set_partiality(Reflection *refl, double p);
 extern void set_lorentz(Reflection *refl, double L);
 extern void set_intensity(Reflection *refl, double intensity);