New partiality model from Ginn et al.

12 files changed, 324 insertions, 287 deletions

diff --git a/libcrystfel/src/geometry.c b/libcrystfel/src/geometry.c
index 576655de..beae9a16 100644
--- a/libcrystfel/src/geometry.c
+++ b/libcrystfel/src/geometry.c
@@ -238,76 +238,37 @@ static double random_partiality(signed int h, signed int k, signed int l,
 }
 
 
-static double partiality(PartialityModel pmodel,
-                         signed int h, signed int k, signed int l,
-                         int serial,
-                         double rlow, double rhigh, double pr)
-{
-	double D = rlow - rhigh;
-
-	/* Convert to partiality */
-	switch ( pmodel ) {
-
-		default:
-		case PMODEL_UNITY:
-		return 1.0;
-
-		case PMODEL_SCSPHERE:
-		return 4.0*sphere_fraction(rlow, rhigh, pr)*pr / (3.0*D);
-
-		case PMODEL_SCGAUSSIAN:
-		return 4.0*gaussian_fraction(rlow, rhigh, pr)*pr / (3.0*D);
-
-		case PMODEL_RANDOM:
-		return random_partiality(h, k, l, serial);
-
-	}
-}
-
-
 static Reflection *check_reflection(struct image *image, Crystal *cryst,
                                     signed int h, signed int k, signed int l,
                                     double xl, double yl, double zl,
                                     Reflection *updateme)
 {
-	double tl;
-	double rlow, rhigh;    /* "Excitation error" */
-	double klow, khigh;    /* Wavenumber */
 	Reflection *refl;
-	double cet, cez;  /* Centre of Ewald sphere */
-	double pr;
-	double del;
+	double R, top;
+	double kmin, kmax, k0, knom, k1;
+	double dcs, exerr;
 
 	/* Don't predict 000 */
 	if ( (updateme == NULL) && (abs(h)+abs(k)+abs(l) == 0) ) return NULL;
 
-	pr = crystal_get_profile_radius(cryst);
-	del = image->div + crystal_get_mosaicity(cryst);
-
-	/* "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);
-
-	/* If the point is looking "backscattery", reject it straight away */
-	if ( (updateme == NULL) && (zl < -khigh/2.0) ) return NULL;
-
-	tl = sqrt(xl*xl + yl*yl);
-
-	cet = -sin(del/2.0) * khigh;
-	cez = -cos(del/2.0) * khigh;
-	rhigh = khigh - distance(cet, cez, tl, zl);  /* Loss of precision */
-
-	cet =  sin(del/2.0) * klow;
-	cez = -cos(del/2.0) * klow;
-	rlow = klow - distance(cet, cez, tl, zl);  /* Loss of precision */
-
-	/* Condition for reflection to be excited at all */
-	if ( (updateme == NULL)
-	     && (signbit(rlow) == signbit(rhigh))
-	     && (fabs(rlow) > pr)
-	     && (fabs(rhigh) > pr) ) return NULL;
+	/* Calculate the limiting wavelengths, lambda0 and lambda1
+	 * = 1/k0 and 1/k1 respectively */
+	R = crystal_get_profile_radius(cryst);
+	top = R*R - xl*xl - yl*yl - zl*zl;
+	k0 = top/(2.0*(zl+R));
+	k1 = top/(2.0*(zl-R));
+
+	/* The reflection is excited if any of the reflection is within 2sigma
+	 * of the nominal * wavelength of the X-ray beam
+	 * (NB image->bw is full width) */
+	kmin = 1.0/(image->lambda + image->lambda*image->bw);
+	knom = 1.0/image->lambda;
+	kmax = 1.0/(image->lambda - image->lambda*image->bw);
+	if ( (k1>kmax) || (k0<kmin) ) return NULL;
+
+	/* Calculate excitation error */
+	dcs = distance3d(0.0, 0.0, -knom, xl, yl, zl);
+	exerr = 1.0/image->lambda - dcs;  /* Loss of precision */
 
 	if ( updateme == NULL ) {
 		refl = reflection_new(h, k, l);
@@ -321,7 +282,7 @@ static Reflection *check_reflection(struct image *image, Crystal *cryst,
 	if ( (image->det != NULL) && (updateme != NULL) ) {
 
 		double fs, ss;
-		locate_peak_on_panel(xl, yl, zl, 1.0/image->lambda,
+		locate_peak_on_panel(xl, yl, zl, knom,
 		                     get_panel(updateme), &fs, &ss);
 		set_detector_pos(refl, fs, ss);
 
@@ -333,7 +294,7 @@ static Reflection *check_reflection(struct image *image, Crystal *cryst,
 
 		double fs, ss;        /* Position on detector */
 		signed int p;         /* Panel number */
-		p = locate_peak(xl, yl, zl, 1.0/image->lambda,
+		p = locate_peak(xl, yl, zl, knom,
 		                image->det, &fs, &ss);
 		if ( p == -1 ) {
 			reflection_free(refl);
@@ -344,20 +305,8 @@ static Reflection *check_reflection(struct image *image, Crystal *cryst,
 
 	}
 
-	if ( unlikely(rlow < rhigh) ) {
-		ERROR("Reflection with rlow < rhigh!\n");
-		ERROR("%3i %3i %3i  rlow = %e, rhigh = %e\n",
-		      h, k, l, rlow, rhigh);
-		ERROR("div + m = %e, R = %e, bw = %e\n", del, pr, image->bw);
-		/* If we are updating, this is (kind of) OK */
-		if ( updateme == NULL ) {
-			reflection_free(refl);
-			return NULL;
-		}
-	}
-
-	set_partial(refl, rlow, rhigh, 1.0); /* Actual partiality set by
-	                                      * calculate_partialities() */
+	set_kpred(refl, knom);
+	set_exerr(refl, exerr);
 	set_lorentz(refl, 1.0);
 	set_symmetric_indices(refl, h, k, l);
 	set_redundancy(refl, 1);
@@ -368,18 +317,12 @@ static Reflection *check_reflection(struct image *image, Crystal *cryst,
 
 double r_gradient(UnitCell *cell, int k, Reflection *refl, struct image *image)
 {
-	double azi;
 	double asx, asy, asz;
 	double bsx, bsy, bsz;
 	double csx, csy, csz;
 	double xl, yl, zl;
 	signed int hs, ks, ls;
-	double rlow, rhigh, p;
-	double philow, phihigh, phi;
-	double khigh, klow;
-	double tl, cet, cez;
-
-	get_partial(refl, &rlow, &rhigh, &p);
+	double tl, phi, azi;
 
 	get_symmetric_indices(refl, &hs, &ks, &ls);
 
@@ -390,26 +333,9 @@ double r_gradient(UnitCell *cell, int k, Reflection *refl, struct image *image)
 	yl = hs*asy + ks*bsy + ls*csy;
 	zl = hs*asz + ks*bsz + ls*csz;
 
-	/* "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);
-
 	tl = sqrt(xl*xl + yl*yl);
-
-	cet = -sin(image->div/2.0) * klow;
-	cez = -cos(image->div/2.0) * klow;
-	philow = angle_between_2d(tl-cet, zl-cez, 0.0, 1.0);
-
-	cet = -sin(image->div/2.0) * khigh;
-	cez = -cos(image->div/2.0) * khigh;
-	phihigh = angle_between_2d(tl-cet, zl-cez, 0.0, 1.0);
-
-	/* Approximation: philow and phihigh are very similar */
-	phi = (philow + phihigh) / 2.0;
-
-	azi = atan2(yl, xl);
+	phi = angle_between_2d(tl, zl+1.0/image->lambda, 0.0, 1.0); /* 2theta */
+	azi = atan2(yl, xl); /* azimuth */
 
 	switch ( k ) {
 
@@ -531,11 +457,17 @@ RefList *predict_to_res(Crystal *cryst, double max_res)
 }
 
 
-static void set_unity_partialities(RefList *list)
+static void set_unity_partialities(Crystal *cryst)
 {
+	RefList *list;
 	Reflection *refl;
 	RefListIterator *iter;
 
+	list = crystal_get_reflections(cryst);
+	if ( list == NULL ) {
+		ERROR("No reflections for partiality calculation!\n");
+		return;
+	}
 	for ( refl = first_refl(list, &iter);
 	      refl != NULL;
 	      refl = next_refl(refl, iter) )
@@ -546,25 +478,11 @@ static void set_unity_partialities(RefList *list)
 }
 
 
-/**
- * calculate_partialities:
- * @cryst: A %Crystal
- * @pmodel: A %PartialityModel
- *
- * Calculates the partialities for the reflections in @cryst, given the current
- * crystal and image parameters.  The crystal's image and reflection lists
- * must be set.  The specified %PartialityModel will be used.
- *
- * You must not have changed the crystal or image parameters since you last
- * called predict_to_res() or update_predictions(), because this function
- * relies on the limiting wavelength values calculated by those functions.
- */
-void calculate_partialities(Crystal *cryst, PartialityModel pmodel)
+static void set_random_partialities(Crystal *cryst)
 {
 	RefList *list;
 	Reflection *refl;
 	RefListIterator *iter;
-	double pr;
 	struct image *image;
 
 	list = crystal_get_reflections(cryst);
@@ -573,8 +491,109 @@ void calculate_partialities(Crystal *cryst, PartialityModel pmodel)
 		return;
 	}
 
-	if ( pmodel == PMODEL_UNITY ) {
-		set_unity_partialities(list);
+	image = crystal_get_image(cryst);
+	if ( image == NULL ) {
+		ERROR("No image structure for partiality calculation!\n");
+		return;
+	}
+
+	for ( refl = first_refl(list, &iter);
+	      refl != NULL;
+	      refl = next_refl(refl, iter) )
+	{
+		signed int h, k, l;
+		get_symmetric_indices(refl, &h, &k, &l);
+		set_partiality(refl, random_partiality(h, k, l, image->serial));
+		set_lorentz(refl, 1.0);
+	}
+}
+
+
+static double do_integral(double q2, double zl, double R,
+                          double lambda, double sig, int verbose)
+{
+	int i;
+	double kmin, kmax, kstart, kfinis;
+	double inc;
+	double total = 0.0;
+	double k0, k1;
+	const int SAMPLES = 50;  /* Number of samples for integration */
+	const double N = 1.5;  /* Pointiness of spectrum */
+	FILE *fh = NULL;
+
+	k0 = (R*R - q2)/(2.0*(zl+R));
+	k1 = (R*R - q2)/(2.0*(zl-R));
+
+	/* Range over which E is significantly different from zero */
+	kmin = 1.0 / (lambda + 5.0*sig);
+	kmax = 1.0 / (lambda - 5.0*sig);
+
+	kstart = kmin > k1 ? kmin : k1;
+	kfinis = (k0 < 0.0) || (kmax < k0) ? kmax : k0;
+	inc = (kfinis - kstart) / SAMPLES;
+
+	if ( verbose ) {
+		char fn[64];
+		snprintf(fn, 63, "partial%i.graph", verbose);
+		fh = fopen(fn, "w");
+		fprintf(fh, "  n    p      wavelength   E           P\n");
+		STATUS("Nominal k = %e m^-1\n", 1.0/lambda);
+		STATUS(" (wavelength %e m)\n", lambda);
+		STATUS("Bandwidth %e m\n", sig);
+		STATUS("k1/2 = %e m^-1\n", -q2/(2.0*zl));
+		STATUS(" (wavelength %e m)\n", 1.0/(-q2/(2.0*zl)));
+		STATUS("Reflection k goes from %e to %e m^-1\n", k1, k0);
+		STATUS(" (wavelengths from %e to %e m\n", 1.0/k1, 1.0/k0);
+		STATUS("Beam goes from %e to %e m^-1\n", kmin, kmax);
+		STATUS(" (wavelengths from %e to %e m\n", 1.0/kmin, 1.0/kmax);
+		STATUS("Integration goes from %e to %e m^-1\n", kstart, kfinis);
+		STATUS(" (wavelengths from %e to %e m\n", 1.0/kstart, 1.0/kfinis);
+	}
+
+	for ( i=0; i<SAMPLES; i++ ) {
+
+		double p, kp, lrel;
+		double E, P;
+
+		kp = kstart + i*inc;
+		double pref = sqrt(q2 + kp*kp + 2.0*zl*kp)/(2.0*R);
+		p = pref + 0.5 - kp/(2.0*R);
+
+		/* Spectral energy term */
+		lrel = fabs(1.0/kp - lambda);
+		E = exp(-0.5 * pow(lrel / sig, N));
+		E /= sqrt(2.0 * M_PI * sig);
+
+		/* RLP profile term */
+		P = 4.0*p * (1.0 - p);
+
+		total += E*P*inc;
+
+		if ( fh != NULL ) {
+			fprintf(fh, "%3i %f %e %e %e\n", i, p, 1.0/kp, E, P);
+		}
+	}
+
+	if ( fh != NULL ) fclose(fh);
+
+	return total;
+}
+
+
+
+static void ginn_spectrum_partialities(Crystal *cryst)
+{
+	RefList *list;
+	Reflection *refl;
+	RefListIterator *iter;
+	double r0, m, lambda, sig;
+	struct image *image;
+	UnitCell *cell;
+	double asx, asy, asz, bsx, bsy, bsz, csx, csy, csz;
+
+	list = crystal_get_reflections(cryst);
+	if ( list == NULL ) {
+		ERROR("No reflections for partiality calculation!\n");
 		return;
 	}
 
@@ -584,31 +603,89 @@ void calculate_partialities(Crystal *cryst, PartialityModel pmodel)
 		return;
 	}
 
-	pr = crystal_get_profile_radius(cryst);
+	cell = crystal_get_cell(cryst);
+	if ( cell == NULL ) {
+		ERROR("No unit cell for partiality calculation!\n");
+		return;
+	}
+	cell_get_reciprocal(cell, &asx, &asy, &asz,
+	                          &bsx, &bsy, &bsz,
+	                          &csx, &csy, &csz);
+
+	r0 = crystal_get_profile_radius(cryst);
+	m = crystal_get_mosaicity(cryst);
+	lambda = image->lambda;
+	sig = image->bw * lambda;
 
 	for ( refl = first_refl(crystal_get_reflections(cryst), &iter);
 	      refl != NULL;
 	      refl = next_refl(refl, iter) )
 	{
-		double rlow, rhigh, part;
+		double R;
 		signed int h, k, l;
+		double xl, yl, zl;
+		double q2;
+		double total, norm;
 
 		get_symmetric_indices(refl, &h, &k, &l);
-		get_partial(refl, &rlow, &rhigh, &part);
+		xl = h*asx + k*bsx + l*csx;
+		yl = h*asy + k*bsy + l*csy;
+		zl = h*asz + k*bsz + l*csz;
 
-		part = partiality(pmodel, h, k, l, image->serial,
-		                  rlow, rhigh, pr);
+		/* Radius of rlp profile */
+		q2 = xl*xl + yl*yl + zl*zl;
 
-		if ( isnan(part) && ((h!=0) || (k!=0) || (l!=0)) ) {
-			ERROR("Assigning NAN partiality!\n");
-			ERROR("%3i %3i %3i  rlow = %e, rhigh = %e\n",
-			      h, k, l, rlow, rhigh);
-			ERROR("R = %e, bw = %e\n", pr, image->bw);
-			ERROR("D = %e\n", rlow - rhigh);
+		R = r0 + m * sqrt(q2);
+
+		total = do_integral(q2, zl, R, lambda, sig, 0);
+		norm = do_integral(q2, -0.5*q2*lambda, R, lambda, sig, 0);
+
+	        set_partiality(refl, total/norm);
+		set_lorentz(refl, 1.0);
+
+		if ( total > 2.0*norm ) {
+			/* Error! */
+			do_integral(q2, zl, R, lambda, sig, 1);
+			do_integral(q2, -0.5*q2*lambda, R, lambda, sig, 2);
 			abort();
-	        }
+		}
 
-	        set_partiality(refl, part);
+	}
+}
+
+
+/**
+ * calculate_partialities:
+ * @cryst: A %Crystal
+ * @pmodel: A %PartialityModel
+ *
+ * Calculates the partialities for the reflections in @cryst, given the current
+ * crystal and image parameters.  The crystal's image and reflection lists
+ * must be set.  The specified %PartialityModel will be used.
+ *
+ * You must not have changed the crystal or image parameters since you last
+ * called predict_to_res() or update_predictions(), because this function
+ * relies on the limiting wavelength values calculated by those functions.
+ */
+void calculate_partialities(Crystal *cryst, PartialityModel pmodel)
+{
+	switch ( pmodel ) {
+
+		case PMODEL_UNITY :
+		set_unity_partialities(cryst);
+		break;
+
+		case PMODEL_XSPHERE :
+		ginn_spectrum_partialities(cryst);
+		break;
+
+		case PMODEL_RANDOM :
+		set_random_partialities(cryst);
+		break;
+
+		default :
+		ERROR("Unknown partiality model %i\n", pmodel);
+		break;
 
 	}
 }
@@ -694,29 +771,30 @@ void polarisation_correction(RefList *list, UnitCell *cell, struct image *image)
 
 
 /* Returns dx_h/dP, where P = any parameter */
-double x_gradient(int param, Reflection *refl, UnitCell *cell,
-                  struct panel *p, double lambda)
+double x_gradient(int param, Reflection *refl, UnitCell *cell, struct panel *p)
 {
 	signed int h, k, l;
-	double x, z, wn;
-	double ax, ay, az, bx, by, bz, cx, cy, cz;
+	double xl, zl, kpred;
+	double asx, asy, asz, bsx, bsy, bsz, csx, csy, csz;
 
 	get_indices(refl, &h, &k, &l);
-	wn = 1.0 / lambda;
-	cell_get_reciprocal(cell, &ax, &ay, &az, &bx, &by, &bz, &cx, &cy, &cz);
-	x = h*ax + k*bx + l*cx;
-	z = h*az + k*bz + l*cz;
+	kpred = get_kpred(refl);
+	cell_get_reciprocal(cell, &asx, &asy, &asz,
+	                          &bsx, &bsy, &bsz,
+	                          &csx, &csy, &csz);
+	xl = h*asx + k*bsx + l*csx;
+	zl = h*asz + k*bsz + l*csz;
 
 	switch ( param ) {
 
 		case GPARAM_ASX :
-		return h * p->clen / (wn+z);
+		return h * p->clen / (kpred + zl);
 
 		case GPARAM_BSX :
-		return k * p->clen / (wn+z);
+		return k * p->clen / (kpred + zl);
 
 		case GPARAM_CSX :
-		return l * p->clen / (wn+z);
+		return l * p->clen / (kpred + zl);
 
 		case GPARAM_ASY :
 		return 0.0;
@@ -728,13 +806,13 @@ double x_gradient(int param, Reflection *refl, UnitCell *cell,
 		return 0.0;
 
 		case GPARAM_ASZ :
-		return -h * x * p->clen / (wn*wn + 2*wn*z + z*z);
+		return -h * xl * p->clen / (kpred*kpred + 2.0*kpred*zl + zl*zl);
 
 		case GPARAM_BSZ :
-		return -k * x * p->clen / (wn*wn + 2*wn*z + z*z);
+		return -k * xl * p->clen / (kpred*kpred + 2.0*kpred*zl + zl*zl);
 
 		case GPARAM_CSZ :
-		return -l * x * p->clen / (wn*wn + 2*wn*z + z*z);
+		return -l * xl * p->clen / (kpred*kpred + 2.0*kpred*zl + zl*zl);
 
 		case GPARAM_DETX :
 		return -1;
@@ -743,7 +821,7 @@ double x_gradient(int param, Reflection *refl, UnitCell *cell,
 		return 0;
 
 		case GPARAM_CLEN :
-		return x / (wn+z);
+		return xl / (kpred+zl);
 
 	}
 
@@ -753,18 +831,19 @@ double x_gradient(int param, Reflection *refl, UnitCell *cell,
 
 
 /* Returns dy_h/dP, where P = any parameter */
-double y_gradient(int param, Reflection *refl, UnitCell *cell,
-                  struct panel *p, double lambda)
+double y_gradient(int param, Reflection *refl, UnitCell *cell, struct panel *p)
 {
 	signed int h, k, l;
-	double y, z, wn;
-	double ax, ay, az, bx, by, bz, cx, cy, cz;
+	double yl, zl, kpred;
+	double asx, asy, asz, bsx, bsy, bsz, csx, csy, csz;
 
 	get_indices(refl, &h, &k, &l);
-	wn = 1.0 / lambda;
-	cell_get_reciprocal(cell, &ax, &ay, &az, &bx, &by, &bz, &cx, &cy, &cz);
-	y = h*ay + k*by + l*cy;
-	z = h*az + k*bz + l*cz;
+	kpred = get_kpred(refl);
+	cell_get_reciprocal(cell, &asx, &asy, &asz,
+	                          &bsx, &bsy, &bsz,
+	                          &csx, &csy, &csz);
+	yl = h*asy + k*bsy + l*csy;
+	zl = h*asz + k*bsz + l*csz;
 
 	switch ( param ) {
 
@@ -778,22 +857,22 @@ double y_gradient(int param, Reflection *refl, UnitCell *cell,
 		return 0.0;
 
 		case GPARAM_ASY :
-		return h * p->clen / (wn+z);
+		return h * p->clen / (kpred + zl);
 
 		case GPARAM_BSY :
-		return k * p->clen / (wn+z);
+		return k * p->clen / (kpred + zl);
 
 		case GPARAM_CSY :
-		return l * p->clen / (wn+z);
+		return l * p->clen / (kpred + zl);
 
 		case GPARAM_ASZ :
-		return -h * y * p->clen / (wn*wn + 2*wn*z + z*z);
+		return -h * yl * p->clen / (kpred*kpred + 2.0*kpred*zl + zl*zl);
 
 		case GPARAM_BSZ :
-		return -k * y * p->clen / (wn*wn + 2*wn*z + z*z);
+		return -k * yl * p->clen / (kpred*kpred + 2.0*kpred*zl + zl*zl);
 
 		case GPARAM_CSZ :
-		return -l * y * p->clen / (wn*wn + 2*wn*z + z*z);
+		return -l * yl * p->clen / (kpred*kpred + 2.0*kpred*zl + zl*zl);
 
 		case GPARAM_DETX :
 		return 0;
@@ -802,7 +881,7 @@ double y_gradient(int param, Reflection *refl, UnitCell *cell,
 		return -1;
 
 		case GPARAM_CLEN :
-		return y / (wn+z);
+		return yl / (kpred+zl);
 
 	}
 
diff --git a/libcrystfel/src/geometry.h b/libcrystfel/src/geometry.h
index 4fbd0bca..c2f2899d 100644
--- a/libcrystfel/src/geometry.h
+++ b/libcrystfel/src/geometry.h
@@ -3,7 +3,7 @@
  *
  * Geometry of diffraction
  *
- * Copyright © 2013-2016 Deutsches Elektronen-Synchrotron DESY,
+ * Copyright © 2013-2017 Deutsches Elektronen-Synchrotron DESY,
  *                       a research centre of the Helmholtz Association.
  * Copyright © 2012 Richard Kirian
  *
@@ -47,8 +47,7 @@ extern "C" {
 /**
  * PartialityModel:
  * @PMODEL_UNITY : Set all all partialities and Lorentz factors to 1.
- * @PMODEL_SCSPHERE : Sphere model with source coverage factor included
- * @PMODEL_SCGAUSSIAN : Gaussian model with source coverage factor included
+ * @PMODEL_XSPHERE : Flat sphere model with super-Gaussian spectrum
  * @PMODEL_RANDOM : Randomly assigned partialities
  *
  * A %PartialityModel describes a geometrical model which can be used to
@@ -57,8 +56,7 @@ extern "C" {
 typedef enum {
 
 	PMODEL_UNITY,
-	PMODEL_SCSPHERE,
-	PMODEL_SCGAUSSIAN,
+	PMODEL_XSPHERE,
 	PMODEL_RANDOM,
 
 } PartialityModel;
@@ -99,9 +97,9 @@ extern double sphere_fraction(double rlow, double rhigh, double pr);
 extern double gaussian_fraction(double rlow, double rhigh, double pr);
 
 extern double x_gradient(int param, Reflection *refl, UnitCell *cell,
-                         struct panel *p, double lambda);
+                         struct panel *p);
 extern double y_gradient(int param, Reflection *refl, UnitCell *cell,
-                         struct panel *p, double lambda);
+                         struct panel *p);
 
 #ifdef __cplusplus
 }
diff --git a/libcrystfel/src/integration.c b/libcrystfel/src/integration.c
index c18ae110..fcbdf658 100644
--- a/libcrystfel/src/integration.c
+++ b/libcrystfel/src/integration.c
@@ -1760,7 +1760,7 @@ void integrate_all_2(struct image *image, IntegrationMethod meth,
                      IntDiag int_diag,
                      signed int idh, signed int idk, signed int idl)
 {
-	integrate_all_3(image, meth, PMODEL_SCSPHERE, 0.0,
+	integrate_all_3(image, meth, PMODEL_XSPHERE, 0.0,
 	                ir_inn, ir_mid, ir_out,
 	                int_diag, idh, idk, idl);
 }
diff --git a/libcrystfel/src/predict-refine.c b/libcrystfel/src/predict-refine.c
index 3246dbec..14a60bc7 100644
--- a/libcrystfel/src/predict-refine.c
+++ b/libcrystfel/src/predict-refine.c
@@ -91,9 +91,7 @@ static void twod_mapping(double fs, double ss, double *px, double *py,
 static double r_dev(struct reflpeak *rp)
 {
 	/* Excitation error term */
-	double rlow, rhigh, p;
-	get_partial(rp->refl, &rlow, &rhigh, &p);
-	return (rlow+rhigh)/2.0;
+	return get_exerr(rp->refl);
 }
 
 
@@ -425,7 +423,7 @@ static int iterate(struct reflpeak *rps, int n, UnitCell *cell,
 		/* Positional x terms */
 		for ( k=0; k<num_params; k++ ) {
 			gradients[k] = x_gradient(rv[k], rps[i].refl, cell,
-			                          rps[i].panel, image->lambda);
+			                          rps[i].panel);
 		}
 
 		for ( k=0; k<num_params; k++ ) {
@@ -457,7 +455,7 @@ static int iterate(struct reflpeak *rps, int n, UnitCell *cell,
 		/* Positional y terms */
 		for ( k=0; k<num_params; k++ ) {
 			gradients[k] = y_gradient(rv[k], rps[i].refl, cell,
-			                          rps[i].panel, image->lambda);
+			                          rps[i].panel);
 		}
 
 		for ( k=0; k<num_params; k++ ) {
diff --git a/libcrystfel/src/reflist.c b/libcrystfel/src/reflist.c
index 707c802e..69fb29be 100644
--- a/libcrystfel/src/reflist.c
+++ b/libcrystfel/src/reflist.c
@@ -67,10 +67,10 @@ struct _refldata {
 	signed int ls;
 
 	/* Partiality and related geometrical stuff */
-	double rlow;  /* Low excitation error */
-	double rhigh;  /* High excitation error */
-	double p;   /* Partiality */
-	double L;   /* Lorentz factor */
+	double kpred; /* Wavenumber of middle of reflection */
+	double exerr; /* Excitation error */
+	double p;     /* Partiality */
+	double L;     /* Lorentz factor */
 
 	/* Location in image */
 	double fs;
@@ -422,22 +422,28 @@ double get_intensity(const Reflection *refl)
 
 
 /**
- * get_partial:
+ * get_kpred:
  * @refl: A %Reflection
- * @rlow: Location at which to store the "low" excitation error
- * @rhigh: Location at which to store the "high" excitation error
- * @p: Location at which to store the partiality
  *
- * This function is used during post refinement (in conjunction with
- * set_partial()) to get access to the details of the partiality calculation.
+ * Returns: the wavenumber at the centre of the reflection
  *
  **/
-void get_partial(const Reflection *refl, double *rlow, double *rhigh,
-                 double *p)
+double get_kpred(const Reflection *refl)
 {
-	*rlow = refl->data.rlow;
-	*rhigh = refl->data.rhigh;
-	*p = get_partiality(refl);
+	return refl->data.kpred;
+}
+
+
+/**
+ * get_exerr:
+ * @refl: A %Reflection
+ *
+ * Returns: the excitation error (in m^-1) for this reflection
+ *
+ **/
+double get_exerr(const Reflection *refl)
+{
+	return refl->data.exerr;
 }
 
 
@@ -614,21 +620,28 @@ void set_panel(Reflection *refl, struct panel *p)
 
 
 /**
- * set_partial:
+ * set_kpred:
  * @refl: A %Reflection
- * @rlow: The "low" excitation error
- * @rhigh: The "high" excitation error
- * @p: The partiality
+ * @kpred: The wavenumber at which the reflection should be predicted
  *
- * This function is used during post refinement (in conjunction with
- * get_partial()) to get access to the details of the partiality calculation.
+ * Sets the wavenumber at the centre of the reflection.  Used by
+ * predict_to_res() and update_predictions()
+ **/
+void set_kpred(Reflection *refl, double kpred)
+{
+	refl->data.kpred = kpred;
+}
+
+
+/**
+ * set_exerr:
+ * @refl: A %Reflection
+ * @exerr: The excitation error for the reflection
  *
  **/
-void set_partial(Reflection *refl, double rlow, double rhigh, double p)
+void set_exerr(Reflection *refl, double exerr)
 {
-	refl->data.rlow = rlow;
-	refl->data.rhigh = rhigh;
-	refl->data.p = p;
+	refl->data.exerr = exerr;
 }
 
 
diff --git a/libcrystfel/src/reflist.h b/libcrystfel/src/reflist.h
index f2126ac9..4dd1e56f 100644
--- a/libcrystfel/src/reflist.h
+++ b/libcrystfel/src/reflist.h
@@ -87,6 +87,8 @@ extern Reflection *next_found_refl(Reflection *refl);
 extern void get_detector_pos(const Reflection *refl, double *fs, double *ss);
 extern struct panel *get_panel(const Reflection *refl);
 extern double get_partiality(const Reflection *refl);
+extern double get_kpred(const Reflection *refl);
+extern double get_exerr(const Reflection *refl);
 extern double get_lorentz(const Reflection *refl);
 extern void get_indices(const Reflection *refl,
                         signed int *h, signed int *k, signed int *l);
@@ -94,8 +96,6 @@ extern void get_symmetric_indices(const Reflection *refl,
                                   signed int *hs, signed int *ks,
                                   signed int *ls);
 extern double get_intensity(const Reflection *refl);
-extern void get_partial(const Reflection *refl, double *rlow, double *rhigh,
-                        double *p);
 extern int get_redundancy(const Reflection *refl);
 extern double get_temp1(const Reflection *refl);
 extern double get_temp2(const Reflection *refl);
@@ -109,7 +109,8 @@ extern int get_flag(const Reflection *refl);
 extern void copy_data(Reflection *to, const Reflection *from);
 extern void set_detector_pos(Reflection *refl, double fs, double ss);
 extern void set_panel(Reflection *refl, struct panel *p);
-extern void set_partial(Reflection *refl, double rlow, double rhigh, double p);
+extern void set_kpred(Reflection *refl, double kpred);
+extern void set_exerr(Reflection *refl, double exerr);
 extern void set_partiality(Reflection *refl, double p);
 extern void set_lorentz(Reflection *refl, double L);
 extern void set_intensity(Reflection *refl, double intensity);
diff --git a/src/partial_sim.c b/src/partial_sim.c
index 4a348742..b2581096 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;
diff --git a/src/partialator.c b/src/partialator.c
index 428b0d68..5ad106c7 100644
--- a/src/partialator.c
+++ b/src/partialator.c
@@ -770,7 +770,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;
@@ -975,10 +975,8 @@ int main(int argc, char *argv[])
 	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 {
diff --git a/src/post-refinement.c b/src/post-refinement.c
index bf30d299..74075503 100644
--- a/src/post-refinement.c
+++ b/src/post-refinement.c
@@ -85,20 +85,6 @@ const char *str_prflag(enum prflag flag)
 }
 
 
-/* Returns dp(gauss)/dr at "r" */
-static double gaussian_fraction_gradient(double r, double R)
-{
-	const double ng = 2.6;
-	const double sig = R/ng;
-
-	/* 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;
-
-	return exp(-pow(r/sig, 2.0)/2.0) / (sig*sqrt(2.0*M_PI));
-}
-
-
 /* Returns dp(sph)/dr at "r" */
 static double sphere_fraction_gradient(double r, double pr)
 {
@@ -130,14 +116,10 @@ static double partiality_gradient(double r, double pr,
 		case PMODEL_UNITY:
 		return 0.0;
 
-		case PMODEL_SCSPHERE:
+		case PMODEL_XSPHERE:
 		A = sphere_fraction_gradient(r, pr)/D;
 		return 4.0*pr*A/3.0;
 
-		case PMODEL_SCGAUSSIAN:
-		A = gaussian_fraction_gradient(r, pr)/D;
-		return 4.0*pr*A/3.0;
-
 	}
 }
 
@@ -152,19 +134,6 @@ static double sphere_fraction_rgradient(double r, double R)
 }
 
 
-static double gaussian_fraction_rgradient(double r, double R)
-{
-	const double ng = 2.6;
-	const double sig = R/ng;
-
-	/* 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;
-
-	return -(ng*r/(sqrt(2.0*M_PI)*R*R))*exp(-r*r/(2.0*sig*sig));
-}
-
-
 static double volume_fraction_rgradient(double r, double pr,
                                        PartialityModel pmodel)
 {
@@ -173,12 +142,9 @@ static double volume_fraction_rgradient(double r, double pr,
 		case PMODEL_UNITY :
 		return 1.0;
 
-		case PMODEL_SCSPHERE :
+		case PMODEL_XSPHERE :
 		return sphere_fraction_rgradient(r, pr);
 
-		case PMODEL_SCGAUSSIAN :
-		return gaussian_fraction_rgradient(r, pr);
-
 		default :
 		ERROR("No pmodel in volume_fraction_rgradient!\n");
 		return 1.0;
@@ -194,12 +160,9 @@ static double volume_fraction(double rlow, double rhigh, double pr,
 		case PMODEL_UNITY :
 		return 1.0;
 
-		case PMODEL_SCSPHERE :
+		case PMODEL_XSPHERE :
 		return sphere_fraction(rlow, rhigh, pr);
 
-		case PMODEL_SCGAUSSIAN :
-		return gaussian_fraction(rlow, rhigh, pr);
-
 		default :
 		ERROR("No pmodel in volume_fraction!\n");
 		return 1.0;
@@ -212,12 +175,14 @@ static double volume_fraction(double rlow, double rhigh, double pr,
 double gradient(Crystal *cr, int k, Reflection *refl, PartialityModel pmodel)
 {
 	double glow, ghigh;
-	double rlow, rhigh, p;
+	double rlow, rhigh;
 	struct image *image = crystal_get_image(cr);
 	double R = crystal_get_profile_radius(cr);
 	double gr;
 
-	get_partial(refl, &rlow, &rhigh, &p);
+	/* FIXME ! */
+	rlow = 0.0;
+	rhigh = 0.0;
 
 	if ( k == GPARAM_R ) {
 
@@ -607,7 +572,7 @@ static void write_residual_graph(Crystal *cr, const RefList *full)
 		                          bsx, bsy, bsz,
 		                          csx, csy, csz);
 		update_predictions(cr);
-		calculate_partialities(cr, PMODEL_SCSPHERE);
+		calculate_partialities(cr, PMODEL_XSPHERE);
 		res = residual(cr, full, 0, &n, NULL);
 		fprintf(fh, "%i %e %e %i\n", i, asx, res, n);
 	}
@@ -616,7 +581,7 @@ static void write_residual_graph(Crystal *cr, const RefList *full)
 	                          bsx, bsy, bsz,
 	                          csx, csy, csz);
 	update_predictions(cr);
-	calculate_partialities(cr, PMODEL_SCSPHERE);
+	calculate_partialities(cr, PMODEL_XSPHERE);
 	fclose(fh);
 }
 
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/tests/pr_p_gradient_check.c b/tests/pr_p_gradient_check.c
index 5322fcca..b9da05bc 100644
--- a/tests/pr_p_gradient_check.c
+++ b/tests/pr_p_gradient_check.c
@@ -59,7 +59,7 @@ static void scan_partialities(RefList *reflections, RefList *compare,
 	{
 		signed int h, k, l;
 		Reflection *refl2;
-		double rlow, rhigh, p;
+		double p;
 
 		get_indices(refl, &h, &k, &l);
 		refl2 = find_refl(compare, h, k, l);
@@ -69,7 +69,7 @@ static void scan_partialities(RefList *reflections, RefList *compare,
 			continue;
 		}
 
-		get_partial(refl2, &rlow, &rhigh, &p);
+		p = get_partiality(refl2);
 		vals[idx][i] = p;
 		if ( unlikely(p < 0.0) ) {
 			ERROR("Negative partiality! %3i %3i %3i  %f\n",
@@ -306,8 +306,6 @@ static double test_gradients(Crystal *cr, double incr_val, int refine,
 			i++;
 		} else {
 
-			double r1, r2, p;
-
 			grad1 = (vals[1][i] - vals[0][i]) / incr_val;
 			grad2 = (vals[2][i] - vals[1][i]) / incr_val;
 			grad = (grad1 + grad2) / 2.0;
@@ -315,8 +313,6 @@ static double test_gradients(Crystal *cr, double incr_val, int refine,
 
 			cgrad = gradient(cr, refine, refl, pmodel);
 
-			get_partial(refl, &r1, &r2, &p);
-
 			if ( isnan(cgrad) ) {
 				n_nan++;
 				continue;
@@ -343,11 +339,11 @@ static double test_gradients(Crystal *cr, double incr_val, int refine,
 			{
 
 				if ( !quiet ) {
-					STATUS("!- %s %3i %3i %3i"
-					       " %10.2Le %10.2e ratio = %5.2Lf"
-					       " %10.2e %10.2e\n",
+					STATUS("!- %s %3i %3i %3i "
+					       "%10.2Le %10.2e "
+                                               "ratio = %5.2Lf\n",
 					       str, h, k, l, grad, cgrad,
-					       cgrad/grad, r1, r2);
+					       cgrad/grad);
 				}
 				n_bad++;
 
@@ -447,18 +443,15 @@ int main(int argc, char *argv[])
 
 	rng = gsl_rng_alloc(gsl_rng_mt19937);
 
-	for ( i=0; i<2; i++ ) {
+	for ( i=0; i<1; i++ ) {
 
 		UnitCell *rot;
 		double val;
 		PartialityModel pmodel;
 
 		if ( i == 0 ) {
-			pmodel = PMODEL_SCSPHERE;
-			STATUS("Testing SCSphere model:\n");
-		} else if ( i == 1 ) {
-			pmodel = PMODEL_SCGAUSSIAN;
-			STATUS("Testing SCGaussian model.\n");
+			pmodel = PMODEL_XSPHERE;
+			STATUS("Testing XSphere model:\n");
 		} else {
 			ERROR("WTF?\n");
 			return 1;
diff --git a/tests/prediction_gradient_check.c b/tests/prediction_gradient_check.c
index 0ffc07ca..85d61ead 100644
--- a/tests/prediction_gradient_check.c
+++ b/tests/prediction_gradient_check.c
@@ -73,7 +73,6 @@ static void scan(RefList *reflections, RefList *compare,
 	{
 		signed int h, k, l;
 		Reflection *refl2;
-		double rlow, rhigh, p;
 		double fs, ss, xh, yh;
 		struct panel *panel;
 
@@ -85,7 +84,6 @@ static void scan(RefList *reflections, RefList *compare,
 			continue;
 		}
 
-		get_partial(refl2, &rlow, &rhigh, &p);
 		get_detector_pos(refl2, &fs, &ss);
 		panel = get_panel(refl2);
 		twod_mapping(fs, ss, &xh, &yh, panel);
@@ -93,7 +91,7 @@ static void scan(RefList *reflections, RefList *compare,
 		switch ( checkrxy ) {
 
 			case 0 :
-			vals[idx][i] = (rlow + rhigh)/2.0;
+			vals[idx][i] = get_exerr(refl2);
 			break;
 
 			case 1 :
@@ -279,8 +277,6 @@ static double test_gradients(Crystal *cr, double incr_val, int refine,
 			i++;
 		} else {
 
-			double r1, r2, p;
-
 			grad1 = (vals[1][i] - vals[0][i]) / incr_val;
 			grad2 = (vals[2][i] - vals[1][i]) / incr_val;
 			grad = (grad1 + grad2) / 2.0;
@@ -300,18 +296,14 @@ static double test_gradients(Crystal *cr, double incr_val, int refine,
 				if ( checkrxy == 1 ) {
 					cgrad = x_gradient(refine, refl,
 					         crystal_get_cell(cr),
-					         &image->det->panels[0],
-					         crystal_get_image(cr)->lambda);
+					         &image->det->panels[0]);
 				} else {
 					cgrad = y_gradient(refine, refl,
 					         crystal_get_cell(cr),
-					         &image->det->panels[0],
-					         crystal_get_image(cr)->lambda);
+					         &image->det->panels[0]);
 				}
 			}
 
-			get_partial(refl, &r1, &r2, &p);
-
 			if ( isnan(cgrad) ) {
 				n_nan++;
 				continue;
@@ -338,11 +330,11 @@ static double test_gradients(Crystal *cr, double incr_val, int refine,
 			{
 
 				if ( !quiet ) {
-					STATUS("!- %s %3i %3i %3i"
-					       " %10.2Le %10.2e ratio = %5.2Lf"
-					       " %10.2e %10.2e\n",
+					STATUS("!- %s %3i %3i %3i "
+					       "%10.2Le %10.2e "
+					       "ratio = %5.2Lf\n",
 					       str, h, k, l, grad, cgrad,
-					       cgrad/grad, r1, r2);
+					       cgrad/grad);
 				}
 				n_bad++;