Merge branch 'tom/partials'

6 files changed, 293 insertions, 122 deletions

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