Merge branch 'tom/partials'

19 files changed, 872 insertions, 196 deletions

diff --git a/Makefile.am b/Makefile.am
index 5632fb53..c0bf3c59 100644
--- a/Makefile.am
+++ b/Makefile.am
@@ -11,14 +11,19 @@ noinst_PROGRAMS = tests/list_check tests/integration_check \
                   tests/centering_check tests/transformation_check \
                   tests/cell_check
 
-TESTS = tests/list_check tests/first_merge_check tests/second_merge_check \
-        tests/third_merge_check tests/fourth_merge_check \
+MERGE_CHECKS = tests/first_merge_check tests/second_merge_check \
+               tests/third_merge_check tests/fourth_merge_check
+
+PARTIAL_CHECKS = tests/partialator_merge_check_1 \
+                 tests/partialator_merge_check_2 \
+                 tests/partialator_merge_check_3
+
+TESTS = tests/list_check $(MERGE_CHECKS) $(PARTIAL_CHECKS) \
         tests/integration_check tests/pr_gradient_check \
         tests/symmetry_check tests/centering_check tests/transformation_check \
         tests/cell_check
 
-EXTRA_DIST += tests/first_merge_check tests/second_merge_check \
-              tests/third_merge_check tests/fourth_merge_check
+EXTRA_DIST += $(MERGE_CHECKS) $(PARTIAL_CHECKS)
 
 if BUILD_HDFSEE
 bin_PROGRAMS += src/hdfsee
diff --git a/doc/man/partial_sim.1 b/doc/man/partial_sim.1
index 5b417978..7321b9ff 100644
--- a/doc/man/partial_sim.1
+++ b/doc/man/partial_sim.1
@@ -1,7 +1,7 @@
 .\"
 .\" partial_sim man page
 .\"
-.\" Copyright © 2012 Thomas White <taw@physics.org>
+.\" Copyright © 2012-2013 Thomas White <taw@physics.org>
 .\"
 .\" Part of CrystFEL - crystallography with a FEL
 .\"
@@ -77,6 +77,24 @@ Add random values with a flat distribution to the components of the reciprocal l
 .PD
 Save a table of values to \fIfilename\fR containing, in resolution shells, the following columns: resolution shell centre in nm^-1, number of reflections in shell, mean partiality, maximum partiality.
 
+.PD 0
+.B
+.IP "\fB--osf-stddev=\fR\fIval\fR
+.PD
+Set the standard deviation of the distribution of overall scaling factors to \fIval\fR.  The distribution will be cut at zero, i.e. negative or zero scaling factors are not allowed.  The distribution will be Gaussian centered on 1.  The default is \fB--osf-stddev=2.0\fR.
+
+.PD 0
+.B
+.IP "\fB--full-stddev=\fR\fIval\fR
+.PD
+Set the standard deviation of the distribution of randomly generated full intensities to \fIval\fR.  The distribution will be Gaussian, centered on zero, and the absolute value will be taken (i.e. there will be no negative values).  The default is \fB--full-stddev=1000.0\fR.  This option has no effect if you also use \fB-i\fR or \fR--input\fB.
+
+.PD 0
+.B
+.IP "\fB--noise-stddev=\fR\fIval\fR
+.PD
+Set the standard deviation of the noise added to the partial intensities to \fIval\fR.  The noise will be Gaussian, and the same for all reflections.  The default is \fB--noise-stddev=20.0\fR.
+
 .SH AUTHOR
 This page was written by Thomas White.
 
@@ -84,7 +102,7 @@ This page was written by Thomas White.
 Report bugs to <taw@physics.org>, or visit <http://www.desy.de/~twhite/crystfel>.
 
 .SH COPYRIGHT AND DISCLAIMER
-Copyright © 2012 Deutsches Elektronen-Synchrotron DESY, a research centre of the Helmholtz Association.
+Copyright © 2012-2013 Deutsches Elektronen-Synchrotron DESY, a research centre of the Helmholtz Association.
 .P
 partial_sim, and this manaul, are part of CrystFEL.
 .P
diff --git a/doc/reference/libcrystfel/CrystFEL-sections.txt b/doc/reference/libcrystfel/CrystFEL-sections.txt
index 720cabff..09267538 100644
--- a/doc/reference/libcrystfel/CrystFEL-sections.txt
+++ b/doc/reference/libcrystfel/CrystFEL-sections.txt
@@ -21,10 +21,11 @@ next_found_refl
 get_excitation_error
 get_detector_pos
 get_partiality
+get_lorentz
+get_partial
 get_indices
 get_symmetric_indices
 get_intensity
-get_partial
 get_scalable
 get_refinable
 get_redundancy
@@ -34,6 +35,8 @@ get_temp1
 get_temp2
 <SUBSECTION>
 set_detector_pos
+set_partiality
+set_lorentz
 set_partial
 set_intensity
 set_scalable
@@ -380,6 +383,7 @@ crystal_set_user_flag
 
 <SECTION>
 <FILE>geometry</FILE>
+PartialityModel
 find_intersections
 select_intersections
 update_partialities
diff --git a/libcrystfel/src/geometry.c b/libcrystfel/src/geometry.c
index 621c9484..cdc62ab5 100644
--- a/libcrystfel/src/geometry.c
+++ b/libcrystfel/src/geometry.c
@@ -134,6 +134,10 @@ static Reflection *check_reflection(struct image *image, Crystal *cryst,
 	Reflection *refl;
 	double cet, cez;
 	double pr;
+	double L;
+
+	/* Don't predict 000 */
+	if ( abs(h)+abs(k)+abs(l) == 0 ) return NULL;
 
 	pr = crystal_get_profile_radius(cryst);
 
@@ -160,6 +164,19 @@ static Reflection *check_reflection(struct image *image, Crystal *cryst,
 	     && (fabs(rlow) > pr)
 	     && (fabs(rhigh) > pr) ) return NULL;
 
+	if ( rlow < rhigh ) {
+		ERROR("Reflection with rlow < rhigh!\n");
+		ERROR("%3i %3i %3i  rlow = %e, rhigh = %e\n",
+		      h, k, l, rlow, rhigh);
+		ERROR("div = %e\n", image->div);
+		return NULL;
+	}
+
+	/* Lorentz factor is determined direction from the r values, before
+	 * clamping.  The multiplication by the profile radius is to make the
+	 * correction factor vaguely near 1. */
+	L = pr / (rlow - rhigh);
+
 	/* If the "lower" Ewald sphere is a long way away, use the
 	 * position at which the Ewald sphere would just touch the
 	 * reflection.
@@ -185,7 +202,6 @@ static Reflection *check_reflection(struct image *image, Crystal *cryst,
 		rhigh = +pr;
 		clamp_high = +1;
 	}
-	assert(clamp_low >= clamp_high);
 
 	/* Calculate partiality */
 	part = partiality(rlow, rhigh, pr);
@@ -198,6 +214,7 @@ static Reflection *check_reflection(struct image *image, Crystal *cryst,
 	refl = reflection_new(h, k, l);
 	set_detector_pos(refl, 0.0, xda, yda);
 	set_partial(refl, rlow, rhigh, part, clamp_low, clamp_high);
+	set_lorentz(refl, L);
 	set_symmetric_indices(refl, h, k, l);
 	set_redundancy(refl, 1);
 
@@ -345,8 +362,23 @@ RefList *select_intersections(struct image *image, Crystal *cryst)
 }
 
 
+static void set_unity_partialities(Crystal *cryst)
+{
+	Reflection *refl;
+	RefListIterator *iter;
+
+	for ( refl = first_refl(crystal_get_reflections(cryst), &iter);
+	      refl != NULL;
+	      refl = next_refl(refl, iter) )
+	{
+		set_partiality(refl, 1.0);
+		set_lorentz(refl, 1.0);
+	}
+}
+
+
 /* Calculate partialities and apply them to the image's reflections */
-void update_partialities(Crystal *cryst)
+void update_partialities(Crystal *cryst, PartialityModel pmodel)
 {
 	Reflection *refl;
 	RefListIterator *iter;
@@ -356,6 +388,14 @@ void update_partialities(Crystal *cryst)
 	double csx, csy, csz;
 	struct image *image = crystal_get_image(cryst);
 
+	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;
+	}
+
 	cell_get_reciprocal(crystal_get_cell(cryst), &asx, &asy, &asz,
 	                    &bsx, &bsy, &bsz, &csx, &csy, &csz);
 
@@ -367,7 +407,7 @@ void update_partialities(Crystal *cryst)
 	      refl = next_refl(refl, iter) )
 	{
 		Reflection *vals;
-		double r1, r2, p, x, y;
+		double r1, r2, L, p, x, y;
 		double xl, yl, zl;
 		signed int h, k, l;
 		int clamp1, clamp2;
@@ -390,6 +430,8 @@ void update_partialities(Crystal *cryst)
 		/* Transfer partiality stuff */
 		get_partial(vals, &r1, &r2, &p, &clamp1, &clamp2);
 		set_partial(refl, r1, r2, p, clamp1, clamp2);
+		L = get_lorentz(vals);
+		set_lorentz(refl, L);
 
 		/* Transfer detector location */
 		get_detector_pos(vals, &x, &y);
diff --git a/libcrystfel/src/geometry.h b/libcrystfel/src/geometry.h
index aecdc28a..2ac5982b 100644
--- a/libcrystfel/src/geometry.h
+++ b/libcrystfel/src/geometry.h
@@ -43,10 +43,26 @@
 extern "C" {
 #endif
 
+/**
+ * PartialityModel:
+ * @PMODEL_SPHERE : Intersection of sphere with excited volume of reciprocal
+ *   space.
+ * @PMODEL_UNITY : Set all all partialities and Lorentz factors to 1.
+ *
+ * A %PartialityModel describes a geometrical model which can be used to
+ * calculate spot partialities and Lorentz correction factors.
+ **/
+typedef enum {
+
+	PMODEL_SPHERE,
+	PMODEL_UNITY,
+
+} PartialityModel;
+
 extern RefList *find_intersections(struct image *image, Crystal *cryst);
 extern RefList *select_intersections(struct image *image, Crystal *cryst);
 
-extern void update_partialities(Crystal *cryst);
+extern void update_partialities(Crystal *cryst, PartialityModel pmodel);
 
 #ifdef __cplusplus
 }
diff --git a/libcrystfel/src/reflist.c b/libcrystfel/src/reflist.c
index 64154774..b3b9f85b 100644
--- a/libcrystfel/src/reflist.c
+++ b/libcrystfel/src/reflist.c
@@ -3,11 +3,11 @@
  *
  * Fast reflection/peak list
  *
- * Copyright © 2012 Deutsches Elektronen-Synchrotron DESY,
- *                  a research centre of the Helmholtz Association.
+ * Copyright © 2012-2013 Deutsches Elektronen-Synchrotron DESY,
+ *                       a research centre of the Helmholtz Association.
  *
  * Authors:
- *   2011-2012 Thomas White <taw@physics.org>
+ *   2011-2013 Thomas White <taw@physics.org>
  *
  * This file is part of CrystFEL.
  *
@@ -70,6 +70,7 @@ struct _refldata {
 	double r1;  /* First excitation error */
 	double r2;  /* Second excitation error */
 	double p;   /* Partiality */
+	double L;   /* Lorentz factor */
 	int clamp1; /* Clamp status for r1 */
 	int clamp2; /* Clamp status for r2 */
 
@@ -388,7 +389,7 @@ void get_symmetric_indices(const Reflection *refl,
  * get_partiality:
  * @refl: A %Reflection
  *
- * Returns: The partiality of the reflection.
+ * Returns: The partiality of the reflection.  See get_lorentz().
  **/
 double get_partiality(const Reflection *refl)
 {
@@ -397,6 +398,19 @@ double get_partiality(const Reflection *refl)
 
 
 /**
+ * get_lorentz:
+ * @refl: A %Reflection
+ *
+ * Returns: The Lorentz factor for the reflection.  To "scale up" a partial
+ * reflection, divide by this multiplied by the partiality.
+ **/
+double get_lorentz(const Reflection *refl)
+{
+	return refl->data.L;
+}
+
+
+/**
  * get_intensity:
  * @refl: A %Reflection
  *
@@ -602,10 +616,35 @@ void set_partial(Reflection *refl, double r1, double r2, double p,
 /**
  * set_intensity:
  * @refl: A %Reflection
+ * @p: The partiality for the reflection.
+ *
+ * Set the partiality for the reflection.  See set_lorentz().
+ **/
+void set_partiality(Reflection *refl, double p)
+{
+	refl->data.p = p;
+}
+
+/**
+ * set_lorentz:
+ * @refl: A %Reflection
+ * @L: The Lorentz factor for the reflection.
+ *
+ * Set the Lorentz factor for the reflection.  To "scale up" a partial
+ * reflection, divide by this multiplied by the partiality.
+ **/
+void set_lorentz(Reflection *refl, double L)
+{
+	refl->data.L = L;
+}
+
+
+/**
+ * set_intensity:
+ * @refl: A %Reflection
  * @intensity: The intensity for the reflection.
  *
- * Set the intensity for the reflection.  Note that retrieval is done with
- * get_intensity().
+ * Set the intensity for the reflection.
  **/
 void set_intensity(Reflection *refl, double intensity)
 {
diff --git a/libcrystfel/src/reflist.h b/libcrystfel/src/reflist.h
index 246ef885..7969235c 100644
--- a/libcrystfel/src/reflist.h
+++ b/libcrystfel/src/reflist.h
@@ -3,11 +3,11 @@
  *
  * Fast reflection/peak list
  *
- * Copyright © 2012 Deutsches Elektronen-Synchrotron DESY,
- *                  a research centre of the Helmholtz Association.
+ * Copyright © 2012-2013 Deutsches Elektronen-Synchrotron DESY,
+ *                       a research centre of the Helmholtz Association.
  *
  * Authors:
- *   2011-2012 Thomas White <taw@physics.org>
+ *   2011-2013 Thomas White <taw@physics.org>
  *
  * This file is part of CrystFEL.
  *
@@ -78,6 +78,7 @@ extern Reflection *next_found_refl(Reflection *refl);
 extern double get_excitation_error(const Reflection *refl);
 extern void get_detector_pos(const Reflection *refl, double *fs, double *ss);
 extern double get_partiality(const Reflection *refl);
+extern double get_lorentz(const Reflection *refl);
 extern void get_indices(const Reflection *refl,
                         signed int *h, signed int *k, signed int *l);
 extern void get_symmetric_indices(const Reflection *refl,
@@ -100,6 +101,8 @@ 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_partiality(Reflection *refl, double p);
+extern void set_lorentz(Reflection *refl, double L);
 extern void set_intensity(Reflection *refl, double intensity);
 extern void set_scalable(Reflection *refl, int scalable);
 extern void set_refinable(Reflection *refl, int refinable);
diff --git a/libcrystfel/src/utils.h b/libcrystfel/src/utils.h
index b75693db..1adb69e6 100644
--- a/libcrystfel/src/utils.h
+++ b/libcrystfel/src/utils.h
@@ -139,6 +139,11 @@ static inline double modulus(double x, double y, double z)
 	return sqrt(x*x + y*y + z*z);
 }
 
+static inline double modulus2d(double x, double y)
+{
+	return sqrt(x*x + y*y);
+}
+
 static inline double modulus_squared(double x, double y, double z) {
 	return x*x + y*y + z*z;
 }
@@ -165,6 +170,21 @@ static inline double angle_between(double x1, double y1, double z1,
 	return acos(cosine);
 }
 
+/* Answer in radians */
+static inline double angle_between_2d(double x1, double y1,
+                                      double x2, double y2)
+{
+	double mod1 = modulus2d(x1, y1);
+	double mod2 = modulus2d(x2, y2);
+	double cosine = (x1*x2 + y1*y2) / (mod1*mod2);
+
+	/* Fix domain if outside due to rounding */
+	if ( cosine > 1.0 ) cosine = 1.0;
+	if ( cosine < -1.0 ) cosine = -1.0;
+
+	return acos(cosine);
+}
+
 static inline int within_tolerance(double a, double b, double percent)
 {
 	double tol = fabs(a) * (percent/100.0);
diff --git a/src/hrs-scaling.c b/src/hrs-scaling.c
index 0498a532..c249ef46 100644
--- a/src/hrs-scaling.c
+++ b/src/hrs-scaling.c
@@ -3,11 +3,11 @@
  *
  * Intensity scaling using generalised HRS target function
  *
- * Copyright © 2012 Deutsches Elektronen-Synchrotron DESY,
- *                  a research centre of the Helmholtz Association.
+ * Copyright © 2012-2013 Deutsches Elektronen-Synchrotron DESY,
+ *                       a research centre of the Helmholtz Association.
  *
  * Authors:
- *   2010-2012 Thomas White <taw@physics.org>
+ *   2010-2013 Thomas White <taw@physics.org>
  *
  * This file is part of CrystFEL.
  *
@@ -60,6 +60,7 @@ struct scale_queue_args
 	Crystal **crystals;
 	int n_started;
 	double max_shift;
+	PartialityModel pmodel;
 };
 
 
@@ -68,6 +69,7 @@ struct scale_worker_args
 	Crystal *crystal;
 	double shift;
 	RefList *reference;
+	PartialityModel pmodel;
 };
 
 
@@ -78,6 +80,7 @@ static void *create_scale_job(void *vqargs)
 
 	wargs = malloc(sizeof(struct scale_worker_args));
 	wargs->reference = qargs->reference;
+	wargs->pmodel = qargs->pmodel;
 
 	wargs->crystal = qargs->crystals[qargs->n_started++];
 
@@ -94,8 +97,8 @@ static void run_scale_job(void *vwargs, int cookie)
 	RefListIterator *iter;
 	double num = 0.0;
 	double den = 0.0;
-	double corr;
-	const double osf = crystal_get_osf(cr);
+	double g;
+	const double G = crystal_get_osf(cr);
 
 	if ( crystal_get_user_flag(cr) ) {
 		wargs->shift = 0.0;
@@ -109,6 +112,7 @@ static void run_scale_job(void *vwargs, int cookie)
 		signed int h, k, l;
 		double Ih, Ihl, esd;
 		Reflection *r;
+		double corr;
 
 		if ( !get_scalable(refl) ) continue;
 
@@ -125,24 +129,41 @@ static void run_scale_job(void *vwargs, int cookie)
 			Ih = get_intensity(r);
 		}
 
-		/* FIXME: This isn't correct */
-		Ihl = get_intensity(refl) / get_partiality(refl);
-		esd = get_esd_intensity(refl) / get_partiality(refl);
+		/* If you change this, be sure to also change
+		 * run_merge_job() and run_esd_job(). */
+		switch ( wargs->pmodel ) {
 
-		num += Ih * (Ihl/osf) / pow(esd/osf, 2.0);
-		den += pow(Ih, 2.0)/pow(esd/osf, 2.0);
+			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;
+
+		}
+
+		Ihl = get_intensity(refl) / corr;
+		esd = get_esd_intensity(refl) / corr;
+
+		num += Ih * (Ihl/G) / pow(esd/G, 2.0);
+		den += pow(Ih, 2.0)/pow(esd/G, 2.0);
 
 	}
 
 	//num += image->osf / pow(SCALING_RESTRAINT, 2.0);
 	//den += pow(image->osf, 2.0)/pow(SCALING_RESTRAINT, 2.0);
 
-	corr = num / den;
-	if ( !isnan(corr) && !isinf(corr) ) {
-		crystal_set_osf(cr, osf*corr);
+	g = num / den;
+	if ( !isnan(g) && !isinf(g) ) {
+		crystal_set_osf(cr, g*G);
 	}
-	wargs->shift = fabs(corr-1.0);
-
+	wargs->shift = fabs(g-1.0);
 }
 
 
@@ -157,7 +178,7 @@ static void finalise_scale_job(void *vqargs, void *vwargs)
 
 
 static double iterate_scale(Crystal **crystals, int n, RefList *reference,
-                            int n_threads)
+                            int n_threads, PartialityModel pmodel)
 {
 	struct scale_queue_args qargs;
 
@@ -167,6 +188,7 @@ static double iterate_scale(Crystal **crystals, int n, RefList *reference,
 	qargs.n_started = 0;
 	qargs.crystals = crystals;
 	qargs.max_shift = 0.0;
+	qargs.pmodel = pmodel;
 
 	run_threads(n_threads, run_scale_job, create_scale_job,
 	            finalise_scale_job, &qargs, n, 0, 0, 0);
@@ -181,6 +203,7 @@ struct merge_queue_args
 	pthread_mutex_t full_lock;
 	Crystal **crystals;
 	int n_started;
+	PartialityModel pmodel;
 };
 
 
@@ -189,6 +212,7 @@ struct merge_worker_args
 	Crystal *crystal;
 	RefList *full;
 	pthread_mutex_t *full_lock;
+	PartialityModel pmodel;
 };
 
 
@@ -200,6 +224,7 @@ static void *create_merge_job(void *vqargs)
 	wargs = malloc(sizeof(struct merge_worker_args));
 	wargs->full = qargs->full;
 	wargs->full_lock = &qargs->full_lock;
+	wargs->pmodel = qargs->pmodel;
 
 	wargs->crystal = qargs->crystals[qargs->n_started++];
 
@@ -228,7 +253,7 @@ static void run_merge_job(void *vwargs, int cookie)
 		signed int h, k, l;
 		double num, den;
 		int red;
-		double Ihl, esd, pcalc;
+		double Ihl, esd, corr;
 
 		if ( !get_scalable(refl) ) continue;
 
@@ -251,9 +276,27 @@ static void run_merge_job(void *vwargs, int cookie)
 			red = get_redundancy(f);
 		}
 
-		pcalc = get_partiality(refl);
-		Ihl = get_intensity(refl) / pcalc;
-		esd = get_esd_intensity(refl) / pcalc;
+		/* 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;
+
+		}
+
+		Ihl = get_intensity(refl) / corr;
+		esd = get_esd_intensity(refl) / corr;
 
 		num += (Ihl/G) / pow(esd/G, 2.0);
 		den += 1.0 / pow(esd/G, 2.0);
@@ -273,7 +316,8 @@ static void finalise_merge_job(void *vqargs, void *vwargs)
 }
 
 
-static RefList *lsq_intensities(Crystal **crystals, int n, int n_threads)
+static RefList *lsq_intensities(Crystal **crystals, int n, int n_threads,
+                                PartialityModel pmodel)
 {
 	RefList *full;
 	struct merge_queue_args qargs;
@@ -285,6 +329,7 @@ static RefList *lsq_intensities(Crystal **crystals, int n, int n_threads)
 	qargs.full = full;
 	qargs.n_started = 0;
 	qargs.crystals = crystals;
+	qargs.pmodel = pmodel;
 	pthread_mutex_init(&qargs.full_lock, NULL);
 
 	run_threads(n_threads, run_merge_job, create_merge_job,
@@ -311,6 +356,7 @@ struct esd_queue_args
 	RefList *full;
 	Crystal **crystals;
 	int n_started;
+	PartialityModel pmodel;
 };
 
 
@@ -318,6 +364,7 @@ struct esd_worker_args
 {
 	Crystal *crystal;
 	RefList *full;
+	PartialityModel pmodel;
 };
 
 
@@ -328,6 +375,7 @@ static void *create_esd_job(void *vqargs)
 
 	wargs = malloc(sizeof(struct esd_worker_args));
 	wargs->full = qargs->full;
+	wargs->pmodel = qargs->pmodel;
 
 	wargs->crystal = qargs->crystals[qargs->n_started++];
 
@@ -355,7 +403,7 @@ static void run_esd_job(void *vwargs, int cookie)
 		Reflection *f;
 		signed int h, k, l;
 		double num;
-		double Ihl, Ih;
+		double Ihl, Ih, corr;
 
 		if ( !get_scalable(refl) ) continue;
 
@@ -367,8 +415,27 @@ 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;
+
+		}
+
 		Ih = get_intensity(f);
-		Ihl = get_intensity(refl) / (get_partiality(refl) * G);
+		Ihl = get_intensity(refl) / (corr * G);
 
 		num += pow(Ihl - Ih, 2.0);
 
@@ -385,7 +452,7 @@ static void finalise_esd_job(void *vqargs, void *vwargs)
 
 
 static void calculate_esds(Crystal **crystals, int n, RefList *full,
-                           int n_threads, int min_red)
+                           int n_threads, int min_red, PartialityModel pmodel)
 {
 	struct esd_queue_args qargs;
 	Reflection *refl;
@@ -394,6 +461,7 @@ static void calculate_esds(Crystal **crystals, int n, RefList *full,
 	qargs.full = full;
 	qargs.n_started = 0;
 	qargs.crystals = crystals;
+	qargs.pmodel = pmodel;
 
 	for ( refl = first_refl(full, &iter);
 	      refl != NULL;
@@ -425,24 +493,25 @@ static void calculate_esds(Crystal **crystals, int n, RefList *full,
 
 /* Scale the stack of images */
 RefList *scale_intensities(Crystal **crystals, int n, RefList *gref,
-                           int n_threads, int noscale)
+                           int n_threads, int noscale, PartialityModel pmodel,
+                           int min_redundancy)
 {
 	int i;
 	double max_corr;
 	RefList *full = NULL;
-	const int min_redundancy = 3;
 
 	for ( i=0; i<n; i++ ) crystal_set_osf(crystals[i], 1.0);
 
 	if ( noscale ) {
-		full = lsq_intensities(crystals, n, n_threads);
-		calculate_esds(crystals, n, full, n_threads, min_redundancy);
+		full = lsq_intensities(crystals, n, n_threads, pmodel);
+		calculate_esds(crystals, n, full, n_threads, min_redundancy,
+		               pmodel);
 		return full;
 	}
 
 	/* No reference -> create an initial list to refine against */
 	if ( gref == NULL ) {
-		full = lsq_intensities(crystals, n, n_threads);
+		full = lsq_intensities(crystals, n, n_threads, pmodel);
 	}
 
 	/* Iterate */
@@ -458,27 +527,26 @@ RefList *scale_intensities(Crystal **crystals, int n, RefList *gref,
 			reference = full;
 		}
 
-		max_corr = iterate_scale(crystals, n, reference, n_threads);
+		max_corr = iterate_scale(crystals, n, reference, n_threads,
+		                         pmodel);
 		//STATUS("Scaling iteration %2i: max correction = %5.2f\n",
 		//       i+1, max_corr);
 
 		/* No reference -> generate list for next iteration */
 		if ( gref == NULL ) {
 			reflist_free(full);
-			full = lsq_intensities(crystals, n, n_threads);
+			full = lsq_intensities(crystals, n, n_threads, pmodel);
 		}
 
-		//show_scale_factors(images, n);
-
 		i++;
 
 	} while ( (max_corr > 0.01) && (i < MAX_CYCLES) );
 
 	if ( gref != NULL ) {
-		full = lsq_intensities(crystals, n, n_threads);
+		full = lsq_intensities(crystals, n, n_threads, pmodel);
 	} /* else we already did it */
 
-	calculate_esds(crystals, n, full, n_threads, min_redundancy);
+	calculate_esds(crystals, n, full, n_threads, min_redundancy, pmodel);
 
 	return full;
 }
diff --git a/src/hrs-scaling.h b/src/hrs-scaling.h
index 80940347..0979bb98 100644
--- a/src/hrs-scaling.h
+++ b/src/hrs-scaling.h
@@ -3,11 +3,11 @@
  *
  * Intensity scaling using generalised HRS target function
  *
- * Copyright © 2012 Deutsches Elektronen-Synchrotron DESY,
- *                  a research centre of the Helmholtz Association.
+ * Copyright © 2012-2013 Deutsches Elektronen-Synchrotron DESY,
+ *                       a research centre of the Helmholtz Association.
  *
  * Authors:
- *   2010-2012 Thomas White <taw@physics.org>
+ *   2010-2013 Thomas White <taw@physics.org>
  *
  * This file is part of CrystFEL.
  *
@@ -37,10 +37,12 @@
 
 #include "crystal.h"
 #include "reflist.h"
+#include "geometry.h"
 
 extern RefList *scale_intensities(Crystal **crystals, int n,
                                   RefList *reference, int n_threads,
-                                  int noscale);
+                                  int noscale, PartialityModel pmodel,
+                                  int min_redundancy);
 
 
 #endif	/* HRS_SCALING_H */
diff --git a/src/partial_sim.c b/src/partial_sim.c
index 7c93da1b..87d30f39 100644
--- a/src/partial_sim.c
+++ b/src/partial_sim.c
@@ -3,11 +3,11 @@
  *
  * Generate partials for testing scaling
  *
- * Copyright © 2012 Deutsches Elektronen-Synchrotron DESY,
- *                  a research centre of the Helmholtz Association.
+ * Copyright © 2012-2013 Deutsches Elektronen-Synchrotron DESY,
+ *                       a research centre of the Helmholtz Association.
  *
  * Authors:
- *   2011-2012 Thomas White <taw@physics.org>
+ *   2011-2013 Thomas White <taw@physics.org>
  *
  * This file is part of CrystFEL.
  *
@@ -88,7 +88,8 @@ static void calculate_partials(Crystal *cr,
                                int random_intensities,
                                pthread_mutex_t *full_lock,
                                unsigned long int *n_ref, double *p_hist,
-                               double *p_max, double max_q)
+                               double *p_max, double max_q, double full_stddev,
+                               double noise_stddev)
 {
 	Reflection *refl;
 	RefListIterator *iter;
@@ -100,12 +101,13 @@ static void calculate_partials(Crystal *cr,
 	{
 		signed int h, k, l;
 		Reflection *rfull;
-		double p, Ip, If;
+		double L, p, Ip, If;
 		int bin;
 
 		get_indices(refl, &h, &k, &l);
 		get_asymm(sym, h, k, l, &h, &k, &l);
 		p = get_partiality(refl);
+		L = get_lorentz(refl);
 
 		pthread_mutex_lock(full_lock);
 		rfull = find_refl(full, h, k, l);
@@ -120,7 +122,7 @@ static void calculate_partials(Crystal *cr,
 				 * thing under lock. */
 				pthread_mutex_lock(full_lock);
 				rfull = add_refl(full, h, k, l);
-				If = fabs(gaussian_noise(0.0, 1000.0));
+				If = fabs(gaussian_noise(0.0, full_stddev));
 				set_intensity(rfull, If);
 				set_redundancy(rfull, 1);
 				pthread_mutex_unlock(full_lock);
@@ -139,7 +141,7 @@ static void calculate_partials(Crystal *cr,
 			}
 		}
 
-		Ip = crystal_get_osf(cr) * p * If;
+		Ip = crystal_get_osf(cr) * L * p * If;
 
 		res = resolution(crystal_get_cell(cr), h, k, l);
 		bin = NBINS*2.0*res/max_q;
@@ -152,10 +154,10 @@ static void calculate_partials(Crystal *cr,
 			       res, bin,  p);
 		}
 
-		Ip = gaussian_noise(Ip, 100.0);
+		Ip = gaussian_noise(Ip, noise_stddev);
 
 		set_intensity(refl, Ip);
-		set_esd_intensity(refl, 100.0);
+		set_esd_intensity(refl, noise_stddev);
 	}
 }
 
@@ -183,6 +185,10 @@ static void show_help(const char *s)
 " -c, --cnoise=<val>       Add random noise, with a flat distribution, to the\n"
 "                          reciprocal lattice vector components given in the\n"
 "                          stream, with maximum error +/- <val> percent.\n"
+"     --osf-stddev=<val>   Set the standard deviation of the scaling factors.\n"
+"     --full-stddev=<val>  Set the standard deviation of the randomly\n"
+"                           generated full intensities, if not using -i.\n"
+"     --noise-stddev=<val>  Set the standard deviation of the noise.\n"
 "\n"
 );
 }
@@ -201,6 +207,9 @@ struct queue_args
 	int random_intensities;
 	UnitCell *cell;
 	double cnoise;
+	double osf_stddev;
+	double full_stddev;
+	double noise_stddev;
 
 	struct image *template_image;
 	double max_q;
@@ -254,6 +263,7 @@ static void run_job(void *vwargs, int cookie)
 	int i;
 	Crystal *cr;
 	RefList *reflections;
+	double osf;
 
 	cr = crystal_new();
 	if ( cr == NULL ) {
@@ -263,7 +273,10 @@ static void run_job(void *vwargs, int cookie)
 	wargs->crystal = cr;
 	crystal_set_image(cr, &wargs->image);
 
-	crystal_set_osf(cr, gaussian_noise(1.0, 0.3));
+	do {
+		osf = gaussian_noise(1.0, qargs->osf_stddev);
+	} while ( osf <= 0.0 );
+	crystal_set_osf(cr, osf);
 	crystal_set_profile_radius(cr, wargs->image.beam->profile_radius);
 
 	/* Set up a random orientation */
@@ -284,7 +297,8 @@ static void run_job(void *vwargs, int cookie)
 	                   qargs->sym, qargs->random_intensities,
 	                   &qargs->full_lock,
 	                   wargs->n_ref, wargs->p_hist, wargs->p_max,
-	                   qargs->max_q);
+	                   qargs->max_q, qargs->full_stddev,
+	                   qargs->noise_stddev);
 
 	/* Give a slightly incorrect cell in the stream */
 	mess_up_cell(cr, qargs->cnoise);
@@ -343,6 +357,9 @@ int main(int argc, char *argv[])
 	int i;
 	FILE *fh;
 	char *phist_file = NULL;
+	double osf_stddev = 2.0;
+	double full_stddev = 1000.0;
+	double noise_stddev = 20.0;
 
 	/* Long options */
 	const struct option longopts[] = {
@@ -354,8 +371,13 @@ int main(int argc, char *argv[])
 		{"geometry",           1, NULL,               'g'},
 		{"symmetry",           1, NULL,               'y'},
 		{"save-random",        1, NULL,               'r'},
-		{"pgraph",             1, NULL,                2},
 		{"cnoise",             1, NULL,               'c'},
+
+		{"pgraph",             1, NULL,                2},
+		{"osf-stddev",         1, NULL,                3},
+		{"full-stddev",        1, NULL,                4},
+		{"noise-stddev",       1, NULL,                5},
+
 		{0, 0, NULL, 0}
 	};
 
@@ -417,6 +439,45 @@ int main(int argc, char *argv[])
 			phist_file = strdup(optarg);
 			break;
 
+			case 3 :
+			osf_stddev = strtod(optarg, &rval);
+			if ( *rval != '\0' ) {
+				ERROR("Invalid OSF standard deviation.\n");
+				return 1;
+			}
+			if ( osf_stddev < 0.0 ) {
+				ERROR("Invalid OSF standard deviation.");
+				ERROR(" (must be positive).\n");
+				return 1;
+			}
+			break;
+
+			case 4 :
+			full_stddev = strtod(optarg, &rval);
+			if ( *rval != '\0' ) {
+				ERROR("Invalid full standard deviation.\n");
+				return 1;
+			}
+			if ( full_stddev < 0.0 ) {
+				ERROR("Invalid full standard deviation.");
+				ERROR(" (must be positive).\n");
+				return 1;
+			}
+			break;
+
+			case 5 :
+			noise_stddev = strtod(optarg, &rval);
+			if ( *rval != '\0' ) {
+				ERROR("Invalid noise standard deviation.\n");
+				return 1;
+			}
+			if ( noise_stddev < 0.0 ) {
+				ERROR("Invalid noise standard deviation.");
+				ERROR(" (must be positive).\n");
+				return 1;
+			}
+			break;
+
 			case 0 :
 			break;
 
@@ -548,6 +609,9 @@ int main(int argc, char *argv[])
 	qargs.template_image = &image;
 	qargs.stream = stream;
 	qargs.cnoise = cnoise;
+	qargs.osf_stddev = osf_stddev;
+	qargs.full_stddev = full_stddev;
+	qargs.noise_stddev = noise_stddev;
 	qargs.max_q = largest_q(&image);
 
 	for ( i=0; i<NBINS; i++ ) {
diff --git a/src/partialator.c b/src/partialator.c
index a4af3c18..cfa3a6ba 100644
--- a/src/partialator.c
+++ b/src/partialator.c
@@ -79,6 +79,9 @@ static void show_help(const char *s)
 "  -r, --reference=<file>     Refine images against reflections in <file>,\n"
 "                              instead of taking the mean of the intensity\n"
 "                              estimates.\n"
+"  -m, --model=<model>        Specify partiality model.\n"
+"      --min-measurements=<n> Require at least <n> measurements before a\n"
+"                             reflection appears in the output.  Default: 2\n"
 "\n"
 "  -j <n>                     Run <n> analyses in parallel.\n");
 }
@@ -88,6 +91,7 @@ struct refine_args
 {
 	RefList *full;
 	Crystal *crystal;
+	PartialityModel pmodel;
 };
 
 
@@ -106,7 +110,7 @@ static void refine_image(void *task, int id)
 	struct refine_args *pargs = task;
 	Crystal *cr = pargs->crystal;
 
-	pr_refine(cr, pargs->full);
+	pr_refine(cr, pargs->full, pargs->pmodel);
 }
 
 
@@ -139,13 +143,18 @@ static void done_image(void *vqargs, void *task)
 
 static void refine_all(Crystal **crystals, int n_crystals,
                        struct detector *det,
-                       RefList *full, int nthreads)
+                       RefList *full, int nthreads, PartialityModel pmodel)
 {
 	struct refine_args task_defaults;
 	struct queue_args qargs;
 
+	/* If the partiality model is "p=1", this refinement is really, really
+	 * easy... */
+	if ( pmodel == PMODEL_UNITY ) return;
+
 	task_defaults.full = full;
 	task_defaults.crystal = NULL;
+	task_defaults.pmodel = pmodel;
 
 	qargs.task_defaults = task_defaults;
 	qargs.n_started = 0;
@@ -314,6 +323,10 @@ int main(int argc, char *argv[])
 	int noscale = 0;
 	Stream *st;
 	Crystal **crystals;
+	char *pmodel_str = NULL;
+	PartialityModel pmodel = PMODEL_SPHERE;
+	int min_measurements = 2;
+	char *rval;
 
 	/* Long options */
 	const struct option longopts[] = {
@@ -326,6 +339,9 @@ int main(int argc, char *argv[])
 		{"iterations",         1, NULL,               'n'},
 		{"no-scale",           0, &noscale,            1},
 		{"reference",          1, NULL,               'r'},
+		{"model",              1, NULL,               'm'},
+		{"min-measurements",   1, NULL,                2},
+
 		{0, 0, NULL, 0}
 	};
 
@@ -336,7 +352,7 @@ int main(int argc, char *argv[])
 	}
 
 	/* Short options */
-	while ((c = getopt_long(argc, argv, "hi:o:g:b:y:n:r:j:",
+	while ((c = getopt_long(argc, argv, "hi:o:g:b:y:n:r:j:m:",
 	                        longopts, NULL)) != -1)
 	{
 
@@ -370,6 +386,10 @@ int main(int argc, char *argv[])
 			n_iter = atoi(optarg);
 			break;
 
+			case 'm' :
+			pmodel_str = strdup(optarg);
+			break;
+
 			case 'b' :
 			beam = get_beam_parameters(optarg);
 			if ( beam == NULL ) {
@@ -383,6 +403,15 @@ int main(int argc, char *argv[])
 			reference_file = strdup(optarg);
 			break;
 
+			case 2 :
+			errno = 0;
+			min_measurements = strtod(optarg, &rval);
+			if ( *rval != '\0' ) {
+				ERROR("Invalid value for --min-measurements.\n");
+				return 1;
+			}
+			break;
+
 			case 0 :
 			break;
 
@@ -434,6 +463,17 @@ int main(int argc, char *argv[])
 		return 1;
 	}
 
+	if ( pmodel_str != NULL ) {
+		if ( strcmp(pmodel_str, "sphere") == 0 ) {
+			pmodel = PMODEL_SPHERE;
+		} else if ( strcmp(pmodel_str, "unity") == 0 ) {
+			pmodel = PMODEL_UNITY;
+		} else {
+			ERROR("Unknown partiality model '%s'.\n", pmodel_str);
+			return 1;
+		}
+	}
+
 	if ( reference_file != NULL ) {
 
 		RefList *list;
@@ -528,6 +568,7 @@ int main(int argc, char *argv[])
 		display_progress(n_images, n_crystals);
 
 	} while ( 1 );
+	fprintf(stderr, "\n");
 
 	close_stream(st);
 
@@ -544,18 +585,19 @@ int main(int argc, char *argv[])
 			crystal_set_image(cryst, &images[i]);
 
 			/* Now it's safe to do the following */
-			update_partialities(cryst);
+			update_partialities(cryst, pmodel);
 			as = crystal_get_reflections(cryst);
 			nobs += select_scalable_reflections(as, reference);
 
 		}
 	}
+	STATUS("%i scalable observations.\n", nobs);
 
 	/* Make initial estimates */
-	STATUS("\nPerforming initial scaling.\n");
+	STATUS("Performing initial scaling.\n");
 	if ( noscale ) STATUS("Scale factors fixed at 1.\n");
 	full = scale_intensities(crystals, n_crystals, reference,
-	                         nthreads, noscale);
+	                         nthreads, noscale, pmodel, min_measurements);
 
 	sr = sr_titlepage(crystals, n_crystals, "scaling-report.pdf",
 	                  infile, cmdline);
@@ -578,7 +620,7 @@ int main(int argc, char *argv[])
 		/* Refine the geometry of all patterns to get the best fit */
 		select_reflections_for_refinement(crystals, n_crystals,
 		                                  comp, have_reference);
-		refine_all(crystals, n_crystals, det, comp, nthreads);
+		refine_all(crystals, n_crystals, det, comp, nthreads, pmodel);
 
 		nobs = 0;
 		for ( j=0; j<n_crystals; j++ ) {
@@ -593,7 +635,8 @@ int main(int argc, char *argv[])
 		/* Re-estimate all the full intensities */
 		reflist_free(full);
 		full = scale_intensities(crystals, n_crystals,
-		                         reference, nthreads, noscale);
+		                         reference, nthreads, noscale, pmodel,
+		                         min_measurements);
 
 		sr_iteration(sr, i+1, crystals, n_crystals, full);
 
diff --git a/src/post-refinement.c b/src/post-refinement.c
index 1439b148..9e4649a2 100644
--- a/src/post-refinement.c
+++ b/src/post-refinement.c
@@ -81,27 +81,27 @@ static double partiality_rgradient(double r, double profile_radius)
 	dpdq = 6.0*(q-pow(q, 2.0));
 
 	/* dq/drad */
-	dqdrad = 0.5 * (1.0 - r * pow(profile_radius, -2.0));
+	dqdrad = -0.5 * r * pow(profile_radius, -2.0);
 
 	return dpdq * dqdrad;
 }
 
 
 /* Return the gradient of parameter 'k' given the current status of 'image'. */
-double gradient(Crystal *cr, int k, Reflection *refl)
+double gradient(Crystal *cr, int k, Reflection *refl, PartialityModel pmodel)
 {
-	double ds, azix, aziy;
-	double ttlow, tthigh, tt;
-	double nom, den;
-	double g;
+	double ds, azi;
+	double glow, ghigh;
 	double asx, asy, asz;
 	double bsx, bsy, bsz;
 	double csx, csy, csz;
 	double xl, yl, zl;
 	signed int hs, ks, ls;
-	double r1, r2, p;
+	double rlow, rhigh, p;
 	int clamp_low, clamp_high;
-	double klow, khigh;
+	double philow, phihigh, phi;
+	double khigh, klow;
+	double tl, cet, cez;
 	double gr;
 	struct image *image = crystal_get_image(cr);
 	double r = crystal_get_profile_radius(cr);
@@ -116,33 +116,40 @@ double gradient(Crystal *cr, int k, Reflection *refl)
 	zl = hs*asz + ks*bsz + ls*csz;
 
 	ds = 2.0 * resolution(crystal_get_cell(cr), hs, ks, ls);
-	get_partial(refl, &r1, &r2, &p, &clamp_low, &clamp_high);
+	get_partial(refl, &rlow, &rhigh, &p, &clamp_low, &clamp_high);
 
+	/* "low" gives the largest Ewald sphere (wavelength short => k large)
+	 * "high" gives the smallest Ewald sphere (wavelength long => k small)
+	 */
 	klow = 1.0/(image->lambda - image->lambda*image->bw/2.0);
 	khigh = 1.0/(image->lambda + image->lambda*image->bw/2.0);
-	ttlow = angle_between(0.0, 0.0, 1.0, xl, yl, zl+klow);
-	tthigh = angle_between(0.0, 0.0, 1.0, xl, yl, zl+khigh);
-	if ( (clamp_low == 0) && (clamp_high == 0) ) {
-		tt = (ttlow+tthigh)/2.0;
-	} else if ( clamp_high == 0 ) {
-		tt = tthigh + image->div;
-	} else if ( clamp_low == 0 ) {
-		tt = ttlow - image->div;
-	} else {
-		tt = 0.0;
-		/* Gradient should come out as zero in this case */
-	}
 
-	azix = angle_between(1.0, 0.0, 0.0, xl, yl, 0.0);
-	aziy = angle_between(0.0, 1.0, 0.0, xl, yl, 0.0);
+	tl = sqrt(xl*xl + yl*yl);
+	ds = modulus(xl, yl, zl);
+
+	cet = -sin(image->div/2.0) * klow;
+	cez = -cos(image->div/2.0) * klow;
+	philow = M_PI_2 - angle_between_2d(tl-cet, zl-cez, 1.0, 0.0);
+
+	cet = -sin(image->div/2.0) * khigh;
+	cez = -cos(image->div/2.0) * khigh;
+	phihigh = M_PI_2 - angle_between_2d(tl-cet, zl-cez, 1.0, 0.0);
+
+	/* Approximation: philow and phihigh are very similar */
+	phi = (philow + phihigh) / 2.0;
+
+	azi = atan2(yl, xl);
 
 	/* Calculate the gradient of partiality wrt excitation error. */
-	g = 0.0;
 	if ( clamp_low == 0 ) {
-		g -= partiality_gradient(r1, r);
+		glow = partiality_gradient(rlow, r);
+	} else {
+		glow = 0.0;
 	}
 	if ( clamp_high == 0 ) {
-		g += partiality_gradient(r2, r);
+		ghigh = partiality_gradient(rhigh, r);
+	} else {
+		ghigh = 0.0;
 	}
 
 	/* For many gradients, just multiply the above number by the gradient
@@ -150,57 +157,41 @@ double gradient(Crystal *cr, int k, Reflection *refl)
 	switch ( k ) {
 
 		case REF_DIV :
-		gr = 0.0;
-		if ( clamp_low == 0 ) {
-			nom = sqrt(2.0) * ds * sin(image->div/2.0);
-			den = sqrt(1.0 - cos(image->div/2.0));
-			gr -= (nom/den) * g;
-		}
-		if ( clamp_high == 0 ) {
-			nom = sqrt(2.0) * ds * sin(image->div/2.0);
-			den = sqrt(1.0 - cos(image->div/2.0));
-			gr += (nom/den) * g;
-		}
-		if ( isnan(gr) ) gr = 0.0;  /* FIXME: This isn't true (?) */
-		return gr / 4.0;  /* FIXME: Shameless fudge factor */
+		/* Small angle approximation */
+		return (ds*glow + ds*ghigh) / 2.0;
 
 		case REF_R :
-		g = 0.0;
-		if ( clamp_low == 0 ) {
-			g += partiality_rgradient(r1, r);
-		}
-		if ( clamp_high == 0 ) {
-			g += partiality_rgradient(r2, r);
-		}
-		return g;
+		gr  = partiality_rgradient(rlow, r);
+		gr -= partiality_rgradient(rhigh, r);
+		return gr;
 
 		/* Cell parameters and orientation */
 		case REF_ASX :
-		return hs * sin(tt) * cos(azix) * g;
+		return hs * sin(phi) * cos(azi) * (ghigh-glow);
 
 		case REF_BSX :
-		return ks * sin(tt) * cos(azix) * g;
+		return ks * sin(phi) * cos(azi) * (ghigh-glow);
 
 		case REF_CSX :
-		return ls * sin(tt) * cos(azix) * g;
+		return ls * sin(phi) * cos(azi) * (ghigh-glow);
 
 		case REF_ASY :
-		return hs * sin(tt) * cos(aziy) * g;
+		return hs * sin(phi) * sin(azi) * (ghigh-glow);
 
 		case REF_BSY :
-		return ks * sin(tt) * cos(aziy) * g;
+		return ks * sin(phi) * sin(azi) * (ghigh-glow);
 
 		case REF_CSY :
-		return ls * sin(tt) * cos(aziy) * g;
+		return ls * sin(phi) * sin(azi) * (ghigh-glow);
 
 		case REF_ASZ :
-		return hs * cos(tt) * g;
+		return hs * cos(phi) * (ghigh-glow);
 
 		case REF_BSZ :
-		return ks * cos(tt) * g;
+		return ks * cos(phi) * (ghigh-glow);
 
 		case REF_CSZ :
-		return ls * cos(tt) * g;
+		return ls * cos(phi) * (ghigh-glow);
 
 	}
 
@@ -364,7 +355,8 @@ static gsl_vector *solve_svd(gsl_vector *v, gsl_matrix *M)
 
 
 /* Perform one cycle of post refinement on 'image' against 'full' */
-static double pr_iterate(Crystal *cr, const RefList *full)
+static double pr_iterate(Crystal *cr, const RefList *full,
+                         PartialityModel pmodel)
 {
 	gsl_matrix *M;
 	gsl_vector *v;
@@ -420,7 +412,7 @@ static double pr_iterate(Crystal *cr, const RefList *full)
 		/* Calculate all gradients for this reflection */
 		for ( k=0; k<NUM_PARAMS; k++ ) {
 			double gr;
-			gr = gradient(cr, k, refl);
+			gr = gradient(cr, k, refl, pmodel);
 			gradients[k] = gr;
 		}
 
@@ -532,7 +524,7 @@ static double guide_dev(Crystal *cr, const RefList *full)
 }
 
 
-void pr_refine(Crystal *cr, const RefList *full)
+void pr_refine(Crystal *cr, const RefList *full, PartialityModel pmodel)
 {
 	double max_shift, dev;
 	int i;
@@ -557,9 +549,9 @@ void pr_refine(Crystal *cr, const RefList *full)
 		cell_get_reciprocal(crystal_get_cell(cr), &asx, &asy, &asz,
 			               &bsx, &bsy, &bsz, &csx, &csy, &csz);
 
-		max_shift = pr_iterate(cr, full);
+		max_shift = pr_iterate(cr, full, pmodel);
 
-		update_partialities(cr);
+		update_partialities(cr, pmodel);
 
 		if ( verbose ) {
 			dev = guide_dev(cr, full);
diff --git a/src/post-refinement.h b/src/post-refinement.h
index fe171882..7b822938 100644
--- a/src/post-refinement.h
+++ b/src/post-refinement.h
@@ -40,6 +40,7 @@
 #include "image.h"
 #include "utils.h"
 #include "crystal.h"
+#include "geometry.h"
 
 
 /* Refineable parameters */
@@ -59,10 +60,11 @@ enum {
 };
 
 
-extern void pr_refine(Crystal *cr, const RefList *full);
+extern void pr_refine(Crystal *cr, const RefList *full, PartialityModel pmodel);
 
 /* Exported so it can be poked by tests/pr_gradient_check */
-extern double gradient(Crystal *cr, int k, Reflection *refl);
+extern double gradient(Crystal *cr, int k, Reflection *refl,
+                       PartialityModel pmodel);
 
 
 #endif	/* POST_REFINEMENT_H */
diff --git a/tests/partialator_merge_check_1 b/tests/partialator_merge_check_1
new file mode 100755
index 00000000..9125b75b
--- /dev/null
+++ b/tests/partialator_merge_check_1
@@ -0,0 +1,78 @@
+#!/bin/sh
+
+cat > partialator_merge_check_1.stream << EOF
+CrystFEL stream format 2.1
+Command line: indexamajig -i dummy.lst -o dummy.stream --kraken=prawn
+----- Begin chunk -----
+Image filename: dummy.h5
+photon_energy_eV = 2000.0
+--- Begin crystal
+Cell parameters 27.74398 27.84377 16.90346 nm, 88.53688 91.11774 118.75944 deg
+astar = -0.0283891 +0.0149254 -0.0257273 nm^-1
+bstar = -0.0068281 +0.0403989 -0.0005196 nm^-1
+cstar = +0.0406926 +0.0052233 -0.0426520 nm^-1
+Reflections measured after indexing
+  h   k   l          I    phase   sigma(I)  counts  fs/px  ss/px
+  1   0   0     100.00        -       1.00       1  938.0  629.0
+End of reflections
+--- End crystal
+--- Begin crystal
+Cell parameters 27.74398 27.84377 16.90346 nm, 88.53688 91.11774 118.75944 deg
+astar = -0.0283891 +0.0149254 -0.0257273 nm^-1
+bstar = -0.0068281 +0.0403989 -0.0005196 nm^-1
+cstar = +0.0406926 +0.0052233 -0.0426520 nm^-1
+Reflections measured after indexing
+  h   k   l          I    phase   sigma(I)  counts  fs/px  ss/px
+  1   0   0     200.00        -       1.00       1  938.0  629.0
+End of reflections
+--- End crystal
+----- End chunk -----
+EOF
+
+# We merge two patterns, without scaling or partiality, the result should just
+# be an average.
+cat > partialator_merge_check_1_ans.hkl << EOF
+  h   k   l          I    phase   sigma(I)  counts  fs/px  ss/px
+  1   0   0     150.00        -      35.36       2    0.0    0.0
+End of reflections
+EOF
+
+cat > partialator_merge_check_1.beam << EOF
+beam/fluence = 2.0e15
+beam/radius = 1.5e-6
+beam/photon_energy = 6000.0
+beam/bandwidth = 0.0005
+beam/divergence = 0.001
+profile_radius = 0.005e9
+EOF
+
+cat > partialator_merge_check_1.geom << EOF
+0/min_fs = 0
+0/max_fs = 1023
+0/min_ss = 0
+0/max_ss = 1023
+0/corner_x = -512.00
+0/corner_y = -512.00
+0/fs = x
+0/ss = y
+0/clen = 70.0e-3
+0/res = 13333.3  ; 75 micron pixel size
+0/badrow_direction = y
+0/adu_per_eV = 1.0
+EOF
+
+
+src/partialator -i partialator_merge_check_1.stream \
+                -o partialator_merge_check_1.hkl  \
+                -g partialator_merge_check_1.geom \
+                -b partialator_merge_check_1.beam \
+                --model=unity --iterations=0 --no-scale
+
+diff partialator_merge_check_1.hkl partialator_merge_check_1_ans.hkl
+if [ $? -ne 0 ]; then
+	exit 1
+fi
+rm -f partialator_merge_check_1.stream partialator_merge_check_1.hkl \
+      partialator_merge_check_1_ans.hkl partialator_merge_check_1.beam \
+      partialator_merge_check_1.geom
+exit 0
diff --git a/tests/partialator_merge_check_2 b/tests/partialator_merge_check_2
new file mode 100755
index 00000000..50b0c9bf
--- /dev/null
+++ b/tests/partialator_merge_check_2
@@ -0,0 +1,79 @@
+#!/bin/sh
+
+cat > partialator_merge_check_2.stream << EOF
+CrystFEL stream format 2.1
+Command line: indexamajig -i dummy.lst -o dummy.stream --kraken=prawn
+----- Begin chunk -----
+Image filename: dummy.h5
+photon_energy_eV = 2000.0
+--- Begin crystal
+Cell parameters 27.74398 27.84377 16.90346 nm, 88.53688 91.11774 118.75944 deg
+astar = -0.0283891 +0.0149254 -0.0257273 nm^-1
+bstar = -0.0068281 +0.0403989 -0.0005196 nm^-1
+cstar = +0.0406926 +0.0052233 -0.0426520 nm^-1
+Reflections measured after indexing
+  h   k   l          I    phase   sigma(I)  counts  fs/px  ss/px
+  1   0   0     100.00        -       1.00       1  938.0  629.0
+End of reflections
+--- End crystal
+--- Begin crystal
+Cell parameters 27.74398 27.84377 16.90346 nm, 88.53688 91.11774 118.75944 deg
+astar = -0.0283891 +0.0149254 -0.0257273 nm^-1
+bstar = -0.0068281 +0.0403989 -0.0005196 nm^-1
+cstar = +0.0406926 +0.0052233 -0.0426520 nm^-1
+Reflections measured after indexing
+  h   k   l          I    phase   sigma(I)  counts  fs/px  ss/px
+  1   0   0     200.00        -       1.00       1  938.0  629.0
+End of reflections
+--- End crystal
+----- End chunk -----
+EOF
+
+# We merge two patterns, without partiality but with scaling, the result should
+# be the mean but with the standard deviation should be zero because the scaling
+# factor can absorb the difference.
+cat > partialator_merge_check_2_ans.hkl << EOF
+  h   k   l          I    phase   sigma(I)  counts  fs/px  ss/px
+  1   0   0     150.00        -       0.00       2    0.0    0.0
+End of reflections
+EOF
+
+cat > partialator_merge_check_2.beam << EOF
+beam/fluence = 2.0e15
+beam/radius = 1.5e-6
+beam/photon_energy = 6000.0
+beam/bandwidth = 0.0005
+beam/divergence = 0.001
+profile_radius = 0.005e9
+EOF
+
+cat > partialator_merge_check_2.geom << EOF
+0/min_fs = 0
+0/max_fs = 1023
+0/min_ss = 0
+0/max_ss = 1023
+0/corner_x = -512.00
+0/corner_y = -512.00
+0/fs = x
+0/ss = y
+0/clen = 70.0e-3
+0/res = 13333.3  ; 75 micron pixel size
+0/badrow_direction = y
+0/adu_per_eV = 1.0
+EOF
+
+
+src/partialator -i partialator_merge_check_2.stream \
+                -o partialator_merge_check_2.hkl  \
+                -g partialator_merge_check_2.geom \
+                -b partialator_merge_check_2.beam \
+                --model=unity --iterations=0
+
+diff partialator_merge_check_2.hkl partialator_merge_check_2_ans.hkl
+if [ $? -ne 0 ]; then
+	exit 1
+fi
+rm -f partialator_merge_check_2.stream partialator_merge_check_2.hkl \
+      partialator_merge_check_2_ans.hkl partialator_merge_check_2.beam \
+      partialator_merge_check_2.geom
+exit 0
diff --git a/tests/partialator_merge_check_3 b/tests/partialator_merge_check_3
new file mode 100755
index 00000000..82b3b787
--- /dev/null
+++ b/tests/partialator_merge_check_3
@@ -0,0 +1,81 @@
+#!/bin/sh
+
+cat > partialator_merge_check_3.stream << EOF
+CrystFEL stream format 2.1
+Command line: indexamajig -i dummy.lst -o dummy.stream --kraken=prawn
+----- Begin chunk -----
+Image filename: dummy.h5
+photon_energy_eV = 2000.0
+--- Begin crystal
+Cell parameters 27.74398 27.84377 16.90346 nm, 88.53688 91.11774 118.75944 deg
+astar = -0.0283891 +0.0149254 -0.0257273 nm^-1
+bstar = -0.0068281 +0.0403989 -0.0005196 nm^-1
+cstar = +0.0406926 +0.0052233 -0.0426520 nm^-1
+Reflections measured after indexing
+  h   k   l          I    phase   sigma(I)  counts  fs/px  ss/px
+  0   1   0     110.00        -       1.00       1  938.0  629.0
+  2   0   0     190.00        -       1.00       1  938.0  629.0
+End of reflections
+--- End crystal
+--- Begin crystal
+Cell parameters 27.74398 27.84377 16.90346 nm, 88.53688 91.11774 118.75944 deg
+astar = -0.0283891 +0.0149254 -0.0257273 nm^-1
+bstar = -0.0068281 +0.0403989 -0.0005196 nm^-1
+cstar = +0.0406926 +0.0052233 -0.0426520 nm^-1
+Reflections measured after indexing
+  h   k   l          I    phase   sigma(I)  counts  fs/px  ss/px
+  1   0   0      50.00        -       1.00       1  938.0  629.0
+  0   2   0     100.00        -       1.00       1  938.0  629.0
+  0   1   1       1.00        -       1.00       1  938.0  629.0
+End of reflections
+--- End crystal
+----- End chunk -----
+EOF
+
+# W
+cat > partialator_merge_check_3_ans.hkl << EOF
+  h   k   l          I    phase   sigma(I)  counts  fs/px  ss/px
+  1   0   0      81.17        -       3.62       2    0.0    0.0
+  2   0   0     144.36        -       2.04       2    0.0    0.0
+End of reflections
+EOF
+
+cat > partialator_merge_check_3.beam << EOF
+beam/fluence = 2.0e15
+beam/radius = 1.5e-6
+beam/photon_energy = 6000.0
+beam/bandwidth = 0.0005
+beam/divergence = 0.001
+profile_radius = 0.005e9
+EOF
+
+cat > partialator_merge_check_3.geom << EOF
+0/min_fs = 0
+0/max_fs = 1023
+0/min_ss = 0
+0/max_ss = 1023
+0/corner_x = -512.00
+0/corner_y = -512.00
+0/fs = x
+0/ss = y
+0/clen = 70.0e-3
+0/res = 13333.3  ; 75 micron pixel size
+0/badrow_direction = y
+0/adu_per_eV = 1.0
+EOF
+
+
+src/partialator -i partialator_merge_check_3.stream \
+                -o partialator_merge_check_3.hkl  \
+                -g partialator_merge_check_3.geom \
+                -b partialator_merge_check_3.beam \
+                --model=unity --iterations=0 -y 4
+
+diff partialator_merge_check_3.hkl partialator_merge_check_3_ans.hkl
+if [ $? -ne 0 ]; then
+	exit 1
+fi
+rm -f partialator_merge_check_3.stream partialator_merge_check_3.hkl \
+      partialator_merge_check_3_ans.hkl partialator_merge_check_3.beam \
+      partialator_merge_check_3.geom
+exit 0
diff --git a/tests/plot_gradients b/tests/plot_gradients
new file mode 100755
index 00000000..b5800d07
--- /dev/null
+++ b/tests/plot_gradients
@@ -0,0 +1,28 @@
+#!/bin/sh
+
+gnuplot -persist << EOF
+
+set key bottom right
+
+set xlabel "Calculated gradient"
+set ylabel "Observed gradient"
+
+  plot "gradient-test-x.dat" using 1:2 w p lc 1 pt 1 title "x"
+replot "gradient-test-y.dat" using 1:2 w p lc 2 pt 1 title "y"
+replot "gradient-test-z.dat" using 1:2 w p lc 3 pt 1 title "z"
+
+EOF
+
+gnuplot -persist << EOF
+set key bottom right
+set xlabel "Calculated gradient"
+set ylabel "Observed gradient"
+plot "gradient-test-R.dat" using 1:2 w p lc 1 pt 1 title "profile radius"
+EOF
+
+gnuplot -persist << EOF
+set key bottom right
+set xlabel "Calculated gradient"
+set ylabel "Observed gradient"
+plot "gradient-test-div.dat" using 1:2 w p lc 1 pt 1 title "divergence"
+EOF
diff --git a/tests/pr_gradient_check.c b/tests/pr_gradient_check.c
index 91b51536..68a625a7 100644
--- a/tests/pr_gradient_check.c
+++ b/tests/pr_gradient_check.c
@@ -29,6 +29,8 @@
 
 #include <stdlib.h>
 #include <stdio.h>
+#include <gsl/gsl_statistics.h>
+#include <getopt.h>
 
 #include <image.h>
 #include <cell.h>
@@ -155,6 +157,7 @@ static Crystal *new_shifted_crystal(Crystal *cr, int refine, double incr_val)
 
 		case REF_R :
 		cell = cell_new_from_cell(crystal_get_cell(cr));
+		crystal_set_cell(cr_new, cell);
 		crystal_set_profile_radius(cr_new, r + incr_val);
 		break;
 
@@ -218,8 +221,9 @@ static void calc_either_side(Crystal *cr, double incr_val,
 
 
 
-static int test_gradients(Crystal *cr, double incr_val, int refine,
-                          const char *str)
+static double test_gradients(Crystal *cr, double incr_val, int refine,
+                             const char *str, const char *file,
+                             PartialityModel pmodel, int quiet, int plot)
 {
 	Reflection *refl;
 	RefListIterator *iter;
@@ -227,9 +231,16 @@ static int test_gradients(Crystal *cr, double incr_val, int refine,
 	int i;
 	int *valid;
 	int nref;
-	int n_acc, n_valid;
+	int n_good, n_invalid, n_small, n_nan, n_bad;
 	RefList *reflections;
-	//FILE *fh;
+	FILE *fh;
+	int ntot = 0;
+	double total = 0.0;
+	char tmp[32];
+	double *vec1;
+	double *vec2;
+	int n_line;
+	double cc;
 
 	reflections = find_intersections(crystal_get_image(cr), cr);
 	crystal_set_reflections(cr, reflections);
@@ -237,7 +248,7 @@ static int test_gradients(Crystal *cr, double incr_val, int refine,
 	nref = num_reflections(reflections);
 	if ( nref < 10 ) {
 		ERROR("Too few reflections found.  Failing test by default.\n");
-		return -1;
+		return 0.0;
 	}
 
 	vals[0] = malloc(nref*sizeof(long double));
@@ -245,13 +256,13 @@ static int test_gradients(Crystal *cr, double incr_val, int refine,
 	vals[2] = malloc(nref*sizeof(long double));
 	if ( (vals[0] == NULL) || (vals[1] == NULL) || (vals[2] == NULL) ) {
 		ERROR("Couldn't allocate memory.\n");
-		return -1;
+		return 0.0;
 	}
 
 	valid = malloc(nref*sizeof(int));
 	if ( valid == NULL ) {
 		ERROR("Couldn't allocate memory.\n");
-		return -1;
+		return 0.0;
 	}
 	for ( i=0; i<nref; i++ ) valid[i] = 1;
 
@@ -259,9 +270,20 @@ static int test_gradients(Crystal *cr, double incr_val, int refine,
 
 	calc_either_side(cr, incr_val, valid, vals, refine);
 
-	//fh = fopen("wrongness.dat", "a");
+	if ( plot ) {
+		snprintf(tmp, 32, "gradient-test-%s.dat", file);
+		fh = fopen(tmp, "w");
+	}
+
+	vec1 = malloc(nref*sizeof(double));
+	vec2 = malloc(nref*sizeof(double));
+	if ( (vec1 == NULL) || (vec2 == NULL) ) {
+		ERROR("Couldn't allocate memory.\n");
+		return 0.0;
+	}
 
-	n_valid = nref;  n_acc = 0;
+	n_invalid = 0;  n_good = 0;
+	n_nan = 0;  n_small = 0;  n_bad = 0;  n_line = 0;
 	i = 0;
 	for ( refl = first_refl(reflections, &iter);
 	      refl != NULL;
@@ -275,7 +297,8 @@ static int test_gradients(Crystal *cr, double incr_val, int refine,
 		get_indices(refl, &h, &k, &l);
 
 		if ( !valid[i] ) {
-			n_valid--;
+			n_invalid++;
+			i++;
 		} else {
 
 			double r1, r2, p;
@@ -284,20 +307,45 @@ static int test_gradients(Crystal *cr, double incr_val, int refine,
 			grad1 = (vals[1][i] - vals[0][i]) / incr_val;
 			grad2 = (vals[2][i] - vals[1][i]) / incr_val;
 			grad = (grad1 + grad2) / 2.0;
+			i++;
 
-			cgrad = gradient(cr, refine, refl);
+			cgrad = gradient(cr, refine, refl, pmodel);
 
 			get_partial(refl, &r1, &r2, &p, &cl, &ch);
 
-			if ( (fabs(cgrad) > 5e-8) &&
-			     !within_tolerance(grad, cgrad, 10.0) )
+			if ( isnan(cgrad) ) {
+				n_nan++;
+				continue;
+			}
+
+			if ( plot ) {
+				fprintf(fh, "%e %Le\n", cgrad, grad);
+			}
+
+			vec1[n_line] = cgrad;
+			vec2[n_line] = grad;
+			n_line++;
+
+			if ( (fabs(cgrad) < 5e-8) && (fabs(grad) < 5e-8) ) {
+				n_small++;
+				continue;
+			}
+
+			total += fabs(cgrad - grad);
+			ntot++;
+
+			if ( !within_tolerance(grad, cgrad, 5.0)
+			  || !within_tolerance(cgrad, grad, 5.0) )
 			{
 
-				STATUS("!- %s %3i %3i %3i"
-				       " %10.2Le %10.2e ratio = %5.2Lf"
-				       " %10.2e %10.2e\n",
-				       str, h, k, l, grad, cgrad, cgrad/grad,
-				       r1, r2);
+				if ( !quiet ) {
+					STATUS("!- %s %3i %3i %3i"
+					       " %10.2Le %10.2e ratio = %5.2Lf"
+					       " %10.2e %10.2e\n",
+					       str, h, k, l, grad, cgrad,
+					       cgrad/grad, r1, r2);
+				}
+				n_bad++;
 
 			} else {
 
@@ -307,29 +355,24 @@ static int test_gradients(Crystal *cr, double incr_val, int refine,
 				//       str, h, k, l, grad, cgrad, cgrad/grad,
 				//       r1, r2);
 
-				n_acc++;
+				n_good++;
 
 			}
 
-			//fprintf(fh, "%e %f\n",
-			        //resolution(image->indexed_cell, h, k, l),
-			        //rad2deg(tt),
-			//        cgrad,
-			//        fabs((grad-cgrad)/grad));
-
 		}
 
-		i++;
-
 	}
 
-	STATUS("%s: %i out of %i valid gradients were accurate.\n",
-	       str, n_acc, n_valid);
-	//fclose(fh);
+	STATUS("%3s: %3i within 5%%, %3i outside, %3i nan, %3i invalid, "
+	       "%3i small. ", str, n_good, n_bad, n_nan, n_invalid, n_small);
 
-	if ( n_acc != n_valid ) return 1;
+	if ( plot ) {
+		fclose(fh);
+	}
 
-	return 0;
+	cc = gsl_stats_correlation(vec1, 1, vec2, 1, n_line);
+	STATUS("CC = %+f\n", cc);
+	return cc;
 }
 
 
@@ -345,7 +388,34 @@ int main(int argc, char *argv[])
 	Crystal *cr;
 	struct quaternion orientation;
 	int i;
-	int val;
+	const PartialityModel pmodel = PMODEL_SPHERE;
+	int fail = 0;
+	int quiet = 0;
+	int plot = 0;
+	int c;
+
+	const struct option longopts[] = {
+		{"quiet",       0, &quiet,        1},
+		{"plot",        0, &plot,         1},
+		{0, 0, NULL, 0}
+	};
+
+	while ((c = getopt_long(argc, argv, "", longopts, NULL)) != -1) {
+		switch (c) {
+
+			case 0 :
+			break;
+
+			case '?' :
+			break;
+
+			default :
+			ERROR("Unhandled option '%c'\n", c);
+			break;
+
+		}
+
+	}
 
 	image.width = 1024;
 	image.height = 1024;
@@ -373,11 +443,10 @@ int main(int argc, char *argv[])
 	                                deg2rad(90.0),
 	                                deg2rad(90.0));
 
-	val = 0;
-
 	for ( i=0; i<1; i++ ) {
 
 		UnitCell *rot;
+		double val;
 
 		orientation = random_quaternion();
 		rot = cell_rotate(cell, orientation);
@@ -388,32 +457,55 @@ int main(int argc, char *argv[])
 			            &bz, &cx, &cy, &cz);
 
 		incr_val = incr_frac * image.div;
-		val += test_gradients(cr, incr_val, REF_DIV, "div");
+		val =  test_gradients(cr, incr_val, REF_DIV, "div", "div",
+		                      pmodel, quiet, plot);
+		if ( val < 0.99 ) fail = 1;
 
-		incr_val = incr_frac * ax;
-		val += test_gradients(cr, incr_val, REF_ASX, "ax*");
-		incr_val = incr_frac * ay;
-		val += test_gradients(cr, incr_val, REF_ASY, "ay*");
-		incr_val = incr_frac * az;
-		val += test_gradients(cr, incr_val, REF_ASZ, "az*");
+		incr_val = incr_frac * crystal_get_profile_radius(cr);
+		val = test_gradients(cr, incr_val, REF_R, "R", "R", pmodel,
+		                     quiet, plot);
+		if ( val < 0.99 ) fail = 1;
 
+		incr_val = incr_frac * ax;
+		val = test_gradients(cr, incr_val, REF_ASX, "ax*", "x", pmodel,
+		                     quiet, plot);
+		if ( val < 0.99 ) fail = 1;
 		incr_val = incr_frac * bx;
-		val += test_gradients(cr, incr_val, REF_BSX, "bx*");
-		incr_val = incr_frac * by;
-		val += test_gradients(cr, incr_val, REF_BSY, "by*");
-		incr_val = incr_frac * bz;
-		val += test_gradients(cr, incr_val, REF_BSZ, "bz*");
-
+		val = test_gradients(cr, incr_val, REF_BSX, "bx*", "x", pmodel,
+		                     quiet, plot);
+		if ( val < 0.99 ) fail = 1;
 		incr_val = incr_frac * cx;
-		val += test_gradients(cr, incr_val, REF_CSX, "cx*");
+		val = test_gradients(cr, incr_val, REF_CSX, "cx*", "x", pmodel,
+		                     quiet, plot);
+		if ( val < 0.99 ) fail = 1;
+
+		incr_val = incr_frac * ay;
+		val = test_gradients(cr, incr_val, REF_ASY, "ay*", "y", pmodel,
+		                     quiet, plot);
+		if ( val < 0.99 ) fail = 1;
+		incr_val = incr_frac * by;
+		val = test_gradients(cr, incr_val, REF_BSY, "by*", "y", pmodel,
+		                     quiet, plot);
+		if ( val < 0.99 ) fail = 1;
 		incr_val = incr_frac * cy;
-		val += test_gradients(cr, incr_val, REF_CSY, "cy*");
+		val = test_gradients(cr, incr_val, REF_CSY, "cy*", "y", pmodel,
+		                     quiet, plot);
+		if ( val < 0.99 ) fail = 1;
+
+		incr_val = incr_frac * az;
+		val = test_gradients(cr, incr_val, REF_ASZ, "az*", "z", pmodel,
+		                     quiet, plot);
+		if ( val < 0.99 ) fail = 1;
+		incr_val = incr_frac * bz;
+		val = test_gradients(cr, incr_val, REF_BSZ, "bz*", "z", pmodel,
+		                     quiet, plot);
+		if ( val < 0.99 ) fail = 1;
 		incr_val = incr_frac * cz;
-		val += test_gradients(cr, incr_val, REF_CSZ, "cz*");
+		val = test_gradients(cr, incr_val, REF_CSZ, "cz*", "z", pmodel,
+		                     quiet, plot);
+		if ( val < 0.99 ) fail = 1;
 
 	}
 
-	STATUS("Returning 0 by default: gradients only needed for experimental"
-	       " features of CrystFEL.\n");
-	return 0;
+	return fail;
 }