Add Thin Ewald Sphere model

7 files changed, 155 insertions, 80 deletions

diff --git a/libcrystfel/src/geometry.c b/libcrystfel/src/geometry.c
index 5cc46b2e..3586ed0b 100644
--- a/libcrystfel/src/geometry.c
+++ b/libcrystfel/src/geometry.c
@@ -103,7 +103,8 @@ static signed int locate_peak(double x, double y, double z, double k,
 }
 
 
-static double partiality(PartialityModel pmodel, double rlow, double rhigh,
+static double partiality(PartialityModel pmodel,
+                         double rlow, double rmid, double rhigh,
                          double r)
 {
 	double qlow, qhigh;
@@ -119,23 +120,22 @@ static double partiality(PartialityModel pmodel, double rlow, double rhigh,
 
 		default:
 		case PMODEL_UNITY:
-		plow = 1.0;
-		phigh = 0.0;
-		break;
+		return 1.0;
 
 		case PMODEL_SPHERE:
-		plow  = 3.0*pow(qlow,2.0)  - 2.0*pow(qlow,3.0);
-		phigh = 3.0*pow(qhigh,2.0) - 2.0*pow(qhigh,3.0);
-		break;
+		plow  = 3.0*qlow*qlow  - 2.0*qlow*qlow*qlow;
+		phigh = 3.0*qhigh*qhigh - 2.0*qhigh*qhigh*qhigh;
+		return plow - phigh;
 
 		case PMODEL_GAUSSIAN:
 		plow = 0.5 * gsl_sf_erf(ng * rlow / (sqrt(2.0)*r));
 		phigh = 0.5 * gsl_sf_erf(ng * rhigh / (sqrt(2.0)*r));
-		break;
+		return plow - phigh;
 
-	}
+		case PMODEL_THIN:
+		return 1.0 - (rmid*rmid)/(r*r);
 
-	return plow - phigh;
+	}
 }
 
 
@@ -146,12 +146,12 @@ static Reflection *check_reflection(struct image *image, Crystal *cryst,
 {
 	const int output = 0;
 	double tl;
-	double rlow, rhigh;     /* "Excitation error" */
+	double rlow, rmid, rhigh;     /* "Excitation error" */
 	double part;            /* Partiality */
 	int clamp_low, clamp_high;
-	double klow, khigh;    /* Wavenumber */
+	double klow, kmid, khigh;    /* Wavenumber */
 	Reflection *refl;
-	double cet, cez;
+	double cet, cez;  /* Centre of Ewald sphere */
 	double pr;
 	double L;
 	double del;
@@ -166,6 +166,7 @@ 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 */
@@ -177,14 +178,14 @@ 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 */
 
-	if ( (signbit(rlow) == signbit(rhigh))
-	     && (fabs(rlow) > pr)
-	     && (fabs(rhigh) > pr) ) return NULL;
-
 	if ( unlikely(rlow < rhigh) ) {
 		ERROR("Reflection with rlow < rhigh!\n");
 		ERROR("%3i %3i %3i  rlow = %e, rhigh = %e\n",
@@ -193,10 +194,41 @@ static Reflection *check_reflection(struct image *image, Crystal *cryst,
 		return NULL;
 	}
 
+	/* Conditions for reflection to be excited at all */
+	switch ( pmodel ) {
+
+		default:
+		case PMODEL_UNITY:  /* PMODEL_UNITY shouldn't end up here */
+		case PMODEL_SPHERE:
+		case PMODEL_GAUSSIAN:
+		if ( (signbit(rlow) == signbit(rhigh))
+		     && (fabs(rlow) > pr)
+		     && (fabs(rhigh) > pr) ) return NULL;
+		break;
+
+		case PMODEL_THIN:
+		if ( fabs(rmid) > pr ) return NULL;
+		break;
+
+	}
+
 	/* Lorentz factor is determined direction from the r values, before
 	 * clamping.  The multiplication by 0.01e9 to make the
 	 * correction factor vaguely near 1. */
-	L = LORENTZ_SCALE / (rlow - rhigh);
+	switch ( pmodel ) {
+
+		default:
+		case PMODEL_SPHERE:
+		case PMODEL_GAUSSIAN:
+		L = LORENTZ_SCALE / (rlow - rhigh);
+		break;
+
+		case PMODEL_UNITY:  /* PMODEL_UNITY shouldn't end up here */
+		case PMODEL_THIN:
+		L = 1.0;
+		break;
+
+	}
 
 	/* If the "lower" Ewald sphere is a long way away, use the
 	 * position at which the Ewald sphere would just touch the
@@ -225,7 +257,7 @@ static Reflection *check_reflection(struct image *image, Crystal *cryst,
 	}
 
 	/* Calculate partiality */
-	part = partiality(pmodel, rlow, rhigh, pr);
+	part = partiality(pmodel, rlow, rmid, rhigh, pr);
 
 	/* Add peak to list */
 	refl = reflection_new(h, k, l);
@@ -244,7 +276,18 @@ static Reflection *check_reflection(struct image *image, Crystal *cryst,
 		set_detector_pos(refl, 0.0, xda, yda);
 	}
 
-	set_partial(refl, rlow, rhigh, part, clamp_low, clamp_high);
+	if ( pmodel != PMODEL_THIN ) {
+		set_partial(refl, rlow, rhigh, part, clamp_low, clamp_high);
+	} else {
+		/* If we are using the TES (Thin Ewald Sphere) model, we abuse
+		 * the fields as follows:
+		 *  rlow       = the r value for the middle (only) Ewald sphere
+		 *  rhigh      = 0.0
+		 *  clamp_low  =  0
+		 *  clamp_high = +1
+		 */
+		set_partial(refl, rmid, 0.0, part, 0, +1);
+	}
 	set_lorentz(refl, L);
 	set_symmetric_indices(refl, h, k, l);
 	set_redundancy(refl, 1);
diff --git a/libcrystfel/src/geometry.h b/libcrystfel/src/geometry.h
index 1f465167..8041936a 100644
--- a/libcrystfel/src/geometry.h
+++ b/libcrystfel/src/geometry.h
@@ -50,6 +50,7 @@ extern "C" {
  *   space.
  * @PMODEL_UNITY : Set all all partialities and Lorentz factors to 1.
  * @PMODEL_GAUSSIAN : Gaussian profiles in 3D
+ * @PMODEL_THIN : Thin Ewald sphere intersecting sphere
  *
  * A %PartialityModel describes a geometrical model which can be used to
  * calculate spot partialities and Lorentz correction factors.
@@ -59,6 +60,7 @@ typedef enum {
 	PMODEL_SPHERE,
 	PMODEL_UNITY,
 	PMODEL_GAUSSIAN,
+	PMODEL_THIN,
 
 } PartialityModel;
 
diff --git a/src/partialator.c b/src/partialator.c
index 2ef01dba..e5083611 100644
--- a/src/partialator.c
+++ b/src/partialator.c
@@ -348,6 +348,8 @@ int main(int argc, char *argv[])
 			pmodel = PMODEL_UNITY;
 		} else if ( strcmp(pmodel_str, "gaussian") == 0 ) {
 			pmodel = PMODEL_GAUSSIAN;
+		} else if ( strcmp(pmodel_str, "thin") == 0 ) {
+			pmodel = PMODEL_THIN;
 		} else {
 			ERROR("Unknown partiality model '%s'.\n", pmodel_str);
 			return 1;
diff --git a/src/post-refinement.c b/src/post-refinement.c
index 5de8a246..99ddf6b1 100644
--- a/src/post-refinement.c
+++ b/src/post-refinement.c
@@ -55,15 +55,23 @@
 #define MAX_CYCLES (10)
 
 
-static double dpdq(double r, double profile_radius, PartialityModel pmodel)
+static double dpdq(double r, double profile_radius)
 {
 	double q;
-	double ng = 3.0;
 
 	/* Calculate degree of penetration */
 	q = (r + profile_radius)/(2.0*profile_radius);
 
-	/* dp/dq */
+	return 6.0*(q-q*q);
+}
+
+
+/* Returns dp/dr at "r" */
+static double partiality_gradient(double r, double profile_radius,
+                                  PartialityModel pmodel)
+{
+	double dqdr;  /* dq/dr */
+
 	switch ( pmodel ) {
 
 		default:
@@ -71,39 +79,46 @@ static double dpdq(double r, double profile_radius, PartialityModel pmodel)
 		return 0.0;
 
 		case PMODEL_SPHERE:
-		return 6.0*(q-pow(q, 2.0));
+		dqdr = 1.0 / (2.0*profile_radius);
+		return dpdq(r, profile_radius) * dqdr;
 
 		case PMODEL_GAUSSIAN:
-		/* The flat sphere model is close enough */
-		return 6.0*(q-pow(q, 2.0));
+		/* FIXME: Get a proper gradient */
+		dqdr = 1.0 / (2.0*profile_radius);
+		return dpdq(r, profile_radius) * dqdr;
+
+		case PMODEL_THIN:
+		return -(2.0*r)/(profile_radius*profile_radius);
 
 	}
 }
 
 
-/* Returns dp/dr at "r" */
-static double partiality_gradient(double r, double profile_radius,
-                                  PartialityModel pmodel)
+/* Returns dp/drad at "r" */
+static double partiality_rgradient(double r, double profile_radius,
+                                   PartialityModel pmodel)
 {
-	double dqdr;
+	double dqdrad;  /* dq/drad */
 
-	/* dq/dr */
-	dqdr = 1.0 / (2.0*profile_radius);
+	switch ( pmodel ) {
 
-	return dpdq(r, profile_radius, pmodel) * dqdr;
-}
+		default:
+		case PMODEL_UNITY:
+		return 0.0;
 
+		case PMODEL_SPHERE:
+		dqdrad = -0.5 * r / (profile_radius * profile_radius);
+		return dpdq(r, profile_radius) * dqdrad;
 
-/* Returns dp/drad at "r" */
-static double partiality_rgradient(double r, double profile_radius,
-                                   PartialityModel pmodel)
-{
-	double dqdrad;
+		case PMODEL_GAUSSIAN:
+		/* FIXME: Get a proper gradient */
+		dqdrad = -0.5 * r / (profile_radius * profile_radius);
+		return dpdq(r, profile_radius) * dqdrad;
 
-	/* dq/drad */
-	dqdrad = -0.5 * r * pow(profile_radius, -2.0);
+		case PMODEL_THIN:
+		return 2.0*r*r*pow(profile_radius, -3.0);
 
-	return dpdq(r, profile_radius, pmodel) * dqdrad;
+	}
 }
 
 
@@ -173,19 +188,48 @@ double p_gradient(Crystal *cr, int k, Reflection *refl, PartialityModel pmodel)
 		ghigh = 0.0;
 	}
 
+	if ( k == REF_R ) {
+		switch ( pmodel ) {
+
+			default:
+			case PMODEL_UNITY:
+			return 0.0;
+
+			case PMODEL_GAUSSIAN:
+			gr  = partiality_rgradient(rlow, r, pmodel);
+			gr -= partiality_rgradient(rhigh, r, pmodel);
+			return gr;
+
+			case PMODEL_SPHERE:
+			gr  = partiality_rgradient(rlow, r, pmodel);
+			gr -= partiality_rgradient(rhigh, r, pmodel);
+			return gr;
+
+			case PMODEL_THIN:
+			return 2.0*rlow*rlow/(r*r*r);
+		}
+	}
+
+	if ( k == REF_DIV ) {
+		switch ( pmodel ) {
+
+			default:
+			case PMODEL_UNITY:
+			return 0.0;
+
+			case PMODEL_GAUSSIAN:
+			case PMODEL_SPHERE:
+			return (ds*glow + ds*ghigh) / 2.0;
+
+			case PMODEL_THIN:
+			return 0.0;
+		}
+	}
+
 	/* For many gradients, just multiply the above number by the gradient
 	 * of excitation error wrt whatever. */
 	switch ( k ) {
 
-		case REF_DIV :
-		/* Small angle approximation */
-		return (ds*glow + ds*ghigh) / 2.0;
-
-		case REF_R :
-		gr  = partiality_rgradient(rlow, r, pmodel);
-		gr -= partiality_rgradient(rhigh, r, pmodel);
-		return gr;
-
 		/* Cell parameters and orientation */
 		case REF_ASX :
 		return hs * sin(phi) * cos(azi) * (ghigh-glow);
@@ -227,37 +271,16 @@ double l_gradient(Crystal *cr, int k, Reflection *refl, PartialityModel pmodel)
 {
 	double ds;
 	signed int hs, ks, ls;
+	double L;
 
-	switch ( k ) {
-
-		/* Cell parameters do not affect Lorentz factor */
-		case REF_ASX :
-		case REF_BSX :
-		case REF_CSX :
-		case REF_ASY :
-		case REF_BSY :
-		case REF_CSY :
-		case REF_ASZ :
-		case REF_BSZ :
-		case REF_CSZ :
-		return 0.0;
-
-		/* Nor does change of radius */
-		case REF_R :
-		return 0.0;
-
-		default:
-		break;
-
-	}
-
-	assert(k == REF_DIV);
+	if ( k != REF_DIV ) return 0.0;
 
 	get_symmetric_indices(refl, &hs, &ks, &ls);
 
 	ds = 2.0 * resolution(crystal_get_cell(cr), hs, ks, ls);
 
-	return -ds*pow(get_lorentz(refl), 2.0) / LORENTZ_SCALE;
+	L = get_lorentz(refl);
+	return -ds*L*L / LORENTZ_SCALE;
 }
 
 
diff --git a/tests/pr_l_gradient_check.c b/tests/pr_l_gradient_check.c
index b707629f..7342d5f8 100644
--- a/tests/pr_l_gradient_check.c
+++ b/tests/pr_l_gradient_check.c
@@ -350,6 +350,7 @@ int main(int argc, char *argv[])
 			pmodel = PMODEL_GAUSSIAN;
 			STATUS("Testing Gaussian model:\n");
 		}
+		/* No point testing TES model, because it has no Lorentz factor */
 
 		orientation = random_quaternion(rng);
 		rot = cell_rotate(cell, orientation);
diff --git a/tests/pr_p_gradient_check.c b/tests/pr_p_gradient_check.c
index 8f7a66e4..14a2b579 100644
--- a/tests/pr_p_gradient_check.c
+++ b/tests/pr_p_gradient_check.c
@@ -174,6 +174,7 @@ static Crystal *new_shifted_crystal(Crystal *cr, int refine, double incr_val)
 	return cr_new;
 }
 
+
 static void calc_either_side(Crystal *cr, double incr_val,
                              int *valid, long double *vals[3], int refine,
                              PartialityModel pmodel)
@@ -225,7 +226,6 @@ static void calc_either_side(Crystal *cr, double incr_val,
 }
 
 
-
 static double test_gradients(Crystal *cr, double incr_val, int refine,
                              const char *str, const char *file,
                              PartialityModel pmodel, int quiet, int plot)
@@ -450,7 +450,7 @@ int main(int argc, char *argv[])
 
 	rng = gsl_rng_alloc(gsl_rng_mt19937);
 
-	for ( i=0; i<2; i++ ) {
+	for ( i=0; i<3; i++ ) {
 
 		UnitCell *rot;
 		double val;
@@ -459,9 +459,12 @@ int main(int argc, char *argv[])
 		if ( i == 0 ) {
 			pmodel = PMODEL_SPHERE;
 			STATUS("Testing flat sphere model:\n");
-		} else {
+		} else if ( i == 1 ) {
 			pmodel = PMODEL_GAUSSIAN;
 			STATUS("Testing Gaussian model:\n");
+		} else {
+			pmodel = PMODEL_THIN;
+			STATUS("Testing thin Ewald sphere model:\n");
 		}
 
 		orientation = random_quaternion(rng);
diff --git a/tests/pr_pl_gradient_check.c b/tests/pr_pl_gradient_check.c
index f83af335..69087a4b 100644
--- a/tests/pr_pl_gradient_check.c
+++ b/tests/pr_pl_gradient_check.c
@@ -174,6 +174,7 @@ static Crystal *new_shifted_crystal(Crystal *cr, int refine, double incr_val)
 	return cr_new;
 }
 
+
 static void calc_either_side(Crystal *cr, double incr_val,
                              int *valid, long double *vals[3], int refine,
                              PartialityModel pmodel)
@@ -225,7 +226,6 @@ static void calc_either_side(Crystal *cr, double incr_val,
 }
 
 
-
 static double test_gradients(Crystal *cr, double incr_val, int refine,
                              const char *str, const char *file,
                              PartialityModel pmodel, int quiet, int plot)
@@ -465,6 +465,7 @@ int main(int argc, char *argv[])
 			pmodel = PMODEL_GAUSSIAN;
 			STATUS("Testing Gaussian model:\n");
 		}
+		/* No point testing TES model */
 
 		orientation = random_quaternion(rng);
 		rot = cell_rotate(cell, orientation);