Merge branch 'tom/ginn-partials'

10 files changed, 1269 insertions, 633 deletions

diff --git a/src/merge.c b/src/merge.c
index 17e57c90..e60aca15 100644
--- a/src/merge.c
+++ b/src/merge.c
@@ -50,7 +50,7 @@
 
 
 /* Minimum partiality of a reflection for it to be merged */
-#define MIN_PART_MERGE (0.05)
+#define MIN_PART_MERGE (0.3)
 
 
 struct merge_queue_args
@@ -59,7 +59,6 @@ struct merge_queue_args
 	pthread_rwlock_t full_lock;
 	Crystal **crystals;
 	int n_started;
-	PartialityModel pmodel;
 	double push_res;
 	int use_weak;
 	long long int n_reflections;
@@ -159,7 +158,6 @@ static void run_merge_job(void *vwargs, int cookie)
 		signed int h, k, l;
 		double mean, sumweight, M2, temp, delta, R;
 		double corr, res, w;
-		//double esd;
 
 		if ( get_partiality(refl) < MIN_PART_MERGE ) continue;
 
@@ -189,7 +187,7 @@ static void run_merge_job(void *vwargs, int cookie)
 		}
 
 		/* Total (multiplicative) correction factor */
-		corr = exp(-G) * exp(B*res*res) * get_lorentz(refl)
+		corr = G * exp(B*res*res) * get_lorentz(refl)
 		        / get_partiality(refl);
 		if ( isnan(corr) ) {
 			ERROR("NaN corr:\n");
@@ -200,8 +198,8 @@ static void run_merge_job(void *vwargs, int cookie)
 			continue;
 		}
 
-		//esd = get_esd_intensity(refl) * corr;
-		w = 1.0;
+		/* Reflections count less the more they have to be scaled up */
+		w = 1.0 / corr;
 
 		/* Running mean and variance calculation */
 		temp = w + sumweight;
@@ -229,7 +227,7 @@ static void finalise_merge_job(void *vqargs, void *vwargs)
 
 
 RefList *merge_intensities(Crystal **crystals, int n, int n_threads,
-                           PartialityModel pmodel, int min_meas,
+                           int min_meas,
                            double push_res, int use_weak)
 {
 	RefList *full;
@@ -245,7 +243,6 @@ RefList *merge_intensities(Crystal **crystals, int n, int n_threads,
 	qargs.full = full;
 	qargs.n_started = 0;
 	qargs.crystals = crystals;
-	qargs.pmodel = pmodel;
 	qargs.push_res = push_res;
 	qargs.use_weak = use_weak;
 	qargs.n_reflections = 0;
diff --git a/src/merge.h b/src/merge.h
index 1911584f..56229e4c 100644
--- a/src/merge.h
+++ b/src/merge.h
@@ -40,7 +40,6 @@
 #include "geometry.h"
 
 extern RefList *merge_intensities(Crystal **crystals, int n, int n_threads,
-                                  PartialityModel pmodel, int min_meas,
-                                  double push_res, int use_weak);
+                                  int min_meas, double push_res, int use_weak);
 
 #endif	/* MERGE */
diff --git a/src/partial_sim.c b/src/partial_sim.c
index 4a348742..fc7845c1 100644
--- a/src/partial_sim.c
+++ b/src/partial_sim.c
@@ -386,7 +386,7 @@ static void run_job(void *vwargs, int cookie)
 
 	reflections = predict_to_res(cr, largest_q(&wargs->image));
 	crystal_set_reflections(cr, reflections);
-	calculate_partialities(cr, PMODEL_SCSPHERE);
+	calculate_partialities(cr, PMODEL_XSPHERE);
 
 	for ( i=0; i<NBINS; i++ ) {
 		wargs->n_ref[i] = 0;
@@ -862,7 +862,7 @@ int main(int argc, char *argv[])
 		STATUS("                     Background: none (no image "
 		       "output)\n");
 	}
-	STATUS("               Partiality model: scsphere (hardcoded)\n");
+	STATUS("               Partiality model: xsphere (hardcoded)\n");
 	STATUS("       Noise standard deviation: %.2f detector units\n",
 	       noise_stddev);
 	if ( random_intensities ) {
@@ -874,6 +874,7 @@ int main(int argc, char *argv[])
 		       input_file, sym_str);
 	}
 	STATUS("   Max error in cell components: %.2f %%\n", cnoise);
+	STATUS("Scale factor standard deviation: %.2f\n", osf_stddev);
 
 	if ( random_intensities ) {
 		full = reflist_new();
diff --git a/src/partialator.c b/src/partialator.c
index 428b0d68..21b02266 100644
--- a/src/partialator.c
+++ b/src/partialator.c
@@ -3,11 +3,11 @@
  *
  * Scaling and post refinement for coherent nanocrystallography
  *
- * Copyright © 2012-2015 Deutsches Elektronen-Synchrotron DESY,
+ * Copyright © 2012-2018 Deutsches Elektronen-Synchrotron DESY,
  *                       a research centre of the Helmholtz Association.
  *
  * Authors:
- *   2010-2015 Thomas White <taw@physics.org>
+ *   2010-2018 Thomas White <taw@physics.org>
  *
  * This file is part of CrystFEL.
  *
@@ -40,6 +40,7 @@
 #include <assert.h>
 #include <pthread.h>
 #include <gsl/gsl_errno.h>
+#include <sys/stat.h>
 
 #include <image.h>
 #include <utils.h>
@@ -190,7 +191,7 @@ static void write_split(Crystal **crystals, int n_crystals, const char *outfile,
 	}
 	snprintf(tmp, 1024, "%s1", outfile);
 	split = merge_intensities(crystals1, n_crystals1, nthreads,
-		                  pmodel, min_measurements, push_res, 1);
+		                  min_measurements, push_res, 1);
 
 	if ( split == NULL ) {
 		ERROR("Not enough crystals for two way split!\n");
@@ -202,7 +203,7 @@ static void write_split(Crystal **crystals, int n_crystals, const char *outfile,
 	reflist_free(split);
 	snprintf(tmp, 1024, "%s2", outfile);
 	split = merge_intensities(crystals2, n_crystals2, nthreads,
-		                  pmodel, min_measurements, push_res, 1);
+		                  min_measurements, push_res, 1);
 	STATUS("and %s\n", tmp);
 	write_reflist_2(tmp, split, sym);
 	reflist_free(split);
@@ -289,7 +290,7 @@ static void write_custom_split(struct custom_split *csplit, int dsn,
 	STATUS("Writing dataset '%s' to %s (%i crystals)\n",
 	       csplit->dataset_names[dsn], tmp, n_crystalsn);
 	split = merge_intensities(crystalsn, n_crystalsn, nthreads,
-		                  pmodel, min_measurements, push_res, 1);
+		                  min_measurements, push_res, 1);
 	write_reflist_2(tmp, split, sym);
 	reflist_free(split);
 
@@ -325,7 +326,12 @@ static void show_help(const char *s)
 "  -j <n>                     Run <n> analyses in parallel.\n"
 "      --no-free              Disable cross-validation (testing only).\n"
 "      --custom-split         List of files for custom dataset splitting.\n"
-"      --max-rel-B            Maximum allowable relative |B| factor.\n");
+"      --max-rel-B            Maximum allowable relative |B| factor.\n"
+"      --no-logs              Do not write extensive log files.\n"
+"  -w <pg>                    Apparent point group for resolving ambiguities.\n"
+"      --operator=<op>        Indexing ambiguity operator for resolving.\n"
+"      --force-bandwidth=<n>  Set all bandwidths to <n> (fraction).\n"
+"      --force-radius=<n>     Set all profile radii to <n> nm^-1.\n");
 }
 
 
@@ -579,7 +585,7 @@ static int set_initial_params(Crystal *cr, FILE *fh)
 
 	} else {
 
-		crystal_set_osf(cr, 0.0);
+		crystal_set_osf(cr, 1.0);
 		crystal_set_Bfac(cr, 0.0);
 
 	}
@@ -604,13 +610,22 @@ static void select_free_reflections(RefList *list, gsl_rng *rng)
 
 
 static void write_to_pgraph(FILE *fh, RefList *list, RefList *full, Crystal *cr,
-                            int fr)
+                            int fr, signed int inum)
 {
 	Reflection *refl;
 	RefListIterator *iter;
 	double G = crystal_get_osf(cr);
 	double B = crystal_get_Bfac(cr);
 	UnitCell *cell = crystal_get_cell(cr);
+	char ins[5];
+
+	if ( inum >= 0 ) {
+		snprintf(ins, 4, "%i", inum);
+	} else {
+		ins[0] = 'F';
+		ins[1] = '\0';
+	}
+
 
 	for ( refl = first_refl(list, &iter);
 	      refl != NULL;
@@ -623,6 +638,9 @@ static void write_to_pgraph(FILE *fh, RefList *list, RefList *full, Crystal *cr,
 
 		if ( !get_flag(refl) ) continue;  /* Not free-flagged */
 
+		/* Strong reflections only */
+		if ( get_intensity(refl) < 3.0*get_esd_intensity(refl) ) continue;
+
 		get_indices(refl, &h, &k, &l);
 		res = resolution(cell, h, k, l);
 		if ( 2.0*res > crystal_get_resolution_limit(cr) ) continue;
@@ -630,46 +648,55 @@ static void write_to_pgraph(FILE *fh, RefList *list, RefList *full, Crystal *cr,
 		match = find_refl(full, h, k, l);
 		if ( match == NULL ) continue;
 
+		/* Don't calculate pobs if reference reflection is weak */
+		if ( fabs(get_intensity(match)) / get_esd_intensity(match) < 3.0 ) continue;
+
 		/* Calculated partiality */
 		pcalc = get_partiality(refl);
 
 		/* Observed partiality */
-		corr = exp(-G) * exp(B*res*res) * get_lorentz(refl);
+		corr = G * exp(B*res*res) * get_lorentz(refl);
 		Ipart = get_intensity(refl) * corr;
 		pobs = Ipart / get_intensity(match);
 
-		fprintf(fh, "%5i %4i %4i %4i %8.4f %8.3f %8.3f\n",
-		        fr, h, k, l, 2*res/1e9, pcalc, pobs);
+		fprintf(fh, "%5i %4i %4i %4i %e %e %8.3f %8.3f %s\n",
+		        fr, h, k, l, 2*res, Ipart, pcalc, pobs, ins);
 
 	}
 }
 
 
 static void write_pgraph(RefList *full, Crystal **crystals, int n_crystals,
-                         int iter)
+                         signed int iter, const char *suff)
 {
 	FILE *fh;
 	char tmp[256];
 	int i;
 
-	snprintf(tmp, 256, "pgraph-iter%i.dat", iter);
+	snprintf(tmp, 256, "pr-logs/pgraph%s.dat", suff);
+
+	if ( iter == 0 ) {
+		fh = fopen(tmp, "w");
+	} else {
+		fh = fopen(tmp, "a");
+	}
 
-	fh = fopen(tmp, "w");
 	if ( fh == NULL ) {
 		ERROR("Failed to open '%s'\n", tmp);
 		return;
 	}
 
-	fprintf(fh, "   fr    h    k    l  1/d(nm)     pcalc   pobs\n");
+	if ( iter == 0 ) {
+		fprintf(fh, "  Crystal    h    k    l  1/d(m)   Ipart    pcalc   pobs   iteration\n");
+	}
 
 	for ( i=0; i<n_crystals; i++ ) {
 		if ( crystal_get_user_flag(crystals[i]) != 0 ) continue;
 		write_to_pgraph(fh, crystal_get_reflections(crystals[i]), full,
-		                crystals[i], i);
+		                crystals[i], i, iter);
 	}
 
 	fclose(fh);
-
 }
 
 
@@ -690,8 +717,8 @@ static void all_residuals(Crystal **crystals, int n_crystals, RefList *full,
 
 		if ( crystal_get_user_flag(crystals[i]) ) continue;
 
-		r = residual(crystals[i], full, 0, NULL, NULL);
-		free_r = residual(crystals[i], full, 1, NULL, NULL);
+		r = residual(crystals[i], full, 0, NULL, NULL, 1);
+		free_r = residual(crystals[i], full, 1, NULL, NULL, 1);
 		log_r = log_residual(crystals[i], full, 0, NULL, NULL);
 		free_log_r = log_residual(crystals[i], full, 1, NULL, NULL);
 
@@ -720,32 +747,77 @@ static void show_all_residuals(Crystal **crystals, int n_crystals,
 }
 
 
-static void dump_parameters(const char *filename, Crystal **crystals,
-                            int n_crystals)
+struct log_qargs
 {
-	FILE *fh;
-	fh = fopen(filename, "w");
-	if ( fh == NULL ) {
-		ERROR("Couldn't open partialator.params!\n");
-	} else {
-		fprintf(fh, "  cr        OSF       relB         div"
-		            " flag                      filename  event\n");
-		int i;
-		for ( i=0; i<n_crystals; i++ ) {
-			struct image *img;
-                        char *evt_str;
-			img = crystal_get_image(crystals[i]);
-                        evt_str = get_event_string(img->event);
-			fprintf(fh, "%4i %10.5f %10.2f %8.5e %-25s %s %s\n",
-			        i, crystal_get_osf(crystals[i]),
-			        crystal_get_Bfac(crystals[i])*1e20,
-			        crystal_get_image(crystals[i])->div,
-			        str_prflag(crystal_get_user_flag(crystals[i])),
-			        img->filename, evt_str);
-                        free(evt_str);
-		}
-		fclose(fh);
-	}
+	int iter;
+	int next;
+	Crystal **crystals;
+	int n_crystals;
+	RefList *full;
+};
+
+
+struct log_args
+{
+	Crystal *cr;
+	RefList *full;
+	int iter;
+	int cnum;
+};
+
+
+static void *get_log_task(void *vp)
+{
+	struct log_qargs *qargs = vp;
+	struct log_args *task;
+
+	if ( qargs->next >= qargs->n_crystals ) return NULL;
+
+	task = malloc(sizeof(struct log_args));
+	if ( task == NULL ) return NULL;
+
+	task->cr = qargs->crystals[qargs->next];
+	task->full = qargs->full;
+	task->iter = qargs->iter;
+	task->cnum = qargs->next;
+
+	qargs->next += 20;
+	return task;
+}
+
+
+static void write_logs(void *vp, int cookie)
+{
+	struct log_args *args = vp;
+	write_specgraph(args->cr, args->full, args->iter, args->cnum);
+	write_gridscan(args->cr, args->full, args->iter, args->cnum);
+	write_test_logs(args->cr, args->full, args->iter, args->cnum);
+}
+
+
+static void done_log(void *qargs, void *vp)
+{
+	struct log_args *task = vp;
+	free(task);
+}
+
+
+static void write_logs_parallel(Crystal **crystals, int n_crystals,
+                                RefList *full, int iter, int n_threads)
+{
+	struct log_qargs qargs;
+
+
+	qargs.iter = iter;
+	qargs.next = 0;
+	qargs.full = full;
+	qargs.crystals = crystals;
+	qargs.n_crystals = n_crystals;
+
+	STATUS("Writing logs...\n");
+	run_threads(n_threads, write_logs, get_log_task, done_log, &qargs,
+	            n_crystals/20, 0, 0, 0);
+	STATUS("Done.\n");
 }
 
 
@@ -756,9 +828,10 @@ int main(int argc, char *argv[])
 	char *outfile = NULL;
 	char *sym_str = NULL;
 	SymOpList *sym;
+	SymOpList *amb;
+	SymOpList *w_sym;
 	int nthreads = 1;
 	int i;
-	struct image *images;
 	int n_iter = 10;
 	RefList *full;
 	int n_images = 0;
@@ -770,7 +843,7 @@ int main(int argc, char *argv[])
 	Stream *st;
 	Crystal **crystals;
 	char *pmodel_str = NULL;
-	PartialityModel pmodel = PMODEL_SCSPHERE;
+	PartialityModel pmodel = PMODEL_XSPHERE;
 	int min_measurements = 2;
 	char *rval;
 	int polarisation = 1;
@@ -786,6 +859,13 @@ int main(int argc, char *argv[])
 	char *csplit_fn = NULL;
 	struct custom_split *csplit = NULL;
 	double max_B = 1e-18;
+	char *rfile = NULL;
+	RefList *reference = NULL;
+	int no_logs = 0;
+	char *w_sym_str = NULL;
+	char *operator = NULL;
+	double force_bandwidth = -1.0;
+	double force_radius = -1.0;
 
 	/* Long options */
 	const struct option longopts[] = {
@@ -806,8 +886,12 @@ int main(int argc, char *argv[])
 		{"push-res",           1, NULL,                5},
 		{"res-push",           1, NULL,                5}, /* compat */
 		{"custom-split",       1, NULL,                6},
-		{"max-rel-B"   ,       1, NULL,                7},
-		{"max-rel-b"   ,       1, NULL,                7}, /* compat */
+		{"max-rel-B",          1, NULL,                7},
+		{"max-rel-b",          1, NULL,                7}, /* compat */
+		{"reference",          1, NULL,                8}, /* ssshhh! */
+		{"operator",           1, NULL,                9},
+		{"force-bandwidth",    1, NULL,               10},
+		{"force-radius",       1, NULL,               11},
 
 		{"no-scale",           0, &no_scale,           1},
 		{"no-pr",              0, &no_pr,              1},
@@ -817,6 +901,7 @@ int main(int argc, char *argv[])
 		{"polarization",       0, &polarisation,       1}, /* compat */
 		{"no-free",            0, &no_free,            1},
 		{"output-every-cycle", 0, &output_everycycle,  1},
+		{"no-logs",            0, &no_logs,            1},
 
 		{0, 0, NULL, 0}
 	};
@@ -828,7 +913,7 @@ int main(int argc, char *argv[])
 	}
 
 	/* Short options */
-	while ((c = getopt_long(argc, argv, "hi:o:g:b:y:n:j:m:",
+	while ((c = getopt_long(argc, argv, "hi:o:g:b:y:n:j:m:w:",
 	                        longopts, NULL)) != -1)
 	{
 
@@ -887,6 +972,10 @@ int main(int argc, char *argv[])
 			pmodel_str = strdup(optarg);
 			break;
 
+			case 'w' :
+			w_sym_str = strdup(optarg);
+			break;
+
 			case 2 :
 			errno = 0;
 			min_measurements = strtod(optarg, &rval);
@@ -933,6 +1022,33 @@ int main(int argc, char *argv[])
 			max_B = max_B * 1e-20;
 			break;
 
+			case 8 :
+			rfile = strdup(optarg);
+			break;
+
+			case 9 :
+			operator = strdup(optarg);
+			break;
+
+			case 10 :
+			errno = 0;
+			force_bandwidth = strtod(optarg, &rval);
+			if ( *rval != '\0' ) {
+				ERROR("Invalid value for --force-bandwidth.\n");
+				return 1;
+			}
+			break;
+
+			case 11 :
+			errno = 0;
+			force_radius = strtod(optarg, &rval);
+			if ( *rval != '\0' ) {
+				ERROR("Invalid value for --force-radius.\n");
+				return 1;
+			}
+			force_radius *= 1e9;
+			break;
+
 			case 0 :
 			break;
 
@@ -972,13 +1088,49 @@ int main(int argc, char *argv[])
 	free(sym_str);
 	if ( sym == NULL ) return 1;
 
+	if ( (w_sym_str != NULL) && (operator != NULL) ) {
+		ERROR("Specify the apparent symmetry (-w) or the operator, "
+		      "not both.\n");
+		return 1;
+	}
+
+	if ( w_sym_str == NULL ) {
+		w_sym = NULL;
+		amb = NULL;
+	} else {
+		pointgroup_warning(w_sym_str);
+		w_sym = get_pointgroup(w_sym_str);
+		free(w_sym_str);
+		if ( w_sym == NULL ) return 1;
+		amb = get_ambiguities(w_sym, sym);
+		if ( amb == NULL ) {
+			ERROR("Couldn't find ambiguity operator.\n");
+			ERROR("Check that your values for -y and -w are "
+			      "correct.\n");
+			return 1;
+		}
+
+	}
+
+	if ( operator ) {
+		amb = parse_symmetry_operations(operator);
+		if ( amb == NULL ) return 1;
+		set_symmetry_name(amb, "Ambiguity");
+	}
+
+	if ( amb != NULL ) {
+		STATUS("Indexing ambiguity resolution enabled.  "
+		       "The ambiguity operation(s) are:\n");
+		describe_symmetry(amb);
+		/* In contrast to ambigator, partialator can deal with multiple
+		 * ambiguities at once */
+	}
+
 	if ( pmodel_str != NULL ) {
 		if ( strcmp(pmodel_str, "unity") == 0 ) {
 			pmodel = PMODEL_UNITY;
-		} else if ( strcmp(pmodel_str, "scgaussian") == 0 ) {
-			pmodel = PMODEL_SCGAUSSIAN;
-		} else if ( strcmp(pmodel_str, "scsphere") == 0 ) {
-			pmodel = PMODEL_SCSPHERE;
+		} else if ( strcmp(pmodel_str, "xsphere") == 0 ) {
+			pmodel = PMODEL_XSPHERE;
 		} else if ( strcmp(pmodel_str, "random") == 0 ) {
 			pmodel = PMODEL_RANDOM;
 		} else {
@@ -989,7 +1141,22 @@ int main(int argc, char *argv[])
 
 	if ( (pmodel == PMODEL_UNITY) && !no_pr ) {
 		no_pr = 1;
-		STATUS("Setting --no-pr because we are not modelling partialities (--model=unity).\n");
+		STATUS("Setting --no-pr because we are not modelling "
+		       "partialities (--model=unity).\n");
+	}
+
+	if ( !no_logs ) {
+		int r = mkdir("pr-logs", S_IRWXU | S_IRWXG | S_IROTH | S_IXOTH);
+		if ( r ) {
+			if ( errno == EEXIST ) {
+				ERROR("A folder called 'pr-logs' exists in the "
+				      "working directory.\n");
+				ERROR("Please delete or move it first.\n");
+			} else {
+				ERROR("Failed to create pr-logs folder.\n");
+			}
+			return 1;
+		}
 	}
 
 	/* Read the custom split list (if applicable) */
@@ -1002,6 +1169,20 @@ int main(int argc, char *argv[])
 		free(csplit_fn);
 	}
 
+	if ( rfile != NULL ) {
+		RefList *rread;
+		rread = read_reflections(rfile);
+		if ( rread == NULL ) {
+			ERROR("Failed to read reference reflections\n");
+			return 1;
+		}
+		reference = asymmetric_indices(rread, sym);
+		reflist_free(rread);
+		ERROR("WARNING: Using an external reference.\n");
+		ERROR("WARNING: If you publish a structure based on the result,"
+		      " expect to have to retract your paper!\n");
+	}
+
 	gsl_set_error_handler_off();
 	rng = gsl_rng_alloc(gsl_rng_mt19937);
 
@@ -1009,7 +1190,6 @@ int main(int argc, char *argv[])
 	n_images = 0;
 	n_crystals = 0;
 	n_crystals_seen = 0;
-	images = NULL;
 	crystals = NULL;
 	if ( sparams_fn != NULL ) {
 		char line[1024];
@@ -1033,35 +1213,27 @@ int main(int argc, char *argv[])
 
 		RefList *as;
 		int i;
-		struct image *images_new;
-		struct image *cur;
+		struct image cur;
 
-		images_new = realloc(images, (n_images+1)*sizeof(struct image));
-		if ( images_new == NULL ) {
-			ERROR("Failed to allocate memory for image list.\n");
-			return 1;
-		}
-		images = images_new;
-		cur = &images[n_images];
-
-		cur->div = NAN;
-		cur->bw = NAN;
-		cur->det = NULL;
-		if ( read_chunk_2(st, cur, STREAM_READ_REFLECTIONS
-		                           | STREAM_READ_UNITCELL) != 0 ) {
+		cur.div = NAN;
+		cur.bw = NAN;
+		cur.det = NULL;
+		if ( read_chunk_2(st, &cur, STREAM_READ_REFLECTIONS
+		                            | STREAM_READ_UNITCELL) != 0 ) {
 			break;
 		}
 
-		if ( isnan(cur->div) || isnan(cur->bw) ) {
+		if ( isnan(cur.div) || isnan(cur.bw) ) {
 			ERROR("Chunk doesn't contain beam parameters.\n");
 			return 1;
 		}
 
-		for ( i=0; i<cur->n_crystals; i++ ) {
+		for ( i=0; i<cur.n_crystals; i++ ) {
 
 			Crystal *cr;
 			Crystal **crystals_new;
 			RefList *cr_refl;
+			struct image *image;
 
 			n_crystals_seen++;
 			if ( n_crystals_seen <= start_after ) continue;
@@ -1074,11 +1246,19 @@ int main(int argc, char *argv[])
 				return 1;
 			}
 			crystals = crystals_new;
-			crystals[n_crystals] = cur->crystals[i];
+			crystals[n_crystals] = cur.crystals[i];
 			cr = crystals[n_crystals];
 
-			/* Image pointer will change due to later reallocs */
-			crystal_set_image(cr, NULL);
+			image = malloc(sizeof(struct image));
+			if ( image == NULL ) {
+				ERROR("Failed to allocatea memory for image.\n");
+				return 1;
+			}
+
+			crystal_set_image(cr, image);
+			*image = cur;
+			image->n_crystals = 1;
+			image->crystals = &crystals[n_crystals];
 
 			/* This is the raw list of reflections */
 			cr_refl = crystal_get_reflections(cr);
@@ -1088,7 +1268,7 @@ int main(int argc, char *argv[])
 			if ( polarisation ) {
 				polarisation_correction(cr_refl,
 						        crystal_get_cell(cr),
-						        cur);
+						        image);
 			}
 
 			if ( !no_free ) select_free_reflections(cr_refl, rng);
@@ -1124,65 +1304,76 @@ int main(int argc, char *argv[])
 
 	close_stream(st);
 
-	/* Fill in image pointers */
-	for ( i=0; i<n_images; i++ ) {
-		int j;
-		for ( j=0; j<images[i].n_crystals; j++ ) {
-
-			Crystal *cryst;
-
-			cryst = images[i].crystals[j];
-			crystal_set_image(cryst, &images[i]);
-
-			/* Now it's safe to do the following */
-			update_predictions(cryst);
-			calculate_partialities(cryst, pmodel);
-
+	STATUS("Initial partiality calculation...\n");
+	for ( i=0; i<n_crystals; i++ ) {
+		Crystal *cr = crystals[i];
+		if ( force_bandwidth > 0.0 ) {
+			crystal_get_image(cr)->bw = force_bandwidth;
+		}
+		if ( force_radius > 0.0 ) {
+			crystal_set_profile_radius(cr, force_radius);
 		}
+		update_predictions(cr);
+		calculate_partialities(cr, pmodel);
 	}
 
 	if (csplit != NULL) check_csplit(crystals, n_crystals, csplit);
 
 	/* Make a first pass at cutting out crap */
-	STATUS("Checking patterns.\n");
+	STATUS("Early rejection...\n");
 	//early_rejection(crystals, n_crystals);
 
-	/* Initial rejection, figures of merit etc */
-	full = merge_intensities(crystals, n_crystals, nthreads, pmodel,
-	                         min_measurements, push_res, 1);
+	/* Create reference data set if we don't already have one */
+	if ( reference == NULL ) {
+		STATUS("Initial scaling...\n");
+		scale_all(crystals, n_crystals, nthreads);
+		full = merge_intensities(crystals, n_crystals, nthreads,
+		                         min_measurements, push_res, 1);
+	} else {
+		full = reference;
+	}
+
+	/* Scale everything to the reference */
+	STATUS("Initial scaling to reference...\n");
+	if ( !no_scale ) {
+		scale_all_to_reference(crystals, n_crystals, full, nthreads);
+	}
+
+	/* Check rejection and write figures of merit */
 	check_rejection(crystals, n_crystals, full, max_B);
 	show_all_residuals(crystals, n_crystals, full);
-	write_pgraph(full, crystals, n_crystals, 0);
+	if ( !no_pr && !no_logs ) {
+		write_pgraph(full, crystals, n_crystals, 0, "");
+		write_logs_parallel(crystals, n_crystals, full, 0, nthreads);
+	}
 
 	/* Iterate */
 	for ( i=0; i<n_iter; i++ ) {
 
 		STATUS("Scaling and refinement cycle %i of %i\n", i+1, n_iter);
 
-		if ( !no_scale ) {
-			scale_all(crystals, n_crystals, nthreads, pmodel);
-			reflist_free(full);
-			full = merge_intensities(crystals, n_crystals, nthreads,
-			                         pmodel, min_measurements,
-			                         push_res, 1);
+		if ( !no_pr ) {
+			refine_all(crystals, n_crystals, full, nthreads, pmodel,
+			           0, i+1, no_logs, sym, amb);
+		} else if ( !no_scale ) {
+			scale_all_to_reference(crystals, n_crystals, full, nthreads);
 		}
 
-		if ( !no_pr ) {
-			refine_all(crystals, n_crystals, full, nthreads,
-			           pmodel);
+		/* Create new reference if needed */
+		if ( reference == NULL ) {
 			reflist_free(full);
+			scale_all(crystals, n_crystals, nthreads);
 			full = merge_intensities(crystals, n_crystals, nthreads,
-			                         pmodel, min_measurements,
+			                         min_measurements,
 			                         push_res, 1);
-		}
+		} /* else full still equals reference */
 
 		check_rejection(crystals, n_crystals, full, max_B);
-		reflist_free(full);
-		full = merge_intensities(crystals, n_crystals, nthreads,
-		                         pmodel, min_measurements,
-		                         push_res, 1);
 		show_all_residuals(crystals, n_crystals, full);
-		write_pgraph(full, crystals, n_crystals, i+1);
+
+		if ( !no_logs ) {
+			write_pgraph(full, crystals, n_crystals, i+1, "");
+		}
 
 		if ( output_everycycle ) {
 
@@ -1209,17 +1400,29 @@ int main(int argc, char *argv[])
 				}
 			}
 
-			/* Dump parameters */
-			snprintf(tmp, 1024, "iter%.2d_partialator.params", i+1);
-			dump_parameters(tmp, crystals, n_crystals);
-
 		}
 	}
 
-	full = merge_intensities(crystals, n_crystals, nthreads, pmodel,
-	                         min_measurements, push_res, 1);
+	/* Final merge */
+	STATUS("Final merge...\n");
+	if ( reference == NULL ) {
+		reflist_free(full);
+		scale_all(crystals, n_crystals, nthreads);
+		full = merge_intensities(crystals, n_crystals, nthreads,
+		                         min_measurements,
+		                         push_res, 1);
+	} else {
+		scale_all_to_reference(crystals, n_crystals, reference, nthreads);
+		full = merge_intensities(crystals, n_crystals, nthreads,
+		                         min_measurements, push_res, 1);
+	}
+
+	/* Write final figures of merit (no rejection any more) */
 	show_all_residuals(crystals, n_crystals, full);
-	write_pgraph(full, crystals, n_crystals, n_iter+1);
+	if ( !no_pr && !no_logs ) {
+		write_pgraph(full, crystals, n_crystals, -1, "");
+		write_logs_parallel(crystals, n_crystals, full, -1, nthreads);
+	}
 
 	/* Output results */
 	STATUS("Writing overall results to %s\n", outfile);
@@ -1239,9 +1442,6 @@ int main(int argc, char *argv[])
 		}
 	}
 
-	/* Dump parameters */
-	dump_parameters("partialator.params", crystals, n_crystals);
-
 	/* Clean up */
 	gsl_rng_free(rng);
 	for ( i=0; i<n_crystals; i++ ) {
@@ -1249,13 +1449,9 @@ int main(int argc, char *argv[])
 		crystal_free(crystals[i]);
 	}
 	reflist_free(full);
-	free(sym);
+	free_symoplist(sym);
 	free(outfile);
 	free(crystals);
-	for ( i=0; i<n_images; i++ ) {
-		free(images[i].filename);
-	}
-	free(images);
 	free(infile);
 
 	return 0;
diff --git a/src/post-refinement.c b/src/post-refinement.c
index bf30d299..79ea187a 100644
--- a/src/post-refinement.c
+++ b/src/post-refinement.c
@@ -3,11 +3,11 @@
  *
  * Post refinement
  *
- * Copyright © 2012-2015 Deutsches Elektronen-Synchrotron DESY,
+ * Copyright © 2012-2018 Deutsches Elektronen-Synchrotron DESY,
  *                       a research centre of the Helmholtz Association.
  *
  * Authors:
- *   2010-2015 Thomas White <taw@physics.org>
+ *   2010-2018 Thomas White <taw@physics.org>
  *
  * This file is part of CrystFEL.
  *
@@ -33,11 +33,8 @@
 
 #include <stdlib.h>
 #include <assert.h>
-#include <gsl/gsl_matrix.h>
-#include <gsl/gsl_vector.h>
-#include <gsl/gsl_linalg.h>
-#include <gsl/gsl_eigen.h>
-#include <gsl/gsl_blas.h>
+#include <gsl/gsl_multimin.h>
+#include <gsl/gsl_fit.h>
 
 #include "image.h"
 #include "post-refinement.h"
@@ -46,15 +43,13 @@
 #include "geometry.h"
 #include "cell.h"
 #include "cell-utils.h"
+#include "reflist-utils.h"
+#include "scaling.h"
 
 
-/* Maximum number of iterations of NLSq to do for each image per macrocycle. */
-#define MAX_CYCLES (10)
-
 struct prdata
 {
 	int refined;
-	int n_filtered;
 };
 
 const char *str_prflag(enum prflag flag)
@@ -79,625 +74,867 @@ const char *str_prflag(enum prflag flag)
 		case PRFLAG_BIGB :
 		return "B too big";
 
+		case PRFLAG_SCALEBAD :
+		return "bad scaling";
+
 		default :
 		return "Unknown flag";
 	}
 }
 
 
-/* Returns dp(gauss)/dr at "r" */
-static double gaussian_fraction_gradient(double r, double R)
+double residual(Crystal *cr, const RefList *full, int free,
+                int *pn_used, const char *filename, int complain)
 {
-	const double ng = 2.6;
-	const double sig = R/ng;
+	Reflection *refl;
+	RefListIterator *iter;
+	int n_used = 0;
+	double num = 0.0;
+	double den = 0.0;
+	double G = crystal_get_osf(cr);
+	double B = crystal_get_Bfac(cr);
+	UnitCell *cell = crystal_get_cell(cr);
 
-	/* If the Ewald sphere isn't within the profile, the gradient is zero */
-	if ( r < -R ) return 0.0;
-	if ( r > +R ) return 0.0;
+	if ( linear_scale(full, crystal_get_reflections(cr), &G, complain) ) {
+		return GSL_NAN;
+	}
 
-	return exp(-pow(r/sig, 2.0)/2.0) / (sig*sqrt(2.0*M_PI));
-}
+	for ( refl = first_refl(crystal_get_reflections(cr), &iter);
+	      refl != NULL;
+	      refl = next_refl(refl, iter) )
+	{
+		double p, w, corr, res;
+		signed int h, k, l;
+		Reflection *match;
+		double I_full;
+		double int1, int2;
 
+		if ( free && !get_flag(refl) ) continue;
 
-/* Returns dp(sph)/dr at "r" */
-static double sphere_fraction_gradient(double r, double pr)
-{
-	double q, dpdq, dqdr;
+		get_indices(refl, &h, &k, &l);
+		res = resolution(cell, h, k, l);
+		match = find_refl(full, h, k, l);
+		if ( match == NULL ) continue;
+		I_full = get_intensity(match);
 
-	/* If the Ewald sphere isn't within the profile, the gradient is zero */
-	if ( r < -pr ) return 0.0;
-	if ( r > +pr ) return 0.0;
+		if ( get_redundancy(match) < 2 ) continue;
 
-	q = (r + pr)/(2.0*pr);
-	dpdq = 6.0*(q - q*q);
-	dqdr = 1.0 / (2.0*pr);
-	return dpdq * dqdr;
-}
+		p = get_partiality(refl);
+		//if ( p < 0.2 ) continue;
 
+		corr = G * exp(B*res*res) * get_lorentz(refl);
+		int1 = get_intensity(refl) * corr;
+		int2 = p*I_full;
+		w = p;
 
-/* Returns dp/dr at "r" */
-static double partiality_gradient(double r, double pr,
-                                  PartialityModel pmodel,
-                                  double rlow, double rhigh)
-{
-	double A, D;
+		num += fabs(int1 - int2) * w;
+		den += fabs(int1 + int2) * w/2.0;
 
-	D = rlow - rhigh;
+		n_used++;
 
-	switch ( pmodel ) {
+	}
+
+	if ( pn_used != NULL ) *pn_used = n_used;
+	return num/(den*sqrt(2));
+}
 
-		default:
-		case PMODEL_UNITY:
-		return 0.0;
 
-		case PMODEL_SCSPHERE:
-		A = sphere_fraction_gradient(r, pr)/D;
-		return 4.0*pr*A/3.0;
+static UnitCell *rotate_cell_xy(const UnitCell *cell, double ang1, double ang2)
+{
+	UnitCell *o;
+	double asx, asy, asz;
+	double bsx, bsy, bsz;
+	double csx, csy, csz;
+	double xnew, ynew, znew;
+
+	o = cell_new_from_cell(cell);
+
+	cell_get_reciprocal(o, &asx, &asy, &asz,
+	                       &bsx, &bsy, &bsz,
+	                       &csx, &csy, &csz);
+
+	/* "a" around x */
+	xnew = asx;
+	ynew = asy*cos(ang1) + asz*sin(ang1);
+	znew = -asy*sin(ang1) + asz*cos(ang1);
+	asx = xnew;  asy = ynew;  asz = znew;
+
+	/* "b" around x */
+	xnew = bsx;
+	ynew = bsy*cos(ang1) + bsz*sin(ang1);
+	znew = -bsy*sin(ang1) + bsz*cos(ang1);
+	bsx = xnew;  bsy = ynew;  bsz = znew;
+
+	/* "c" around x */
+	xnew = csx;
+	ynew = csy*cos(ang1) + csz*sin(ang1);
+	znew = -csy*sin(ang1) + csz*cos(ang1);
+	csx = xnew;  csy = ynew;  csz = znew;
+
+	/* "a" around y */
+	xnew = asx*cos(ang2) + asz*sin(ang2);
+	ynew = asy;
+	znew = -asx*sin(ang2) + asz*cos(ang2);
+	asx = xnew;  asy = ynew;  asz = znew;
+
+	/* "b" around y */
+	xnew = bsx*cos(ang2) + bsz*sin(ang2);
+	ynew = bsy;
+	znew = -bsx*sin(ang2) + bsz*cos(ang2);
+	bsx = xnew;  bsy = ynew;  bsz = znew;
+
+	/* "c" around y */
+	xnew = csx*cos(ang2) + csz*sin(ang2);
+	ynew = csy;
+	znew = -csx*sin(ang2) + csz*cos(ang2);
+	csx = xnew;  csy = ynew;  csz = znew;
+
+	cell_set_reciprocal(o, asx, asy, asz, bsx, bsy, bsz, csx, csy, csz);
+
+	return o;
+}
 
-		case PMODEL_SCGAUSSIAN:
-		A = gaussian_fraction_gradient(r, pr)/D;
-		return 4.0*pr*A/3.0;
 
+static const char *get_label(enum gparam p)
+{
+	switch ( p ) {
+		case GPARAM_ANG1 : return "First angle (radians)";
+		case GPARAM_ANG2 : return "Second angle (radians)";
+		case GPARAM_R : return "Profile radius (m^-1)";
+		case GPARAM_WAVELENGTH : return "Wavelength (Angstroms)";
+		default : return "unknown";
 	}
 }
 
 
-static double sphere_fraction_rgradient(double r, double R)
+/* We set all the step sizes to 1, then scale them.
+ * This way, the size of the simplex stays meaningful and we possibly also
+ *  avoid some roundoff errors */
+static double get_scale(enum gparam p)
 {
-	/* If the Ewald sphere isn't within the profile, the gradient is zero */
-	if ( r < -R ) return 0.0;
-	if ( r > +R ) return 0.0;
+	switch ( p ) {
+
+		case GPARAM_ANG1 : return deg2rad(0.05);
+		case GPARAM_ANG2 : return deg2rad(0.05);
+		case GPARAM_R : return 0.0005e9;
+		case GPARAM_WAVELENGTH : return 0.001e-10;
 
-	return -(3.0*r/(4.0*R*R)) * (1.0 - r*r/(R*R));
+		default : return 0.0;
+
+	}
 }
 
 
-static double gaussian_fraction_rgradient(double r, double R)
-{
-	const double ng = 2.6;
-	const double sig = R/ng;
+struct rf_priv {
+	const Crystal *cr;
+	const RefList *full;
+	enum gparam rv[32];
+	int verbose;
+	int serial;
+	gsl_vector *initial;
 
-	/* If the Ewald sphere isn't within the profile, the gradient is zero */
-	if ( r < -R ) return 0.0;
-	if ( r > +R ) return 0.0;
+	/* For freeing later */
+	gsl_vector *vals;
+	gsl_vector *step;
 
-	return -(ng*r/(sqrt(2.0*M_PI)*R*R))*exp(-r*r/(2.0*sig*sig));
+	/* So that it stays around until the end of minimisation */
+	gsl_multimin_function f;
+};
+
+
+static double get_actual_val(const gsl_vector *v, const gsl_vector *initial,
+                             enum gparam *rv, int i)
+{
+	return gsl_vector_get(v, i) * get_scale(rv[i])
+	             + gsl_vector_get(initial, i);
 }
 
 
-static double volume_fraction_rgradient(double r, double pr,
-                                       PartialityModel pmodel)
+static void apply_parameters(const gsl_vector *v, const gsl_vector *initial,
+                             enum gparam *rv, Crystal *cr)
 {
-	switch ( pmodel )
-	{
-		case PMODEL_UNITY :
-		return 1.0;
+	int i;
+	double ang1, ang2, R, lambda;
+	UnitCell *cell;
 
-		case PMODEL_SCSPHERE :
-		return sphere_fraction_rgradient(r, pr);
+	/* Default parameters if not used in refinement */
+	ang1 = 0.0;
+	ang2 = 0.0;
+	R = crystal_get_profile_radius(cr);
+	lambda = crystal_get_image(cr)->lambda;
 
-		case PMODEL_SCGAUSSIAN :
-		return gaussian_fraction_rgradient(r, pr);
+	for ( i=0; i<v->size; i++ ) {
 
-		default :
-		ERROR("No pmodel in volume_fraction_rgradient!\n");
-		return 1.0;
-	}
-}
+		double val;
 
+		val = get_actual_val(v, initial, rv, i);
 
-static double volume_fraction(double rlow, double rhigh, double pr,
-                              PartialityModel pmodel)
-{
-	switch ( pmodel )
-	{
-		case PMODEL_UNITY :
-		return 1.0;
+		switch ( rv[i] ) {
 
-		case PMODEL_SCSPHERE :
-		return sphere_fraction(rlow, rhigh, pr);
+			case GPARAM_ANG1 :
+			ang1 = val;
+			break;
 
-		case PMODEL_SCGAUSSIAN :
-		return gaussian_fraction(rlow, rhigh, pr);
+			case GPARAM_ANG2 :
+			ang2 = val;
+			break;
 
-		default :
-		ERROR("No pmodel in volume_fraction!\n");
-		return 1.0;
-	}
-}
+			case GPARAM_R :
+			R = val;
+			break;
 
+			case GPARAM_WAVELENGTH :
+			lambda = val;
+			break;
 
-/* Return the gradient of "fx" wrt parameter 'k' given the current
- * status of the crystal. */
-double gradient(Crystal *cr, int k, Reflection *refl, PartialityModel pmodel)
-{
-	double glow, ghigh;
-	double rlow, rhigh, p;
-	struct image *image = crystal_get_image(cr);
-	double R = crystal_get_profile_radius(cr);
-	double gr;
+			default :
+			ERROR("Don't understand parameter %i\n", rv[i]);
+			break;
 
-	get_partial(refl, &rlow, &rhigh, &p);
+		}
+	}
 
-	if ( k == GPARAM_R ) {
+	cell = rotate_cell_xy(crystal_get_cell_const(cr), ang1, ang2);
+	crystal_set_cell(cr, cell);
 
-		double Rglow, Rghigh;
-		double D, psph;
+	crystal_set_profile_radius(cr, R);
+	crystal_get_image(cr)->lambda = lambda;
+}
 
-		D = rlow - rhigh;
-		psph = volume_fraction(rlow, rhigh, R, pmodel);
 
-		Rglow = volume_fraction_rgradient(rlow, R, pmodel);
-		Rghigh = volume_fraction_rgradient(rhigh, R, pmodel);
+static double residual_f(const gsl_vector *v, void *pp)
+{
+	struct rf_priv *pv = pp;
+	RefList *list;
+	struct image im;
+	Crystal *cr;
+	UnitCell *cell;
+	double res;
+	int i;
 
-		gr = 4.0*psph/(3.0*D) + (4.0*R/(3.0*D))*(Rglow - Rghigh);
-		return gr;
+	for ( i=0; i<v->size; i++ ) {
+		if ( gsl_vector_get(v, i) > 100.0 ) return GSL_NAN;
+	}
 
+	cr = crystal_copy(pv->cr);
+	cell = cell_new_from_cell(crystal_get_cell(cr));
+	if ( cell == NULL ) return GSL_NAN;
+	crystal_set_cell(cr, cell);
+	im = *crystal_get_image(cr);
+	crystal_set_image(cr, &im);
+	apply_parameters(v, pv->initial, pv->rv, cr);
+
+	if ( fabs(crystal_get_profile_radius(cr)) > 5e9 ) {
+		crystal_free(cr);
+		if ( pv->verbose ) STATUS("radius > 5e9\n");
+		return GSL_NAN;
 	}
 
-	/* Calculate the gradient of partiality wrt excitation error. */
-	glow = partiality_gradient(rlow, R, pmodel, rlow, rhigh);
-	ghigh = partiality_gradient(rhigh, R, pmodel, rlow, rhigh);
+	/* Can happen with grid scans and certain --force-radius values */
+	if ( fabs(crystal_get_profile_radius(cr)) < 0.0000001e9 ) {
+		crystal_free(cr);
+		if ( pv->verbose ) STATUS("radius very small\n");
+		return GSL_NAN;
+	}
 
-	if ( k == GPARAM_DIV ) {
+	if ( im.lambda <= 0.0 ) {
+		crystal_free(cr);
+		if ( pv->verbose ) STATUS("lambda < 0\n");
+		return GSL_NAN;
+	}
 
-		double D, psph, ds;
-		signed int hs, ks, ls;
+	list = copy_reflist(crystal_get_reflections(cr));
+	crystal_set_reflections(cr, list);
 
-		D = rlow - rhigh;
-		psph = volume_fraction(rlow, rhigh, R, pmodel);
-		get_symmetric_indices(refl, &hs, &ks, &ls);
-		ds = 2.0 * resolution(crystal_get_cell(cr), hs, ks, ls);
+	update_predictions(cr);
+	calculate_partialities(cr, PMODEL_XSPHERE);
 
-		gr = (ds/2.0)*(glow+ghigh) - 4.0*R*psph*ds/(3.0*D*D);
-		return gr;
+	res = residual(cr, pv->full, 0, NULL, NULL, 0);
 
-	}
+	cell_free(crystal_get_cell(cr));
+	reflist_free(crystal_get_reflections(cr));
+	crystal_free(cr);
 
-	gr = r_gradient(crystal_get_cell(cr), k, refl, image) * (glow-ghigh);
-	return gr;
+	return res;
 }
 
 
-static void apply_cell_shift(UnitCell *cell, int k, double shift)
+static double get_initial_param(Crystal *cr, enum gparam p)
 {
-	double asx, asy, asz;
-	double bsx, bsy, bsz;
-	double csx, csy, csz;
-	double as, bs, cs;
+	switch ( p ) {
 
-	cell_get_reciprocal(cell, &asx, &asy, &asz,
-	                          &bsx, &bsy, &bsz,
-	                          &csx, &csy, &csz);
-	as = modulus(asx, asy, asz);
-	bs = modulus(bsx, bsy, bsz);
-	cs = modulus(csx, csy, csz);
+		case GPARAM_ANG1 : return 0.0;
+		case GPARAM_ANG2 : return 0.0;
+		case GPARAM_R : return crystal_get_profile_radius(cr);
+		case GPARAM_WAVELENGTH : return crystal_get_image(cr)->lambda;
 
-	/* Forbid any step which looks too large */
-	switch ( k )
-	{
-		case GPARAM_ASX :
-		case GPARAM_ASY :
-		case GPARAM_ASZ :
-		if ( fabs(shift) > 0.1*as ) return;
-		break;
-
-		case GPARAM_BSX :
-		case GPARAM_BSY :
-		case GPARAM_BSZ :
-		if ( fabs(shift) > 0.1*bs ) return;
-		break;
-
-		case GPARAM_CSX :
-		case GPARAM_CSY :
-		case GPARAM_CSZ :
-		if ( fabs(shift) > 0.1*cs ) return;
-		break;
-	}
-
-	switch ( k )
-	{
-		case GPARAM_ASX :  asx += shift;  break;
-		case GPARAM_ASY :  asy += shift;  break;
-		case GPARAM_ASZ :  asz += shift;  break;
-		case GPARAM_BSX :  bsx += shift;  break;
-		case GPARAM_BSY :  bsy += shift;  break;
-		case GPARAM_BSZ :  bsz += shift;  break;
-		case GPARAM_CSX :  csx += shift;  break;
-		case GPARAM_CSY :  csy += shift;  break;
-		case GPARAM_CSZ :  csz += shift;  break;
-	}
-
-	cell_set_reciprocal(cell, asx, asy, asz,
-	                          bsx, bsy, bsz,
-	                          csx, csy, csz);
+		default: return 0.0;
+
+	}
 }
 
 
-/* Apply the given shift to the 'k'th parameter of 'image'. */
-static void apply_shift(Crystal *cr, int k, double shift)
+static int check_angle_shifts(gsl_vector *cur, const gsl_vector *initial,
+                              enum gparam *rv, struct rf_priv *residual_f_priv)
 {
-	double t;
-	struct image *image = crystal_get_image(cr);
-
-	switch ( k ) {
-
-		case GPARAM_DIV :
-		if ( shift > 0.1*image->div ) return;
-		image->div += shift;
-		if ( image->div < 0.0 ) image->div = 0.0;
-		break;
-
-		case GPARAM_R :
-		t = crystal_get_profile_radius(cr);
-		if ( shift > 0.1*t ) return;
-		t += shift;
-		crystal_set_profile_radius(cr, t);
-		break;
-
-		case GPARAM_ASX :
-		case GPARAM_ASY :
-		case GPARAM_ASZ :
-		case GPARAM_BSX :
-		case GPARAM_BSY :
-		case GPARAM_BSZ :
-		case GPARAM_CSX :
-		case GPARAM_CSY :
-		case GPARAM_CSZ :
-		apply_cell_shift(crystal_get_cell(cr), k, shift);
-		break;
+	int i = 0;
+	double ang = 0.0;
 
-		default :
-		ERROR("No shift defined for parameter %i\n", k);
-		abort();
+	while ( rv[i] != GPARAM_EOL ) {
+		if ( (rv[i] == GPARAM_ANG1) || (rv[i] == GPARAM_ANG2) ) {
+			ang += fabs(get_actual_val(cur, initial, rv, i));
+		}
+		rv++;
+	}
 
+	if ( rad2deg(ang) > 5.0 ) {
+		ERROR("More than 5 degrees total rotation!\n");
+		residual_f_priv->verbose = 1;
+		double res = residual_f(cur, residual_f_priv);
+		STATUS("residual after rotation = %e\n", res);
+		residual_f_priv->verbose = 2;
+		res = residual_f(initial, residual_f_priv);
+		STATUS("residual before rotation = %e\n", res);
+		return 1;
 	}
+	return 0;
 }
 
 
-/* Perform one cycle of post refinement on 'image' against 'full' */
-static double pr_iterate(Crystal *cr, const RefList *full,
-                         PartialityModel pmodel,
-                         int *n_filtered, int verbose)
+static RefList *reindex_reflections(RefList *input, SymOpList *amb,
+                                    SymOpList *sym, int idx)
 {
-	gsl_matrix *M;
-	gsl_vector *v;
-	gsl_vector *shifts;
-	int param;
+	RefList *n;
 	Reflection *refl;
 	RefListIterator *iter;
-	RefList *reflections;
-	double max_shift;
-	int nref = 0;
-	int num_params = 0;
-	enum gparam rv[32];
-	double G, B;
-
-	if ( n_filtered != NULL ) *n_filtered = 0;
-
-	rv[num_params++] = GPARAM_ASX;
-	rv[num_params++] = GPARAM_ASY;
-	rv[num_params++] = GPARAM_ASZ;
-	rv[num_params++] = GPARAM_BSX;
-	rv[num_params++] = GPARAM_BSY;
-	rv[num_params++] = GPARAM_BSZ;
-	rv[num_params++] = GPARAM_CSX;
-	rv[num_params++] = GPARAM_CSY;
-	rv[num_params++] = GPARAM_CSZ;
-//	rv[num_params++] = GPARAM_R;
-//	rv[num_params++] = GPARAM_DIV;
-
-	M = gsl_matrix_calloc(num_params, num_params);
-	v = gsl_vector_calloc(num_params);
-
-	reflections = crystal_get_reflections(cr);
-	G = crystal_get_osf(cr);
-	B = crystal_get_Bfac(cr);
-
-	/* Post-refinement terms */
-	for ( refl = first_refl(reflections, &iter);
+
+	n = reflist_new();
+
+	for ( refl = first_refl(input, &iter);
 	      refl != NULL;
 	      refl = next_refl(refl, iter) )
 	{
-		signed int ha, ka, la;
-		double I_full, delta_I, esd;
-		double w;
-		double I_partial;
-		int k;
-		double p, L, s;
-		double fx;
-		Reflection *match;
-		double gradients[num_params];
-
-		if ( get_flag(refl) ) continue;
+		signed int h, k, l;
+		Reflection *rn;
 
-		get_indices(refl, &ha, &ka, &la);
-		match = find_refl(full, ha, ka, la);
-		if ( match == NULL ) continue;
-		I_full = get_intensity(match);
+		get_indices(refl, &h, &k, &l);
+		get_equiv(amb, NULL, idx, h, k, l, &h, &k, &l);
+		get_asymm(sym, h, k, l, &h, &k, &l);
+		rn = add_refl(n, h, k, l);
+		copy_data(rn, refl);
+
+		get_symmetric_indices(rn, &h, &k, &l);
+		get_equiv(amb, NULL, idx, h, k, l, &h, &k, &l);
+		set_symmetric_indices(rn, h, k, l);
+	}
 
-		if ( get_redundancy(match) < 2 ) continue;
+	return n;
+}
 
-		p = get_partiality(refl);
-		L = get_lorentz(refl);
-		I_partial = get_intensity(refl);
-		esd = get_esd_intensity(refl);
-		s = resolution(crystal_get_cell(cr), ha, ka, la);
 
-		if ( I_partial < 3.0*esd ) continue;
+static void reindex_cell(UnitCell *cell, SymOpList *amb, int idx)
+{
+	signed int h, k, l;
+	struct rvec na, nb, nc;
+	struct rvec as, bs, cs;
+
+	cell_get_reciprocal(cell, &as.u, &as.v, &as.w,
+	                          &bs.u, &bs.v, &bs.w,
+	                          &cs.u, &cs.v, &cs.w);
+
+	get_equiv(amb, NULL, idx, 1, 0, 0, &h, &k, &l);
+	na.u = as.u*h + bs.u*k + cs.u*l;
+	na.v = as.v*h + bs.v*k + cs.v*l;
+	na.w = as.w*h + bs.w*k + cs.w*l;
+
+	get_equiv(amb, NULL, idx, 0, 1, 0, &h, &k, &l);
+	nb.u = as.u*h + bs.u*k + cs.u*l;
+	nb.v = as.v*h + bs.v*k + cs.v*l;
+	nb.w = as.w*h + bs.w*k + cs.w*l;
+
+	get_equiv(amb, NULL, idx, 0, 0, 1, &h, &k, &l);
+	nc.u = as.u*h + bs.u*k + cs.u*l;
+	nc.v = as.v*h + bs.v*k + cs.v*l;
+	nc.w = as.w*h + bs.w*k + cs.w*l;
+
+	cell_set_reciprocal(cell, na.u, na.v, na.w,
+	                          nb.u, nb.v, nb.w,
+	                          nc.u, nc.v, nc.w);
+}
 
-		/* Calculate the weight for this reflection */
-		w = (s/1e9)*(s/1e9) / (esd*esd);
 
-		/* Calculate all gradients for this reflection */
-		for ( k=0; k<num_params; k++ ) {
-			gradients[k] = gradient(cr, rv[k], refl, pmodel);
-			gradients[k] *= exp(G)*exp(-B*s*s)*I_full/L;
-		}
+void try_reindex(Crystal *crin, const RefList *full,
+                 SymOpList *sym, SymOpList *amb)
+{
+	RefList *list;
+	Crystal *cr;
+	UnitCell *cell;
+	double residual_original, residual_flipped;
+	int idx, n;
 
-		for ( k=0; k<num_params; k++ ) {
+	if ( sym == NULL || amb == NULL ) return;
 
-			int g;
-			double v_c, v_curr;
+	residual_original = residual(crin, full, 0, NULL, NULL, 0);
 
-			for ( g=0; g<num_params; g++ ) {
+	cr = crystal_copy(crin);
 
-				double M_c, M_curr;
+	n = num_equivs(amb, NULL);
 
-				/* Matrix is symmetric */
-				if ( g > k ) continue;
+	for ( idx=0; idx<n; idx++ ) {
 
-				M_c = w * gradients[g] * gradients[k];
+		cell = cell_new_from_cell(crystal_get_cell(crin));
+		if ( cell == NULL ) return;
+		reindex_cell(cell, amb, idx);
+		crystal_set_cell(cr, cell);
 
-				M_curr = gsl_matrix_get(M, k, g);
-				gsl_matrix_set(M, k, g, M_curr + M_c);
-				gsl_matrix_set(M, g, k, M_curr + M_c);
+		list = reindex_reflections(crystal_get_reflections(crin),
+		                           amb, sym, idx);
+		crystal_set_reflections(cr, list);
 
-			}
+		update_predictions(cr);
+		calculate_partialities(cr, PMODEL_XSPHERE);
 
-			fx = exp(G)*p*exp(-B*s*s)*I_full/L;
-			delta_I = I_partial - fx;
-			v_c = w * delta_I * gradients[k];
-			v_curr = gsl_vector_get(v, k);
-			gsl_vector_set(v, k, v_curr + v_c);
+		residual_flipped = residual(cr, full, 0, NULL, NULL, 1);
 
+		if ( residual_flipped < residual_original ) {
+			crystal_set_cell(crin, cell);
+			crystal_set_reflections(crin, list);
+			residual_original = residual_flipped;
+		} else {
+			cell_free(crystal_get_cell(cr));
+			reflist_free(crystal_get_reflections(cr));
 		}
-
-		nref++;
 	}
-	if ( verbose ) {
-		//STATUS("The original equation:\n");
-		//show_matrix_eqn(M, v);
-		STATUS("%i reflections went into the equations.\n", nref);
-	}
-
-	if ( nref < num_params ) {
-		crystal_set_user_flag(cr, PRFLAG_FEWREFL);
-		gsl_matrix_free(M);
-		gsl_vector_free(v);
-		return 0.0;
-	}
-
-	max_shift = 0.0;
-	shifts = solve_svd(v, M, n_filtered, 0);
-	if ( shifts != NULL ) {
-
-		for ( param=0; param<num_params; param++ ) {
-			double shift = gsl_vector_get(shifts, param);
-			if ( verbose ) STATUS("Shift %i: %e\n", param, shift);
-			if ( isnan(shift) ) {
-				//ERROR("NaN shift parameter %i (image ser %i), "
-				//       "%i reflections.\n", rv[param],
-				//       crystal_get_image(cr)->serial,
-				//       nref);
-			} else {
-				apply_shift(cr, rv[param], shift);
-			}
-			if ( fabs(shift) > max_shift ) max_shift = fabs(shift);
-		}
 
+	crystal_free(cr);
+}
+
+
+void write_test_logs(Crystal *crystal, const RefList *full,
+                     signed int cycle, int serial)
+{
+	FILE *fh;
+	struct image *image = crystal_get_image(crystal);
+	char tmp[256];
+	char ins[5];
+
+	snprintf(tmp, 256, "pr-logs/parameters-crystal%i.dat", serial);
+
+	if ( cycle == 0 ) {
+		fh = fopen(tmp, "w");
 	} else {
-		crystal_set_user_flag(cr, PRFLAG_SOLVEFAIL);
+		fh = fopen(tmp, "a");
+	}
+
+	if ( fh == NULL ) {
+		ERROR("Failed to open '%s'\n", tmp);
+		return;
 	}
 
-	gsl_matrix_free(M);
-	gsl_vector_free(v);
-	gsl_vector_free(shifts);
+	if ( cycle == 0 ) {
+		fprintf(fh, "Image: %s %s\n",
+		        image->filename, get_event_string(image->event));
+	}
 
-	return max_shift;
+	if ( cycle >= 0 ) {
+		snprintf(ins, 4, "%i", cycle);
+	} else {
+		ins[0] = 'F';
+		ins[1] = '\0';
+	}
+
+	fprintf(fh, "%s rlp_size = %e\n", ins, crystal_get_profile_radius(crystal)/1e10);
+	fprintf(fh, "%s mosaicity = %e\n", ins, crystal_get_mosaicity(crystal));
+	fprintf(fh, "%s wavelength = %f A\n", ins, image->lambda*1e10);
+	fprintf(fh, "%s bandwidth = %f\n", ins, image->bw);
+	fprintf(fh, "%s my scale factor = %e\n", ins, crystal_get_osf(crystal));
+
+	double asx, asy, asz, bsx, bsy, bsz, csx, csy, csz;
+	cell_get_reciprocal(crystal_get_cell(crystal), &asx, &asy, &asz,
+	                                               &bsx, &bsy, &bsz,
+	                                               &csx, &csy, &csz);
+	fprintf(fh, "%s %f, %f, %f, %f, %f, %f, %f, %f, %f\n", ins,
+	        asx/1e10, bsx/1e10, csx/1e10,
+	        asy/1e10, bsy/1e10, csy/1e10,
+	        asz/1e10, bsz/1e10, csz/1e10);
+	fclose(fh);
 }
 
 
-double residual(Crystal *cr, const RefList *full, int free,
-                int *pn_used, const char *filename)
+void write_specgraph(Crystal *crystal, const RefList *full,
+                     signed int cycle, int serial)
 {
-	double dev = 0.0;
-	double G, B;
+	FILE *fh;
+	char tmp[256];
 	Reflection *refl;
 	RefListIterator *iter;
-	FILE *fh = NULL;
-	int n_used = 0;
+	double G = crystal_get_osf(crystal);
+	double B = crystal_get_Bfac(crystal);
+	UnitCell *cell;
+	struct image *image = crystal_get_image(crystal);
+	char ins[5];
 
-	if ( filename != NULL ) {
-		fh = fopen(filename, "a");
-		if ( fh == NULL ) {
-			ERROR("Failed to open '%s'\n", filename);
-		}
+	snprintf(tmp, 256, "pr-logs/specgraph-crystal%i.dat", serial);
+
+	if ( cycle == 0 ) {
+		fh = fopen(tmp, "w");
+	} else {
+		fh = fopen(tmp, "a");
 	}
 
-	G = crystal_get_osf(cr);
-	B = crystal_get_Bfac(cr);
+	if ( fh == NULL ) {
+		ERROR("Failed to open '%s'\n", tmp);
+		return;
+	}
 
-	for ( refl = first_refl(crystal_get_reflections(cr), &iter);
+	if ( cycle == 0 ) {
+		fprintf(fh, "Image: %s %s\n",
+		        image->filename, get_event_string(image->event));
+		fprintf(fh, "khalf/m   1/d(m)  pcalc    pobs   iteration  h  k  l\n");
+	}
+
+	cell = crystal_get_cell(crystal);
+
+	if ( cycle >= 0 ) {
+		snprintf(ins, 4, "%i", cycle);
+	} else {
+		ins[0] = 'F';
+		ins[1] = '\0';
+	}
+
+	for ( refl = first_refl(crystal_get_reflections(crystal), &iter);
 	      refl != NULL;
 	      refl = next_refl(refl, iter) )
 	{
-		double p, L, s, w;
+		double corr, Ipart, Ifull, pobs, pcalc;
+		double res;
 		signed int h, k, l;
 		Reflection *match;
-		double esd, I_full, I_partial;
-		double fx, dc;
-
-		if ( free && !get_flag(refl) ) continue;
 
 		get_indices(refl, &h, &k, &l);
+		res = resolution(cell, h, k, l);
+
 		match = find_refl(full, h, k, l);
 		if ( match == NULL ) continue;
-		I_full = get_intensity(match);
 
-		if ( get_redundancy(match) < 2 ) continue;
+		/* Don't calculate pobs if reference reflection is weak */
+		if ( fabs(get_intensity(match)) / get_esd_intensity(match) < 3.0 ) continue;
 
-		p = get_partiality(refl);
-		L = get_lorentz(refl);
-		I_partial = get_intensity(refl);
-		esd = get_esd_intensity(refl);
-		s = resolution(crystal_get_cell(cr), h, k, l);
-
-		if ( I_partial < 3.0*esd ) continue;
+		corr = G * exp(B*res*res) * get_lorentz(refl);
+		Ipart = get_intensity(refl) * corr;
+		Ifull = get_intensity(match);
+		pobs = Ipart / Ifull;
+		pcalc = get_partiality(refl);
 
-		fx = exp(G)*p*exp(-B*s*s)*I_full/L;
-		dc = I_partial - fx;
-		w = (s/1e9)*(s/1e9)/(esd*esd);
-		dev += w*dc*dc;
-		n_used++;
+		fprintf(fh, "%e   %e   %f   %f   %s  %4i %4i %4i\n",
+		        get_khalf(refl), 2.0*res, pcalc, pobs, ins, h, k, l);
 
-		if ( fh != NULL ) {
-			fprintf(fh, "%4i %4i %4i %e %e\n",
-			        h, k, l, s, dev);
-		}
 	}
 
-	if ( fh != NULL ) fclose(fh);
-
-	if ( pn_used != NULL ) *pn_used = n_used;
-	return dev;
+	fclose(fh);
 }
 
 
-static void write_residual_graph(Crystal *cr, const RefList *full)
+static gsl_multimin_fminimizer *setup_minimiser(Crystal *cr, const RefList *full,
+                                                int verbose, int serial,
+                                                struct rf_priv *priv)
 {
+	gsl_multimin_fminimizer *min;
+	int n_params = 0;
 	int i;
-	double asx, asy, asz;
-	double bsx, bsy, bsz;
-	double csx, csy, csz;
-	double a;
-	double step;
-	double orig_asx;
-	FILE *fh;
-	UnitCell *cell;
 
-	cell = crystal_get_cell(cr);
+	/* The parameters to be refined */
+	priv->rv[n_params++] = GPARAM_ANG1;
+	priv->rv[n_params++] = GPARAM_ANG2;
+	priv->rv[n_params++] = GPARAM_R;
+	priv->rv[n_params++] = GPARAM_WAVELENGTH;
+	priv->rv[n_params] = GPARAM_EOL;  /* End of list */
+
+	priv->cr = cr;
+	priv->full = full;
+	priv->verbose = verbose;
+	priv->serial = serial;
+
+	priv->f.f = residual_f;
+	priv->f.n = n_params;
+	priv->f.params = priv;
+
+	priv->initial = gsl_vector_alloc(n_params);
+	priv->vals = gsl_vector_alloc(n_params);
+	priv->step = gsl_vector_alloc(n_params);
+
+	for ( i=0; i<n_params; i++ ) {
+		gsl_vector_set(priv->initial, i, get_initial_param(cr, priv->rv[i]));
+		gsl_vector_set(priv->vals, i, 0.0);
+		gsl_vector_set(priv->step, i, 1.0);
+	}
 
-	fh = fopen("residual-plot.dat", "w");
+	min = gsl_multimin_fminimizer_alloc(gsl_multimin_fminimizer_nmsimplex2,
+	                                    n_params);
+	gsl_multimin_fminimizer_set(min, &priv->f, priv->vals, priv->step);
 
-	cell_get_reciprocal(cell, &asx, &asy, &asz,
-	                          &bsx, &bsy, &bsz,
-	                          &csx, &csy, &csz);
-	a = modulus(asx, asy, asz);
-	orig_asx = asx;
+	return min;
+}
 
-	step = a/100.0;
 
-	for ( i=-50; i<=50; i++ ) {
+static void write_grid(Crystal *cr, const RefList *full,
+                       signed int cycle, int serial,
+                       const enum gparam par1, const enum gparam par2,
+                       const char *name)
+{
+	FILE *fh;
+	char fn[64];
+	char ins[5];
+	gsl_multimin_fminimizer *min;
+	struct rf_priv priv;
+	int idx1, idx2;
+	int i;
 
-		double res;
-		int n;
-		asx = orig_asx + (i*step);
-		cell_set_reciprocal(cell, asx, asy, asz,
-		                          bsx, bsy, bsz,
-		                          csx, csy, csz);
-		update_predictions(cr);
-		calculate_partialities(cr, PMODEL_SCSPHERE);
-		res = residual(cr, full, 0, &n, NULL);
-		fprintf(fh, "%i %e %e %i\n", i, asx, res, n);
+	min = setup_minimiser(cr, full, 0, serial, &priv);
+
+	idx1 = 99;
+	idx2 = 99;
+	for ( i=0; i<priv.f.n; i++ ) {
+		if ( priv.rv[i] == par1 ) idx1 = i;
+		if ( priv.rv[i] == par2 ) idx2 = i;
 	}
+	assert(idx1 != 99);
+	assert(idx2 != 99);
 
-	cell_set_reciprocal(cell, orig_asx, asy, asz,
-	                          bsx, bsy, bsz,
-	                          csx, csy, csz);
-	update_predictions(cr);
-	calculate_partialities(cr, PMODEL_SCSPHERE);
+	if ( cycle >= 0 ) {
+		snprintf(ins, 4, "%i", cycle);
+	} else {
+		ins[0] = 'F';
+		ins[1] = '\0';
+	}
+
+	snprintf(fn, 63, "pr-logs/grid-crystal%i-cycle%s-%s.dat",
+	         serial, ins, name);
+	fh = fopen(fn, "w");
+	if ( fh != NULL ) {
+		double v1, v2;
+		fprintf(fh, "%e %e %e %s\n",
+		        -5.0*get_scale(par1)+get_initial_param(cr, par1),
+		         5.0*get_scale(par1)+get_initial_param(cr, par1),
+		         get_initial_param(cr, par1),
+		         get_label(par1));
+		fprintf(fh, "%e %e %e %s\n",
+		        -5.0*get_scale(par2)+get_initial_param(cr, par2),
+		         5.0*get_scale(par2)+get_initial_param(cr, par2),
+		         get_initial_param(cr, par2),
+		        get_label(par2));
+		for ( v2=-5.0; v2<=5.0; v2+=0.25 ) {
+			int first=1;
+			for ( v1=-5.0; v1<=5.0; v1+=0.25 ) {
+				double res;
+				gsl_vector_set(min->x, idx1, v1);
+				gsl_vector_set(min->x, idx2, v2);
+				res = residual_f(min->x, &priv);
+				if ( !first ) fprintf(fh, " ");
+				first = 0;
+				fprintf(fh, "%e", res);
+			}
+			fprintf(fh, "\n");
+		}
+	}
 	fclose(fh);
+
+	gsl_multimin_fminimizer_free(min);
+	gsl_vector_free(priv.initial);
+	gsl_vector_free(priv.vals);
+	gsl_vector_free(priv.step);
+}
+
+
+void write_gridscan(Crystal *cr, const RefList *full,
+                    signed int cycle, int serial)
+{
+	write_grid(cr, full, cycle, serial,
+	           GPARAM_ANG1, GPARAM_ANG2, "ang1-ang2");
+	write_grid(cr, full, cycle, serial,
+	           GPARAM_ANG1, GPARAM_WAVELENGTH, "ang1-wave");
+	write_grid(cr, full, cycle, serial,
+	           GPARAM_R, GPARAM_WAVELENGTH, "R-wave");
 }
 
 
 static void do_pr_refine(Crystal *cr, const RefList *full,
-                         PartialityModel pmodel, int verbose)
+                         PartialityModel pmodel, int verbose, int serial,
+                         int cycle, int write_logs,
+                         SymOpList *sym, SymOpList *amb)
 {
-	int i, done;
-	double old_dev;
-	UnitCell *cell = crystal_get_cell(cr);
+	gsl_multimin_fminimizer *min;
+	struct rf_priv priv;
+	int n_iter = 0;
+	int status;
+	int r;
+	double G;
+	double residual_start, residual_free_start;
+	FILE *fh = NULL;
+
+	try_reindex(cr, full, sym, amb);
 
-	old_dev = residual(cr, full, 0, 0, NULL);
+	residual_start = residual(cr, full, 0, NULL, NULL, 0);
+	residual_free_start = residual(cr, full, 1, NULL, NULL, 0);
 
 	if ( verbose ) {
-		double asx, asy, asz;
-		double bsx, bsy, bsz;
-		double csx, csy, csz;
-		cell_get_reciprocal(cell, &asx, &asy, &asz,
-		                          &bsx, &bsy, &bsz,
-		                          &csx, &csy, &csz);
-		STATUS("Initial asx = %e\n", asx);
-		STATUS("PR initial  dev =  %10.5e, free dev = %10.5e\n",
-		       old_dev, residual(cr, full, 1, NULL, NULL));
-	}
-
-	i = 0;
-	done = 0;
-	do {
+		STATUS("\nPR initial: dev = %10.5e, free dev = %10.5e\n",
+		       residual_start, residual_free_start);
+	}
 
-		double dev;
+	min = setup_minimiser(cr, full, verbose, serial, &priv);
 
-		pr_iterate(cr, full, pmodel, NULL, verbose);
+	if ( verbose ) {
+		double res = residual_f(min->x, &priv);
+		double size = gsl_multimin_fminimizer_size(min);
+		STATUS("At start: %f %f %f %f ----> %f %f %e %f residual = %e size %f\n",
+		       gsl_vector_get(min->x, 0),
+		       gsl_vector_get(min->x, 1),
+		       gsl_vector_get(min->x, 2),
+		       gsl_vector_get(min->x, 3),
+		       rad2deg(get_actual_val(min->x, priv.initial, priv.rv, 0)),
+		       rad2deg(get_actual_val(min->x, priv.initial, priv.rv, 1)),
+		       get_actual_val(min->x, priv.initial, priv.rv, 2),
+		       get_actual_val(min->x, priv.initial, priv.rv, 3)*1e10,
+		       res, size);
+	}
 
-		update_predictions(cr);
-		calculate_partialities(cr, pmodel);
+	if ( write_logs ) {
 
-		dev = residual(cr, full, 0, 0, NULL);
-		if ( fabs(dev - old_dev) < dev*0.0001 ) done = 1;
+		char fn[64];
+
+		snprintf(fn, 63, "pr-logs/crystal%i-cycle%i.log", serial, cycle);
+		fh = fopen(fn, "w");
+		if ( fh != NULL ) {
+			fprintf(fh, "iteration  RtoReference  CCtoReference  nref  "
+			            "ang1     ang2    radius    wavelength");
+			double res = residual_f(min->x, &priv);
+			fprintf(fh, "%5i %10.8f  %10.8f  %5i  %10.8f %10.8f  %e  %e\n",
+			        n_iter, res, 0.0, 0,
+			        rad2deg(get_actual_val(min->x, priv.initial, priv.rv, 0)),
+			        rad2deg(get_actual_val(min->x, priv.initial, priv.rv, 1)),
+			        get_actual_val(min->x, priv.initial, priv.rv, 2),
+			        get_actual_val(min->x, priv.initial, priv.rv, 3)*1e10);
+		}
+
+	}
+
+	do {
+		double res;
+
+		n_iter++;
+
+		status = gsl_multimin_fminimizer_iterate(min);
+		if ( status ) break;
+
+		res = residual_f(min->x, &priv);
+		if ( isnan(res) ) {
+			status = GSL_ENOPROG;
+			break;
+		}
 
 		if ( verbose ) {
-			STATUS("PR iter %2i: dev = %10.5e, free dev = %10.5e\n",
-			       i+1, dev, residual(cr, full, 1, NULL, NULL));
+			double res = residual_f(min->x, &priv);
+			double size = gsl_multimin_fminimizer_size(min);
+			STATUS("%f %f %f %f ----> %f %f %e %f residual = %e size %f\n",
+			       gsl_vector_get(min->x, 0),
+			       gsl_vector_get(min->x, 1),
+			       gsl_vector_get(min->x, 2),
+			       gsl_vector_get(min->x, 3),
+			       rad2deg(get_actual_val(min->x, priv.initial, priv.rv, 0)),
+			       rad2deg(get_actual_val(min->x, priv.initial, priv.rv, 1)),
+			       get_actual_val(min->x, priv.initial, priv.rv, 2),
+			       get_actual_val(min->x, priv.initial, priv.rv, 3)*1e10,
+			       res, size);
 		}
 
-		i++;
-		old_dev = dev;
+		if ( fh != NULL ) {
+			fprintf(fh, "%5i %10.8f  %10.8f  %5i  %10.8f %10.8f  %e  %e\n",
+			        n_iter, res, 0.0, 0,
+			        rad2deg(get_actual_val(min->x, priv.initial, priv.rv, 0)),
+			        rad2deg(get_actual_val(min->x, priv.initial, priv.rv, 1)),
+			        get_actual_val(min->x, priv.initial, priv.rv, 2),
+			        get_actual_val(min->x, priv.initial, priv.rv, 3)*1e10);
+		}
 
-	} while ( i < 30 && !done );
+		status = gsl_multimin_test_size(min->size, 0.005);
+
+	} while ( status == GSL_CONTINUE && n_iter < 1000 );
 
 	if ( verbose ) {
-		double asx, asy, asz;
-		double bsx, bsy, bsz;
-		double csx, csy, csz;
-		cell_get_reciprocal(cell, &asx, &asy, &asz,
-		                          &bsx, &bsy, &bsz,
-		                          &csx, &csy, &csz);
-		STATUS("Final asx = %e\n", asx);
+		STATUS("Done with refinement after %i iter\n", n_iter);
+		STATUS("status = %i (%s)\n", status, gsl_strerror(status));
 	}
-}
 
+	if ( status == GSL_SUCCESS ) {
 
-static struct prdata pr_refine(Crystal *cr, const RefList *full,
-                               PartialityModel pmodel)
-{
-	int verbose = 0;
-	struct prdata prdata;
+		if ( check_angle_shifts(min->x, priv.initial, priv.rv, &priv) ) return;
 
-	prdata.refined = 0;
-	prdata.n_filtered = 0;
+		if ( verbose ) {
 
-	if ( verbose ) {
-		write_residual_graph(cr, full);
+			double res = residual_f(min->x, &priv);
+			double size = gsl_multimin_fminimizer_size(min);
+			STATUS("At end: %f %f %f %f ----> %f %f %e %f residual = %e size %f\n",
+			       gsl_vector_get(min->x, 0),
+			       gsl_vector_get(min->x, 1),
+			       gsl_vector_get(min->x, 2),
+			       gsl_vector_get(min->x, 3),
+			       rad2deg(get_actual_val(min->x, priv.initial, priv.rv, 0)),
+			       rad2deg(get_actual_val(min->x, priv.initial, priv.rv, 1)),
+			       get_actual_val(min->x, priv.initial, priv.rv, 2),
+			       get_actual_val(min->x, priv.initial, priv.rv, 3)*1e10,
+			       res, size);
+
+		}
+
+		if ( fh != NULL ) {
+			double res = residual_f(min->x, &priv);
+			fprintf(fh, "%5i %10.8f  %10.8f  %5i  %10.8f %10.8f  %e  %e\n",
+			        n_iter, res, 0.0, 0,
+			        rad2deg(get_actual_val(min->x, priv.initial, priv.rv, 0)),
+			        rad2deg(get_actual_val(min->x, priv.initial, priv.rv, 1)),
+			        get_actual_val(min->x, priv.initial, priv.rv, 2),
+			        get_actual_val(min->x, priv.initial, priv.rv, 3)*1e10);
+			fclose(fh);
+		}
+
+		/* Apply the final shifts */
+		apply_parameters(min->x, priv.initial, priv.rv, cr);
+		update_predictions(cr);
+		calculate_partialities(cr, PMODEL_XSPHERE);
+		r = linear_scale(full, crystal_get_reflections(cr), &G, 0);
+		if ( r == 0 ) {
+			crystal_set_osf(cr, G);
+		}
+
+		if ( verbose ) {
+
+			STATUS("After applying final shifts:\n");
+			STATUS("PR final: dev = %10.5e, free dev = %10.5e\n",
+			       residual(cr, full, 0, NULL, NULL, 0),
+			       residual(cr, full, 1, NULL, NULL, 0));
+			STATUS("Final R = %e m^-1\n", crystal_get_profile_radius(cr));
+
+		}
+
+	} else {
+		ERROR("Bad refinement: crystal %i\n", serial);
 	}
 
-	do_pr_refine(cr, full, pmodel, verbose);
+	if ( write_logs ) {
+		write_gridscan(cr, full, cycle, serial);
+		write_specgraph(cr, full, cycle, serial);
+		write_test_logs(cr, full, cycle, serial);
+	}
 
-	if ( crystal_get_user_flag(cr) == 0 ) {
-		prdata.refined = 1;
+	if ( crystal_get_profile_radius(cr) > 5e9 ) {
+		ERROR("Very large radius: crystal %i\n", serial);
 	}
 
-	return prdata;
+	gsl_multimin_fminimizer_free(min);
+	gsl_vector_free(priv.initial);
+	gsl_vector_free(priv.vals);
+	gsl_vector_free(priv.step);
 }
 
 
@@ -706,7 +943,13 @@ struct refine_args
 	RefList *full;
 	Crystal *crystal;
 	PartialityModel pmodel;
+	int serial;
 	struct prdata prdata;
+	int verbose;
+	int cycle;
+	int no_logs;
+	SymOpList *sym;
+	SymOpList *amb;
 };
 
 
@@ -724,8 +967,18 @@ static void refine_image(void *task, int id)
 {
 	struct refine_args *pargs = task;
 	Crystal *cr = pargs->crystal;
+	int write_logs = 0;
+
+	write_logs = !pargs->no_logs && (pargs->serial % 20 == 0);
+	pargs->prdata.refined = 0;
 
-	pargs->prdata = pr_refine(cr, pargs->full, pargs->pmodel);
+	do_pr_refine(cr, pargs->full, pargs->pmodel, pargs->verbose,
+	             pargs->serial, pargs->cycle, write_logs,
+	             pargs->sym, pargs->amb);
+
+	if ( crystal_get_user_flag(cr) == 0 ) {
+		pargs->prdata.refined = 1;
+	}
 }
 
 
@@ -738,6 +991,7 @@ static void *get_image(void *vqargs)
 	memcpy(task, &qargs->task_defaults, sizeof(struct refine_args));
 
 	task->crystal = qargs->crystals[qargs->n_started];
+	task->serial = qargs->n_started;
 
 	qargs->n_started++;
 
@@ -757,7 +1011,9 @@ static void done_image(void *vqargs, void *task)
 
 
 void refine_all(Crystal **crystals, int n_crystals,
-                RefList *full, int nthreads, PartialityModel pmodel)
+                RefList *full, int nthreads, PartialityModel pmodel,
+                int verbose, int cycle, int no_logs,
+                SymOpList *sym, SymOpList *amb)
 {
 	struct refine_args task_defaults;
 	struct queue_args qargs;
@@ -766,7 +1022,11 @@ void refine_all(Crystal **crystals, int n_crystals,
 	task_defaults.crystal = NULL;
 	task_defaults.pmodel = pmodel;
 	task_defaults.prdata.refined = 0;
-	task_defaults.prdata.n_filtered = 0;
+	task_defaults.verbose = verbose;
+	task_defaults.cycle = cycle;
+	task_defaults.no_logs = no_logs;
+	task_defaults.sym = sym;
+	task_defaults.amb = amb;
 
 	qargs.task_defaults = task_defaults;
 	qargs.n_started = 0;
diff --git a/src/post-refinement.h b/src/post-refinement.h
index 7f395013..e2416a23 100644
--- a/src/post-refinement.h
+++ b/src/post-refinement.h
@@ -3,11 +3,11 @@
  *
  * Post refinement
  *
- * Copyright © 2012-2015 Deutsches Elektronen-Synchrotron DESY,
+ * Copyright © 2012-2018 Deutsches Elektronen-Synchrotron DESY,
  *                       a research centre of the Helmholtz Association.
  *
  * Authors:
- *   2010-2015 Thomas White <taw@physics.org>
+ *   2010-2018 Thomas White <taw@physics.org>
  *
  * This file is part of CrystFEL.
  *
@@ -41,6 +41,7 @@
 #include "utils.h"
 #include "crystal.h"
 #include "geometry.h"
+#include "symmetry.h"
 
 
 enum prflag
@@ -51,19 +52,31 @@ enum prflag
 	PRFLAG_EARLY = 18,
 	PRFLAG_CC = 19,
 	PRFLAG_BIGB = 20,
+	PRFLAG_SCALEBAD = 21,
 };
 
 
 extern const char *str_prflag(enum prflag flag);
 
 extern void refine_all(Crystal **crystals, int n_crystals,
-                       RefList *full, int nthreads, PartialityModel pmodel);
+                       RefList *full, int nthreads, PartialityModel pmodel,
+                       int verbose, int cycle, int no_logs,
+                       SymOpList *sym, SymOpList *amb);
+
+extern void write_gridscan(Crystal *cr, const RefList *full,
+                           int cycle, int serial);
+
+extern void write_specgraph(Crystal *crystal, const RefList *full,
+                            signed int cycle, int serial);
 
 /* Exported so it can be poked by tests/pr_p_gradient_check */
 extern double gradient(Crystal *cr, int k, Reflection *refl,
                        PartialityModel pmodel);
 
 extern double residual(Crystal *cr, const RefList *full, int free,
-                       int *pn_used, const char *filename);
+                       int *pn_used, const char *filename, int complain);
+
+extern void write_test_logs(Crystal *crystal, const RefList *full,
+                            signed int cycle, int serial);
 
 #endif	/* POST_REFINEMENT_H */
diff --git a/src/process_image.c b/src/process_image.c
index 87089289..dbd4f427 100644
--- a/src/process_image.c
+++ b/src/process_image.c
@@ -246,7 +246,7 @@ void process_image(const struct index_args *iargs, struct pattern_args *pargs,
 	if ( iargs->fix_bandwidth >= 0.0 ) {
 		image.bw = iargs->fix_bandwidth;
 	} else {
-		image.bw = 0.00000001;
+		image.bw = 0.0013;
 	}
 
 	if ( image_feature_count(image.features) < iargs->min_peaks ) {
@@ -301,7 +301,7 @@ void process_image(const struct index_args *iargs, struct pattern_args *pargs,
 	/* Integrate! */
 	time_accounts_set(taccs, TACC_INTEGRATION);
 	sb_shared->pings[cookie]++;
-	integrate_all_5(&image, iargs->int_meth, PMODEL_SCSPHERE,
+	integrate_all_5(&image, iargs->int_meth, PMODEL_XSPHERE,
 	                iargs->push_res,
 	                iargs->ir_inn, iargs->ir_mid, iargs->ir_out,
 	                iargs->int_diag, iargs->int_diag_h,
diff --git a/src/rejection.c b/src/rejection.c
index 10f13f9c..9d41f9b4 100644
--- a/src/rejection.c
+++ b/src/rejection.c
@@ -137,7 +137,7 @@ static void check_cc(Crystal *cr, RefList *full)
 		pcalc = get_partiality(refl);
 
 		/* Observed partiality */
-		corr = exp(-G) * exp(B*res*res) * get_lorentz(refl);
+		corr = G * exp(B*res*res) * get_lorentz(refl);
 		Ipart = get_intensity(refl) * corr;
 		pobs = Ipart / get_intensity(match);
 
diff --git a/src/scaling.c b/src/scaling.c
index 77a64384..f4c93c66 100644
--- a/src/scaling.c
+++ b/src/scaling.c
@@ -3,11 +3,11 @@
  *
  * Scaling
  *
- * Copyright © 2012-2015 Deutsches Elektronen-Synchrotron DESY,
+ * Copyright © 2012-2017 Deutsches Elektronen-Synchrotron DESY,
  *                       a research centre of the Helmholtz Association.
  *
  * Authors:
- *   2010-2015 Thomas White <taw@physics.org>
+ *   2010-2017 Thomas White <taw@physics.org>
  *
  * This file is part of CrystFEL.
  *
@@ -38,16 +38,14 @@
 #include <gsl/gsl_linalg.h>
 #include <gsl/gsl_eigen.h>
 #include <gsl/gsl_blas.h>
+#include <gsl/gsl_fit.h>
 
 #include "merge.h"
 #include "post-refinement.h"
 #include "symmetry.h"
 #include "cell.h"
 #include "cell-utils.h"
-
-
-/* Maximum number of iterations of NLSq to do for each image per macrocycle. */
-#define MAX_CYCLES (10)
+#include "scaling.h"
 
 
 /* Apply the given shift to the 'k'th parameter of 'image'. */
@@ -64,9 +62,9 @@ static void apply_shift(Crystal *cr, int k, double shift)
 		break;
 
 		case GPARAM_OSF :
-		t = crystal_get_osf(cr);
+		t = -log(crystal_get_osf(cr));
 		t += shift;
-		crystal_set_osf(cr, t);
+		crystal_set_osf(cr, exp(-t));
 		break;
 
 		default :
@@ -78,8 +76,7 @@ static void apply_shift(Crystal *cr, int k, double shift)
 
 
 /* Perform one cycle of scaling of 'cr' against 'full' */
-static double scale_iterate(Crystal *cr, const RefList *full,
-                            PartialityModel pmodel, int *nr)
+static double scale_iterate(Crystal *cr, const RefList *full, int *nr)
 {
 	gsl_matrix *M;
 	gsl_vector *v;
@@ -182,7 +179,7 @@ static double scale_iterate(Crystal *cr, const RefList *full,
 
 			}
 
-			fx = G + log(p) - log(L) - B*s*s + log(I_full);
+			fx = -log(G) + log(p) - log(L) - B*s*s + log(I_full);
 			delta_I = log(I_partial) - fx;
 			v_c = w * delta_I * gradients[k];
 			v_curr = gsl_vector_get(v, k);
@@ -271,7 +268,7 @@ double log_residual(Crystal *cr, const RefList *full, int free,
 		if ( I_full <= 0 ) continue;  /* Because log */
 		if ( p <= 0.0 ) continue; /* Because of log */
 
-		fx = G + log(p) - log(L) - B*s*s + log(I_full);
+		fx = -log(G) + log(p) - log(L) - B*s*s + log(I_full);
 		dc = log(I_partial) - fx;
 		w = 1.0;
 		dev += w*dc*dc;
@@ -290,8 +287,7 @@ double log_residual(Crystal *cr, const RefList *full, int free,
 }
 
 
-static void do_scale_refine(Crystal *cr, const RefList *full,
-                            PartialityModel pmodel, int *nr)
+static void do_scale_refine(Crystal *cr, const RefList *full, int *nr)
 {
 	int i, done;
 	double old_dev;
@@ -304,7 +300,7 @@ static void do_scale_refine(Crystal *cr, const RefList *full,
 
 		double dev;
 
-		scale_iterate(cr, full, pmodel, nr);
+		scale_iterate(cr, full, nr);
 
 		dev = log_residual(cr, full, 0, 0, NULL);
 		if ( fabs(dev - old_dev) < dev*0.01 ) done = 1;
@@ -312,7 +308,7 @@ static void do_scale_refine(Crystal *cr, const RefList *full,
 		i++;
 		old_dev = dev;
 
-	} while ( i < MAX_CYCLES && !done );
+	} while ( i < 10 && !done );
 }
 
 
@@ -320,7 +316,6 @@ struct scale_args
 {
 	RefList *full;
 	Crystal *crystal;
-	PartialityModel pmodel;
 	int n_reflections;
 };
 
@@ -339,8 +334,7 @@ struct queue_args
 static void scale_crystal(void *task, int id)
 {
 	struct scale_args *pargs = task;
-	do_scale_refine(pargs->crystal, pargs->full, pargs->pmodel,
-	                &pargs->n_reflections);
+	do_scale_refine(pargs->crystal, pargs->full, &pargs->n_reflections);
 }
 
 
@@ -381,14 +375,12 @@ static double total_log_r(Crystal **crystals, int n_crystals, RefList *full,
 	int n = 0;
 
 	for ( i=0; i<n_crystals; i++ ) {
-
 		double r;
 		if ( crystal_get_user_flag(crystals[i]) ) continue;
 		r = log_residual(crystals[i], full, 0, NULL, NULL);
 		if ( isnan(r) ) continue;
 		total += r;
 		n++;
-
 	}
 
 	if ( ninc != NULL ) *ninc = n;
@@ -397,8 +389,7 @@ static double total_log_r(Crystal **crystals, int n_crystals, RefList *full,
 
 
 /* Perform iterative scaling, all the way to convergence */
-void scale_all(Crystal **crystals, int n_crystals, int nthreads,
-               PartialityModel pmodel)
+void scale_all(Crystal **crystals, int n_crystals, int nthreads)
 {
 	struct scale_args task_defaults;
 	struct queue_args qargs;
@@ -406,7 +397,6 @@ void scale_all(Crystal **crystals, int n_crystals, int nthreads,
 	int niter = 0;
 
 	task_defaults.crystal = NULL;
-	task_defaults.pmodel = pmodel;
 
 	qargs.task_defaults = task_defaults;
 	qargs.n_crystals = n_crystals;
@@ -422,7 +412,7 @@ void scale_all(Crystal **crystals, int n_crystals, int nthreads,
 		double bef_res;
 
 		full = merge_intensities(crystals, n_crystals, nthreads,
-		                         pmodel, 2, INFINITY, 0);
+		                         2, INFINITY, 0);
 		old_res = new_res;
 		bef_res = total_log_r(crystals, n_crystals, full, NULL);
 
@@ -456,3 +446,178 @@ void scale_all(Crystal **crystals, int n_crystals, int nthreads,
 		ERROR("Too many iterations - giving up!\n");
 	}
 }
+
+
+static void scale_crystal_linear(void *task, int id)
+{
+	struct scale_args *pargs = task;
+	int r;
+	double G;
+
+	/* Simple iterative algorithm */
+	r = linear_scale(pargs->full, crystal_get_reflections(pargs->crystal), &G, 0);
+	if ( r == 0 ) {
+		crystal_set_osf(pargs->crystal, G);
+	} else {
+		crystal_set_user_flag(pargs->crystal, PRFLAG_SCALEBAD);
+	}
+}
+
+
+/* Calculates G, by which list2 should be multiplied to fit list1 */
+int linear_scale(const RefList *list1, const RefList *list2, double *G,
+                 int complain_loudly)
+{
+	const Reflection *refl1;
+	const Reflection *refl2;
+	RefListIterator *iter;
+	int max_n = 256;
+	int n = 0;
+	double *x;
+	double *y;
+	double *w;
+	int r;
+	double cov11;
+	double sumsq;
+	int n_esd1 = 0;
+	int n_esd2 = 0;
+	int n_ih1 = 0;
+	int n_ih2 = 0;
+	int n_nan1 = 0;
+	int n_nan2 = 0;
+	int n_inf1 = 0;
+	int n_inf2 = 0;
+	int n_part = 0;
+	int n_nom = 0;
+
+	x = malloc(max_n*sizeof(double));
+	w = malloc(max_n*sizeof(double));
+	y = malloc(max_n*sizeof(double));
+	if ( (x==NULL) || (y==NULL) || (w==NULL) ) {
+		ERROR("Failed to allocate memory for scaling.\n");
+		return 1;
+	}
+
+	int nb = 0;
+	for ( refl2 = first_refl_const(list2, &iter);
+	      refl2 != NULL;
+	      refl2 = next_refl_const(refl2, iter) )
+	{
+		signed int h, k, l;
+		double Ih1, Ih2;
+		double esd1, esd2;
+		nb++;
+
+		get_indices(refl2, &h, &k, &l);
+		refl1 = find_refl(list1, h, k, l);
+		if ( refl1 == NULL ) {
+			n_nom++;
+			continue;
+		}
+
+		Ih1 = get_intensity(refl1);
+		Ih2 = get_intensity(refl2);
+		esd1 = get_esd_intensity(refl1);
+		esd2 = get_esd_intensity(refl2);
+
+		/* Problem cases in approximate descending order of severity */
+		if ( isnan(Ih1) ) { n_nan1++; continue; }
+		if ( isinf(Ih1) ) { n_inf1++; continue; }
+		if ( isnan(Ih2) ) { n_nan2++; continue; }
+		if ( isinf(Ih2) ) { n_inf2++; continue; }
+		if ( get_partiality(refl2) < 0.3 ) { n_part++; continue; }
+		//if ( Ih1 <= 0.0 ) { n_ih1++; continue; }
+		if ( Ih2 <= 0.0 ) { n_ih2++; continue; }
+		//if ( Ih1 <= 3.0*esd1 ) { n_esd1++;  continue; }
+		if ( Ih2 <= 3.0*esd2 ) { n_esd2++; continue; }
+		//if ( get_redundancy(refl1) < 2 ) continue;
+
+		if ( n == max_n ) {
+			max_n *= 2;
+			x = realloc(x, max_n*sizeof(double));
+			y = realloc(y, max_n*sizeof(double));
+			w = realloc(w, max_n*sizeof(double));
+			if ( (x==NULL) || (y==NULL) || (w==NULL) ) {
+				ERROR("Failed to allocate memory for scaling.\n");
+				return 1;
+			}
+		}
+
+		x[n] = Ih2 / get_partiality(refl2);
+		y[n] = Ih1;
+		w[n] = get_partiality(refl2);
+		n++;
+
+	}
+
+	if ( n < 2 ) {
+		if ( complain_loudly ) {
+			ERROR("Not enough reflections for scaling (had %i, but %i remain)\n", nb, n);
+			if ( n_esd1 ) ERROR("%i reference reflection esd\n", n_esd1);
+			if ( n_esd2 ) ERROR("%i subject reflection esd\n", n_esd2);
+			if ( n_ih1 ) ERROR("%i reference reflection intensity\n", n_ih1);
+			if ( n_ih2 ) ERROR("%i subject reflection intensity\n", n_ih2);
+			if ( n_nan1 ) ERROR("%i reference reflection nan\n", n_nan1);
+			if ( n_nan2 ) ERROR("%i subject reflection nan\n", n_nan2);
+			if ( n_inf1 ) ERROR("%i reference reflection inf\n", n_inf1);
+			if ( n_inf2 ) ERROR("%i subject reflection inf\n", n_inf2);
+			if ( n_part ) ERROR("%i subject reflection partiality\n", n_part);
+			if ( n_nom ) ERROR("%i no match in reference list\n", n_nom);
+		}
+		*G = 1.0;
+		return 1;
+	}
+
+	r = gsl_fit_wmul(x, 1, w, 1, y, 1, n, G, &cov11, &sumsq);
+
+	if ( r ) {
+		ERROR("Scaling failed.\n");
+		return 1;
+	}
+
+	if ( isnan(*G) ) {
+
+		if ( complain_loudly ) {
+			ERROR("Scaling gave NaN (%i pairs)\n", n);
+			if ( n < 10 ) {
+				int i;
+				for ( i=0; i<n; i++ ) {
+					STATUS("%i %e %e %e\n", i, x[i], y[i], w[n]);
+					}
+			}
+		}
+
+		*G = 1.0;
+		return 1;
+	}
+
+	free(x);
+	free(y);
+	free(w);
+
+	return 0;
+}
+
+
+void scale_all_to_reference(Crystal **crystals, int n_crystals,
+                            RefList *reference, int nthreads)
+{
+	struct scale_args task_defaults;
+	struct queue_args qargs;
+
+	task_defaults.crystal = NULL;
+
+	qargs.task_defaults = task_defaults;
+	qargs.n_crystals = n_crystals;
+	qargs.crystals = crystals;
+
+	/* Don't have threads which are doing nothing */
+	if ( n_crystals < nthreads ) nthreads = n_crystals;
+
+	qargs.task_defaults.full = reference;
+	qargs.n_started = 0;
+	qargs.n_done = 0;
+	qargs.n_reflections = 0;
+	run_threads(nthreads, scale_crystal_linear, get_crystal, done_crystal,
+	            &qargs, n_crystals, 0, 0, 0);
+}
diff --git a/src/scaling.h b/src/scaling.h
index 534043d9..514bee03 100644
--- a/src/scaling.h
+++ b/src/scaling.h
@@ -3,11 +3,11 @@
  *
  * Scaling
  *
- * Copyright © 2012-2015 Deutsches Elektronen-Synchrotron DESY,
+ * Copyright © 2012-2017 Deutsches Elektronen-Synchrotron DESY,
  *                       a research centre of the Helmholtz Association.
  *
  * Authors:
- *   2010-2015 Thomas White <taw@physics.org>
+ *   2010-2017 Thomas White <taw@physics.org>
  *
  * This file is part of CrystFEL.
  *
@@ -41,7 +41,12 @@
 extern double log_residual(Crystal *cr, const RefList *full, int free,
                            int *pn_used, const char *filename);
 
-extern void scale_all(Crystal **crystals, int n_crystals, int nthreads,
-                      PartialityModel pmodel);
+extern int linear_scale(const RefList *list1, const RefList *list2, double *G,
+                        int complain_loudly);
+
+extern void scale_all(Crystal **crystals, int n_crystals, int nthreads);
+
+extern void scale_all_to_reference(Crystal **crystals, int n_crystals,
+                                   RefList *reference, int nthreads);
 
 #endif	/* SCALING_H */