Add Gaussian partiality model

3 files changed, 41 insertions, 74 deletions

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 */