Add Gaussian partiality model

8 files changed, 107 insertions, 94 deletions

diff --git a/libcrystfel/src/geometry.c b/libcrystfel/src/geometry.c
index 0c0a09b8..357f64f8 100644
--- a/libcrystfel/src/geometry.c
+++ b/libcrystfel/src/geometry.c
@@ -3,11 +3,11 @@
  *
  * Geometry of diffraction
  *
- * Copyright © 2012-2013 Deutsches Elektronen-Synchrotron DESY,
+ * Copyright © 2012-2014 Deutsches Elektronen-Synchrotron DESY,
  *                       a research centre of the Helmholtz Association.
  *
  * Authors:
- *   2010-2013 Thomas White <taw@physics.org>
+ *   2010-2014 Thomas White <taw@physics.org>
  *
  * This file is part of CrystFEL.
  *
@@ -34,6 +34,7 @@
 #include <stdlib.h>
 #include <assert.h>
 #include <fenv.h>
+#include <gsl/gsl_sf_erf.h>
 
 #include "utils.h"
 #include "cell.h"
@@ -102,24 +103,44 @@ static signed int locate_peak(double x, double y, double z, double k,
 }
 
 
-static double partiality(double rlow, double rhigh, double r)
+static double partiality(PartialityModel pmodel, double rlow, double rhigh,
+                         double r)
 {
 	double qlow, qhigh;
 	double plow, phigh;
+	const double ng = 2.6;
 
 	/* Calculate degrees of penetration */
 	qlow  = (rlow + r)/(2.0*r);
 	qhigh = (rhigh + r)/(2.0*r);
 
 	/* Convert to partiality */
-	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);
+	switch ( pmodel ) {
+
+		default:
+		case PMODEL_UNITY:
+		plow = 1.0;
+		phigh = 0.0;
+		break;
+
+		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;
+
+		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;
 }
 
 
 static Reflection *check_reflection(struct image *image, Crystal *cryst,
+                                    PartialityModel pmodel,
                                     signed int h, signed int k, signed int l,
                                     double xl, double yl, double zl)
 {
@@ -204,7 +225,7 @@ static Reflection *check_reflection(struct image *image, Crystal *cryst,
 	}
 
 	/* Calculate partiality */
-	part = partiality(rlow, rhigh, pr);
+	part = partiality(pmodel, rlow, rhigh, pr);
 
 	/* Add peak to list */
 	refl = reflection_new(h, k, l);
@@ -298,7 +319,8 @@ RefList *find_intersections(struct image *image, Crystal *cryst)
 		yl = h*asy + k*bsy + l*csy;
 		zl = h*asz + k*bsz + l*csz;
 
-		refl = check_reflection(image, cryst, h, k, l, xl, yl, zl);
+		refl = check_reflection(image, cryst, PMODEL_SPHERE,
+		                        h, k, l, xl, yl, zl);
 
 		if ( refl != NULL ) {
 			add_refl_to_list(refl, reflections);
@@ -401,9 +423,6 @@ void update_partialities_2(Crystal *cryst, PartialityModel pmodel,
 	int n = 0;
 
 	if ( pmodel == PMODEL_UNITY ) {
-		/* It isn't strictly necessary to set the partialities to 1,
-		 * because the scaling stuff will just a correction factor of
-		 * 1 anyway.  This is just to help things make sense. */
 		set_unity_partialities(cryst);
 		return;
 	}
@@ -430,7 +449,8 @@ void update_partialities_2(Crystal *cryst, PartialityModel pmodel,
 		yl = h*asy + k*bsy + l*csy;
 		zl = h*asz + k*bsz + l*csz;
 
-		vals = check_reflection(image, cryst, h, k, l, xl, yl, zl);
+		vals = check_reflection(image, cryst, pmodel,
+		                        h, k, l, xl, yl, zl);
 
 		if ( vals == NULL ) {
 
diff --git a/libcrystfel/src/geometry.h b/libcrystfel/src/geometry.h
index fa7ed53b..7188f6ed 100644
--- a/libcrystfel/src/geometry.h
+++ b/libcrystfel/src/geometry.h
@@ -3,12 +3,12 @@
  *
  * Geometry of diffraction
  *
- * Copyright © 2013 Deutsches Elektronen-Synchrotron DESY,
- *                  a research centre of the Helmholtz Association.
+ * Copyright © 2013-2014 Deutsches Elektronen-Synchrotron DESY,
+ *                       a research centre of the Helmholtz Association.
  * Copyright © 2012 Richard Kirian
  *
  * Authors:
- *   2010-2013 Thomas White <taw@physics.org>
+ *   2010-2014 Thomas White <taw@physics.org>
  *   2012      Richard Kirian
  *
  * This file is part of CrystFEL.
@@ -49,6 +49,7 @@ extern "C" {
  * @PMODEL_SPHERE : Intersection of sphere with excited volume of reciprocal
  *   space.
  * @PMODEL_UNITY : Set all all partialities and Lorentz factors to 1.
+ * @PMODEL_GAUSSIAN : Gaussian profiles in 3D
  *
  * A %PartialityModel describes a geometrical model which can be used to
  * calculate spot partialities and Lorentz correction factors.
@@ -57,6 +58,7 @@ typedef enum {
 
 	PMODEL_SPHERE,
 	PMODEL_UNITY,
+	PMODEL_GAUSSIAN,
 
 } PartialityModel;
 
diff --git a/src/hrs-scaling.c b/src/hrs-scaling.c
index 76d1ff75..ac9091ed 100644
--- a/src/hrs-scaling.c
+++ b/src/hrs-scaling.c
@@ -119,24 +119,7 @@ static void run_scale_job(void *vwargs, int cookie)
 			Ih = get_intensity(r);
 		}
 
-		/* If you change this, be sure to also change
-		 * run_merge_job() and run_esd_job(). */
-		switch ( wargs->pmodel ) {
-
-			case PMODEL_UNITY :
-			corr = 1.0;
-			break;
-
-			case PMODEL_SPHERE :
-			corr = get_partiality(refl) * get_lorentz(refl);
-			break;
-
-			default :
-			ERROR("Unrecognised partiality model!\n");
-			abort();
-			break;
-
-		}
+		corr = get_partiality(refl) * get_lorentz(refl);
 
 		Ihl = get_intensity(refl) / corr;
 
@@ -285,24 +268,7 @@ static void run_merge_job(void *vwargs, int cookie)
 
 		}
 
-		/* If you change this, be sure to also change
-		 * run_scale_job() and run_esd_job(). */
-		switch ( wargs->pmodel ) {
-
-			case PMODEL_UNITY :
-			corr = 1.0;
-			break;
-
-			case PMODEL_SPHERE :
-			corr = get_partiality(refl) * get_lorentz(refl);
-			break;
-
-			default :
-			ERROR("Unrecognised partiality model!\n");
-			abort();
-			break;
-
-		}
+		corr = get_partiality(refl) * get_lorentz(refl);
 
 		Ihl = get_intensity(refl) / corr;
 
@@ -421,24 +387,7 @@ static void run_esd_job(void *vwargs, int cookie)
 
 		num = get_temp1(f);
 
-		/* If you change this, be sure to also change
-		 * run_scale_job() and run_merge_job(). */
-		switch ( wargs->pmodel ) {
-
-			case PMODEL_UNITY :
-			corr = 1.0;
-			break;
-
-			case PMODEL_SPHERE :
-			corr = get_partiality(refl) * get_lorentz(refl);
-			break;;
-
-			default :
-			ERROR("Unrecognised partiality model!\n");
-			abort();
-			break;
-
-		}
+		corr = get_partiality(refl) * get_lorentz(refl);
 
 		Ih = get_intensity(f);
 		Ihl = G * get_intensity(refl) / corr;
diff --git a/src/partialator.c b/src/partialator.c
index bc8704b0..c1dedad8 100644
--- a/src/partialator.c
+++ b/src/partialator.c
@@ -468,6 +468,8 @@ int main(int argc, char *argv[])
 			pmodel = PMODEL_SPHERE;
 		} else if ( strcmp(pmodel_str, "unity") == 0 ) {
 			pmodel = PMODEL_UNITY;
+		} else if ( strcmp(pmodel_str, "gaussian") == 0 ) {
+			pmodel = PMODEL_GAUSSIAN;
 		} else {
 			ERROR("Unknown partiality model '%s'.\n", pmodel_str);
 			return 1;
diff --git a/src/post-refinement.c b/src/post-refinement.c
index d5d8a4d1..8fd29783 100644
--- a/src/post-refinement.c
+++ b/src/post-refinement.c
@@ -3,11 +3,11 @@
  *
  * Post refinement
  *
- * Copyright © 2012-2013 Deutsches Elektronen-Synchrotron DESY,
+ * Copyright © 2012-2014 Deutsches Elektronen-Synchrotron DESY,
  *                       a research centre of the Helmholtz Association.
  *
  * Authors:
- *   2010-2013 Thomas White <taw@physics.org>
+ *   2010-2014 Thomas White <taw@physics.org>
  *
  * This file is part of CrystFEL.
  *
@@ -52,39 +52,55 @@
 #define MAX_CYCLES (10)
 
 
-/* Returns dp/dr at "r" */
-static double partiality_gradient(double r, double profile_radius)
+static double dpdq(double r, double profile_radius, PartialityModel pmodel)
 {
-	double q, dpdq, dqdr;
+	double q;
+	double ng = 3.0;
 
 	/* Calculate degree of penetration */
 	q = (r + profile_radius)/(2.0*profile_radius);
 
 	/* dp/dq */
-	dpdq = 6.0*(q-pow(q, 2.0));
+	switch ( pmodel ) {
+
+		default:
+		case PMODEL_UNITY:
+		return 0.0;
+
+		case PMODEL_SPHERE:
+		return 6.0*(q-pow(q, 2.0));
+
+		case PMODEL_GAUSSIAN:
+		/* The flat sphere model is close enough */
+		return 6.0*(q-pow(q, 2.0));
+
+	}
+}
+
+
+/* Returns dp/dr at "r" */
+static double partiality_gradient(double r, double profile_radius,
+                                  PartialityModel pmodel)
+{
+	double dqdr;
 
 	/* dq/dr */
 	dqdr = 1.0 / (2.0*profile_radius);
 
-	return dpdq * dqdr;
+	return dpdq(r, profile_radius, pmodel) * dqdr;
 }
 
 
 /* Returns dp/drad at "r" */
-static double partiality_rgradient(double r, double profile_radius)
+static double partiality_rgradient(double r, double profile_radius,
+                                   PartialityModel pmodel)
 {
-	double q, dpdq, dqdrad;
-
-	/* Calculate degree of penetration */
-	q = (r + profile_radius)/(2.0*profile_radius);
-
-	/* dp/dq */
-	dpdq = 6.0*(q-pow(q, 2.0));
+	double dqdrad;
 
 	/* dq/drad */
 	dqdrad = -0.5 * r * pow(profile_radius, -2.0);
 
-	return dpdq * dqdrad;
+	return dpdq(r, profile_radius, pmodel) * dqdrad;
 }
 
 
@@ -144,12 +160,12 @@ double p_gradient(Crystal *cr, int k, Reflection *refl, PartialityModel pmodel)
 
 	/* Calculate the gradient of partiality wrt excitation error. */
 	if ( clamp_low == 0 ) {
-		glow = partiality_gradient(rlow, r);
+		glow = partiality_gradient(rlow, r, pmodel);
 	} else {
 		glow = 0.0;
 	}
 	if ( clamp_high == 0 ) {
-		ghigh = partiality_gradient(rhigh, r);
+		ghigh = partiality_gradient(rhigh, r, pmodel);
 	} else {
 		ghigh = 0.0;
 	}
@@ -163,8 +179,8 @@ double p_gradient(Crystal *cr, int k, Reflection *refl, PartialityModel pmodel)
 		return (ds*glow + ds*ghigh) / 2.0;
 
 		case REF_R :
-		gr  = partiality_rgradient(rlow, r);
-		gr -= partiality_rgradient(rhigh, r);
+		gr  = partiality_rgradient(rlow, r, pmodel);
+		gr -= partiality_rgradient(rhigh, r, pmodel);
 		return gr;
 
 		/* Cell parameters and orientation */
diff --git a/tests/pr_l_gradient_check.c b/tests/pr_l_gradient_check.c
index 4b3894ba..72d353f9 100644
--- a/tests/pr_l_gradient_check.c
+++ b/tests/pr_l_gradient_check.c
@@ -279,7 +279,6 @@ int main(int argc, char *argv[])
 	Crystal *cr;
 	struct quaternion orientation;
 	int i;
-	const PartialityModel pmodel = PMODEL_SPHERE;
 	int fail = 0;
 	int quiet = 0;
 	int plot = 0;
@@ -337,10 +336,19 @@ int main(int argc, char *argv[])
 
 	rng = gsl_rng_alloc(gsl_rng_mt19937);
 
-	for ( i=0; i<1; i++ ) {
+	for ( i=0; i<2; i++ ) {
 
 		UnitCell *rot;
 		double val;
+		PartialityModel pmodel;
+
+		if ( i == 0 ) {
+			pmodel = PMODEL_SPHERE;
+			STATUS("Testing flat sphere model:\n");
+		} else {
+			pmodel = PMODEL_GAUSSIAN;
+			STATUS("Testing Gaussian model:\n");
+		}
 
 		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 c0582036..c3b55cfa 100644
--- a/tests/pr_p_gradient_check.c
+++ b/tests/pr_p_gradient_check.c
@@ -392,7 +392,6 @@ int main(int argc, char *argv[])
 	Crystal *cr;
 	struct quaternion orientation;
 	int i;
-	const PartialityModel pmodel = PMODEL_SPHERE;
 	int fail = 0;
 	int quiet = 0;
 	int plot = 0;
@@ -450,10 +449,19 @@ int main(int argc, char *argv[])
 
 	rng = gsl_rng_alloc(gsl_rng_mt19937);
 
-	for ( i=0; i<1; i++ ) {
+	for ( i=0; i<2; i++ ) {
 
 		UnitCell *rot;
 		double val;
+		PartialityModel pmodel;
+
+		if ( i == 0 ) {
+			pmodel = PMODEL_SPHERE;
+			STATUS("Testing flat sphere model:\n");
+		} else {
+			pmodel = PMODEL_GAUSSIAN;
+			STATUS("Testing Gaussian model:\n");
+		}
 
 		orientation = random_quaternion(rng);
 		rot = cell_rotate(cell, orientation);
diff --git a/tests/pr_pl_gradient_check.c b/tests/pr_pl_gradient_check.c
index 785d9973..ddbc1840 100644
--- a/tests/pr_pl_gradient_check.c
+++ b/tests/pr_pl_gradient_check.c
@@ -394,7 +394,6 @@ int main(int argc, char *argv[])
 	Crystal *cr;
 	struct quaternion orientation;
 	int i;
-	const PartialityModel pmodel = PMODEL_SPHERE;
 	int fail = 0;
 	int quiet = 0;
 	int plot = 0;
@@ -452,10 +451,19 @@ int main(int argc, char *argv[])
 
 	rng = gsl_rng_alloc(gsl_rng_mt19937);
 
-	for ( i=0; i<1; i++ ) {
+	for ( i=0; i<2; i++ ) {
 
 		UnitCell *rot;
 		double val;
+		PartialityModel pmodel;
+
+		if ( i == 0 ) {
+			pmodel = PMODEL_SPHERE;
+			STATUS("Testing flat sphere model:\n");
+		} else {
+			pmodel = PMODEL_GAUSSIAN;
+			STATUS("Testing Gaussian model:\n");
+		}
 
 		orientation = random_quaternion(rng);
 		rot = cell_rotate(cell, orientation);