Remove "sphere", "thin" and "gaussian" partiality models, add "scgaussian"

13 files changed, 45 insertions, 1095 deletions

diff --git a/Makefile.am b/Makefile.am
index 033ea746..78f3e277 100644
--- a/Makefile.am
+++ b/Makefile.am
@@ -10,8 +10,7 @@ bin_PROGRAMS = src/pattern_sim src/process_hkl src/get_hkl src/indexamajig \
 noinst_PROGRAMS = tests/list_check tests/integration_check \
                   tests/pr_p_gradient_check tests/symmetry_check \
                   tests/centering_check tests/transformation_check \
-                  tests/cell_check tests/pr_l_gradient_check \
-                  tests/pr_pl_gradient_check tests/ring_check \
+                  tests/cell_check tests/ring_check \
                   tests/prof2d_check tests/ambi_check
 
 MERGE_CHECKS = tests/first_merge_check tests/second_merge_check \
@@ -20,9 +19,7 @@ MERGE_CHECKS = tests/first_merge_check tests/second_merge_check \
 PARTIAL_CHECKS = tests/partialator_merge_check_1 \
                  tests/partialator_merge_check_2 \
                  tests/partialator_merge_check_3 \
-                 tests/pr_p_gradient_check \
-                 tests/pr_l_gradient_check \
-                 tests/pr_pl_gradient_check
+                 tests/pr_p_gradient_check
 
 TESTS = tests/list_check $(MERGE_CHECKS) $(PARTIAL_CHECKS) \
         tests/integration_check \
@@ -117,12 +114,6 @@ tests_ambi_check_SOURCES = tests/ambi_check.c
 tests_pr_p_gradient_check_SOURCES = tests/pr_p_gradient_check.c \
                                   src/post-refinement.c
 
-tests_pr_pl_gradient_check_SOURCES = tests/pr_pl_gradient_check.c \
-                                  src/post-refinement.c
-
-tests_pr_l_gradient_check_SOURCES = tests/pr_l_gradient_check.c \
-                                  src/post-refinement.c
-
 tests_centering_check_SOURCES = tests/centering_check.c
 
 tests_transformation_check_SOURCES = tests/transformation_check.c
diff --git a/libcrystfel/src/geometry.c b/libcrystfel/src/geometry.c
index aff19198..8952486d 100644
--- a/libcrystfel/src/geometry.c
+++ b/libcrystfel/src/geometry.c
@@ -121,18 +121,10 @@ static double partiality(PartialityModel pmodel,
 		case PMODEL_UNITY:
 		return 1.0;
 
-		case PMODEL_SPHERE:
-		plow  = 3.0*qlow*qlow  - 2.0*qlow*qlow*qlow;
-		phigh = 3.0*qhigh*qhigh - 2.0*qhigh*qhigh*qhigh;
-		return plow - phigh;
-
-		case PMODEL_GAUSSIAN:
+		case PMODEL_SCGAUSSIAN:
 		plow = 0.5 * gsl_sf_erf(ng * rlow / (sqrt(2.0)*r));
 		phigh = 0.5 * gsl_sf_erf(ng * rhigh / (sqrt(2.0)*r));
-		return plow - phigh;
-
-		case PMODEL_THIN:
-		return 1.0 - (rmid*rmid)/(r*r);
+		return 4.0*(plow-phigh)*r / (3.0*D);
 
 		case PMODEL_SCSPHERE:
 		plow  = 3.0*qlow*qlow  - 2.0*qlow*qlow*qlow;
@@ -157,7 +149,7 @@ static Reflection *check_reflection(struct image *image, Crystal *cryst,
 	Reflection *refl;
 	double cet, cez;  /* Centre of Ewald sphere */
 	double pr;
-	double L, D;
+	double D;
 	double del;
 
 	/* Don't predict 000 */
@@ -198,44 +190,13 @@ static Reflection *check_reflection(struct image *image, Crystal *cryst,
 		return NULL;
 	}
 
-	/* Conditions for reflection to be excited at all */
-	switch ( pmodel ) {
-
-		case PMODEL_UNITY:  /* PMODEL_UNITY shouldn't end up here */
-		case PMODEL_SPHERE:
-		case PMODEL_GAUSSIAN:
-		case PMODEL_SCSPHERE:
-		if ( (signbit(rlow) == signbit(rhigh))
-		     && (fabs(rlow) > pr)
-		     && (fabs(rhigh) > pr) ) return NULL;
-		break;
-
-		case PMODEL_THIN:
-		if ( fabs(rmid) > pr ) return NULL;
-		break;
-
-	}
+	/* Condition for reflection to be excited at all */
+	if ( (signbit(rlow) == signbit(rhigh))
+	     && (fabs(rlow) > pr)
+	     && (fabs(rhigh) > pr) ) return NULL;
 
 	D = rlow - rhigh;
 
-	/* Lorentz factor is determined direction from the r values, before
-	 * clamping.  The multiplication by 0.01e9 to make the
-	 * correction factor vaguely near 1. */
-	switch ( pmodel ) {
-
-		case PMODEL_SPHERE:
-		case PMODEL_GAUSSIAN:
-		L = LORENTZ_SCALE / D;
-		break;
-
-		case PMODEL_UNITY:  /* PMODEL_UNITY shouldn't end up here */
-		case PMODEL_THIN:
-		case PMODEL_SCSPHERE:
-		L = 1.0;
-		break;
-
-	}
-
 	/* If the "lower" Ewald sphere is a long way away, use the
 	 * position at which the Ewald sphere would just touch the
 	 * reflection.
@@ -282,19 +243,8 @@ static Reflection *check_reflection(struct image *image, Crystal *cryst,
 		set_detector_pos(refl, 0.0, xda, yda);
 	}
 
-	if ( pmodel != PMODEL_THIN ) {
-		set_partial(refl, rlow, rhigh, part, clamp_low, clamp_high);
-	} else {
-		/* If we are using the TES (Thin Ewald Sphere) model, we abuse
-		 * the fields as follows:
-		 *  rlow       = the r value for the middle (only) Ewald sphere
-		 *  rhigh      = 0.0
-		 *  clamp_low  =  0
-		 *  clamp_high = +1
-		 */
-		set_partial(refl, rmid, 0.0, part, 0, +1);
-	}
-	set_lorentz(refl, L);
+	set_partial(refl, rlow, rhigh, part, clamp_low, clamp_high);
+	set_lorentz(refl, 1.0);
 	set_symmetric_indices(refl, h, k, l);
 	set_redundancy(refl, 1);
 
diff --git a/libcrystfel/src/geometry.h b/libcrystfel/src/geometry.h
index d8d226f0..ccdea1c2 100644
--- a/libcrystfel/src/geometry.h
+++ b/libcrystfel/src/geometry.h
@@ -46,23 +46,18 @@ extern "C" {
 
 /**
  * PartialityModel:
- * @PMODEL_SPHERE : Intersection of sphere with excited volume of reciprocal
- *   space.
  * @PMODEL_UNITY : Set all all partialities and Lorentz factors to 1.
- * @PMODEL_GAUSSIAN : Gaussian profiles in 3D
- * @PMODEL_THIN : Thin Ewald sphere intersecting sphere
  * @PMODEL_SCSPHERE : Sphere model with source coverage factor included
+ * @PMODEL_SCGAUSSIAN : Gaussian model with source coverage factor included
  *
  * A %PartialityModel describes a geometrical model which can be used to
  * calculate spot partialities and Lorentz correction factors.
  **/
 typedef enum {
 
-	PMODEL_SPHERE,
 	PMODEL_UNITY,
-	PMODEL_GAUSSIAN,
-	PMODEL_THIN,
 	PMODEL_SCSPHERE,
+	PMODEL_SCGAUSSIAN,
 
 } PartialityModel;
 
diff --git a/libcrystfel/src/integration.c b/libcrystfel/src/integration.c
index 635a7c74..34bde83a 100644
--- a/libcrystfel/src/integration.c
+++ b/libcrystfel/src/integration.c
@@ -1914,7 +1914,8 @@ 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_SPHERE, 0.0, ir_inn, ir_mid, ir_out,
+	integrate_all_3(image, meth, PMODEL_SCSPHERE, 0.0,
+	                ir_inn, ir_mid, ir_out,
 	                int_diag, idh, idk, idl);
 }
 
diff --git a/src/partial_sim.c b/src/partial_sim.c
index 10f5e58d..77b5a73b 100644
--- a/src/partial_sim.c
+++ b/src/partial_sim.c
@@ -357,7 +357,7 @@ static void run_job(void *vwargs, int cookie)
 		snprintf(wargs->image.filename, 255, "dummy.h5");
 	}
 
-	reflections = find_intersections(&wargs->image, cr, PMODEL_SPHERE);
+	reflections = find_intersections(&wargs->image, cr, PMODEL_SCSPHERE);
 	crystal_set_reflections(cr, reflections);
 
 	for ( i=0; i<NBINS; i++ ) {
diff --git a/src/partialator.c b/src/partialator.c
index 0d532997..bdee29c5 100644
--- a/src/partialator.c
+++ b/src/partialator.c
@@ -219,7 +219,7 @@ int main(int argc, char *argv[])
 	Stream *st;
 	Crystal **crystals;
 	char *pmodel_str = NULL;
-	PartialityModel pmodel = PMODEL_SPHERE;
+	PartialityModel pmodel = PMODEL_SCSPHERE;
 	int min_measurements = 2;
 	char *rval;
 	struct srdata srdata;
@@ -342,14 +342,10 @@ int main(int argc, char *argv[])
 	free(sym_str);
 
 	if ( pmodel_str != NULL ) {
-		if ( strcmp(pmodel_str, "sphere") == 0 ) {
-			pmodel = PMODEL_SPHERE;
-		} else if ( strcmp(pmodel_str, "unity") == 0 ) {
+		if ( strcmp(pmodel_str, "unity") == 0 ) {
 			pmodel = PMODEL_UNITY;
-		} else if ( strcmp(pmodel_str, "gaussian") == 0 ) {
-			pmodel = PMODEL_GAUSSIAN;
-		} else if ( strcmp(pmodel_str, "thin") == 0 ) {
-			pmodel = PMODEL_THIN;
+		} else if ( strcmp(pmodel_str, "scgaussian") == 0 ) {
+			pmodel = PMODEL_SCGAUSSIAN;
 		} else if ( strcmp(pmodel_str, "scsphere") == 0 ) {
 			pmodel = PMODEL_SCSPHERE;
 		} else {
diff --git a/src/post-refinement.c b/src/post-refinement.c
index 4c16f2d8..e971f96a 100644
--- a/src/post-refinement.c
+++ b/src/post-refinement.c
@@ -83,24 +83,17 @@ static double partiality_gradient(double r, double profile_radius,
 		case PMODEL_UNITY:
 		return 0.0;
 
-		case PMODEL_SPHERE:
-		dqdr = 1.0 / (2.0*profile_radius);
-		return dpdq(r, profile_radius) * dqdr;
-
-		case PMODEL_GAUSSIAN:
-		/* FIXME: Get a proper gradient */
-		dqdr = 1.0 / (2.0*profile_radius);
-		return dpdq(r, profile_radius) * dqdr;
-
-		case PMODEL_THIN:
-		return -(2.0*r)/(profile_radius*profile_radius);
-
 		case PMODEL_SCSPHERE:
 		dqdr = 1.0 / (2.0*profile_radius);
 		dpsphdr = dpdq(r, profile_radius) * dqdr;
 		A = (dpsphdr/D) - p*pow(rlow-rhigh, -2.0);
 		return 4.0*profile_radius*A/3.0;
 
+		case PMODEL_SCGAUSSIAN:
+		/* FIXME: Get a proper gradient */
+		dqdr = 1.0 / (2.0*profile_radius);
+		return dpdq(r, profile_radius) * dqdr;
+
 	}
 }
 
@@ -117,18 +110,16 @@ static double partiality_rgradient(double r, double profile_radius,
 		case PMODEL_UNITY:
 		return 0.0;
 
-		case PMODEL_SPHERE:
+		case PMODEL_SCSPHERE:
+		/* FIXME: wrong (?) */
 		dqdrad = -0.5 * r / (profile_radius * profile_radius);
 		return dpdq(r, profile_radius) * dqdrad;
 
-		case PMODEL_GAUSSIAN:
+		case PMODEL_SCGAUSSIAN:
 		/* FIXME: Get a proper gradient */
 		dqdrad = -0.5 * r / (profile_radius * profile_radius);
 		return dpdq(r, profile_radius) * dqdrad;
 
-		case PMODEL_THIN:
-		return 2.0*r*r*pow(profile_radius, -3.0);
-
 	}
 }
 
@@ -206,18 +197,16 @@ double p_gradient(Crystal *cr, int k, Reflection *refl, PartialityModel pmodel)
 			case PMODEL_UNITY:
 			return 0.0;
 
-			case PMODEL_GAUSSIAN:
+			case PMODEL_SCSPHERE:
 			gr  = partiality_rgradient(rlow, r, pmodel);
 			gr -= partiality_rgradient(rhigh, r, pmodel);
 			return gr;
 
-			case PMODEL_SPHERE:
+			case PMODEL_SCGAUSSIAN:
 			gr  = partiality_rgradient(rlow, r, pmodel);
 			gr -= partiality_rgradient(rhigh, r, pmodel);
 			return gr;
 
-			case PMODEL_THIN:
-			return 2.0*rlow*rlow/(r*r*r);
 		}
 	}
 
@@ -226,15 +215,14 @@ double p_gradient(Crystal *cr, int k, Reflection *refl, PartialityModel pmodel)
 
 			default:
 			case PMODEL_UNITY:
-			case PMODEL_THIN:
 			return 0.0;
 
-			case PMODEL_GAUSSIAN:
-			case PMODEL_SPHERE:
-			return (ds*glow + ds*ghigh) / 2.0;
-
 			case PMODEL_SCSPHERE:
 			return 0.0; /* FIXME */
+
+			case PMODEL_SCGAUSSIAN:
+			return (ds*glow + ds*ghigh) / 2.0; /* FIXME: Wrong */
+
 		}
 	}
 
@@ -277,28 +265,6 @@ double p_gradient(Crystal *cr, int k, Reflection *refl, PartialityModel pmodel)
 }
 
 
-/* Return the gradient of Lorentz factor wrt parameter 'k' given the current
- * status of 'image'. */
-double l_gradient(Crystal *cr, int k, Reflection *refl, PartialityModel pmodel)
-{
-	double ds;
-	signed int hs, ks, ls;
-	double L;
-
-	/* L has a non-zero gradient only for div in sphere or gaussian model */
-	if ( (pmodel != PMODEL_SPHERE)
-	  && (pmodel != PMODEL_GAUSSIAN) ) return 0.0;
-	if ( k != REF_DIV ) return 0.0;
-
-	get_symmetric_indices(refl, &hs, &ks, &ls);
-
-	ds = 2.0 * resolution(crystal_get_cell(cr), hs, ks, ls);
-
-	L = get_lorentz(refl);
-	return -ds*L*L / LORENTZ_SCALE;
-}
-
-
 static void apply_cell_shift(UnitCell *cell, int k, double shift)
 {
 	double asx, asy, asz;
@@ -564,10 +530,7 @@ static double pr_iterate(Crystal *cr, const RefList *full,
 
 		/* Calculate all gradients for this reflection */
 		for ( k=0; k<NUM_PARAMS; k++ ) {
-			double gr;
-			gr  = p_gradient(cr, k, refl, pmodel) * l;
-			gr += l_gradient(cr, k, refl, pmodel) * p;
-			gradients[k] = gr;
+			gradients[k] = p_gradient(cr, k, refl, pmodel) * l;
 		}
 
 		for ( k=0; k<NUM_PARAMS; k++ ) {
diff --git a/src/post-refinement.h b/src/post-refinement.h
index e419c51d..0c564690 100644
--- a/src/post-refinement.h
+++ b/src/post-refinement.h
@@ -3,11 +3,11 @@
  *
  * Post refinement
  *
- * Copyright © 2012 Deutsches Elektronen-Synchrotron DESY,
- *                  a research centre of the Helmholtz Association.
+ * Copyright © 2012-2014 Deutsches Elektronen-Synchrotron DESY,
+ *                       a research centre of the Helmholtz Association.
  *
  * Authors:
- *   2010-2012 Thomas White <taw@physics.org>
+ *   2010-2014 Thomas White <taw@physics.org>
  *
  * This file is part of CrystFEL.
  *
@@ -71,12 +71,8 @@ struct prdata
 extern struct prdata pr_refine(Crystal *cr, const RefList *full,
                                PartialityModel pmodel);
 
-/* Exported so it can be poked by tests/pr_gradient_check */
+/* Exported so it can be poked by tests/pr_p_gradient_check */
 extern double p_gradient(Crystal *cr, int k, Reflection *refl,
                          PartialityModel pmodel);
 
-
-extern double l_gradient(Crystal *cr, int k, Reflection *refl,
-                         PartialityModel pmodel);
-
 #endif	/* POST_REFINEMENT_H */
diff --git a/src/process_image.c b/src/process_image.c
index 83973e50..1d1dac82 100644
--- a/src/process_image.c
+++ b/src/process_image.c
@@ -178,7 +178,7 @@ void process_image(const struct index_args *iargs, struct pattern_args *pargs,
 
 	/* Integrate all the crystals at once - need all the crystals so that
 	 * overlaps can be detected. */
-	integrate_all_4(&image, iargs->int_meth, PMODEL_SPHERE, iargs->push_res,
+	integrate_all_4(&image, iargs->int_meth, PMODEL_SCSPHERE, iargs->push_res,
 	                iargs->ir_inn, iargs->ir_mid, iargs->ir_out,
 	                iargs->int_diag, iargs->int_diag_h,
 	                iargs->int_diag_k, iargs->int_diag_l, results_pipe);
diff --git a/tests/pr_l_gradient_check.c b/tests/pr_l_gradient_check.c
deleted file mode 100644
index b18196c9..00000000
--- a/tests/pr_l_gradient_check.c
+++ /dev/null
@@ -1,379 +0,0 @@
-/*
- * pr_l_gradient_check.c
- *
- * Check Lorentz factor gradients for post refinement
- *
- * Copyright © 2012-2014 Deutsches Elektronen-Synchrotron DESY,
- *                       a research centre of the Helmholtz Association.
- *
- * Authors:
- *   2012-2014 Thomas White <taw@physics.org>
- *
- * This file is part of CrystFEL.
- *
- * CrystFEL is free software: you can redistribute it and/or modify
- * it under the terms of the GNU General Public License as published by
- * the Free Software Foundation, either version 3 of the License, or
- * (at your option) any later version.
- *
- * CrystFEL is distributed in the hope that it will be useful,
- * but WITHOUT ANY WARRANTY; without even the implied warranty of
- * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
- * GNU General Public License for more details.
- *
- * You should have received a copy of the GNU General Public License
- * along with CrystFEL.  If not, see <http://www.gnu.org/licenses/>.
- *
- */
-
-#ifdef HAVE_CONFIG_H
-#include <config.h>
-#endif
-
-
-#include <stdlib.h>
-#include <stdio.h>
-#include <gsl/gsl_statistics.h>
-#include <getopt.h>
-
-#include <image.h>
-#include <cell.h>
-#include <cell-utils.h>
-#include <geometry.h>
-#include <reflist.h>
-#include "../src/post-refinement.h"
-
-
-static void scan_partialities(RefList *reflections, RefList *compare,
-                              int *valid, long double *vals[3], int idx)
-{
-	int i;
-	Reflection *refl;
-	RefListIterator *iter;
-
-	i = 0;
-	for ( refl = first_refl(reflections, &iter);
-	      refl != NULL;
-	      refl = next_refl(refl, iter) )
-	{
-		signed int h, k, l;
-		Reflection *refl2;
-
-		get_indices(refl, &h, &k, &l);
-		refl2 = find_refl(compare, h, k, l);
-		if ( refl2 == NULL ) {
-			valid[i] = 0;
-			i++;
-			continue;
-		}
-
-		vals[idx][i] = get_lorentz(refl2);
-		i++;
-	}
-}
-
-
-static void shift_parameter(struct image *image, int k, double shift)
-{
-	switch ( k )
-	{
-		case REF_DIV :
-		image->div += shift;
-		break;
-
-		default:
-		ERROR("This test can't check parameter %i\n", k);
-		exit(1);
-	}
-}
-
-
-static void calc_either_side(Crystal *cr, double incr_val,
-                             int *valid, long double *vals[3], int refine,
-                             PartialityModel pmodel)
-{
-	RefList *compare;
-	struct image *image = crystal_get_image(cr);
-	struct image im_moved;
-
-	im_moved = *image;
-	shift_parameter(&im_moved, refine, -incr_val);
-	compare = find_intersections(&im_moved, cr, pmodel);
-	scan_partialities(crystal_get_reflections(cr), compare,
-	                  valid, vals, 0);
-	reflist_free(compare);
-
-	im_moved = *image;
-	shift_parameter(&im_moved, refine, +incr_val);
-	compare = find_intersections(&im_moved, cr, pmodel);
-	scan_partialities(crystal_get_reflections(cr), compare,
-	                  valid, vals, 2);
-	reflist_free(compare);
-}
-
-
-
-static double test_gradients(Crystal *cr, double incr_val, int refine,
-                             const char *str, const char *file,
-                             PartialityModel pmodel, int quiet, int plot)
-{
-	Reflection *refl;
-	RefListIterator *iter;
-	long double *vals[3];
-	int i;
-	int *valid;
-	int nref;
-	int n_good, n_invalid, n_small, n_nan, n_bad;
-	RefList *reflections;
-	FILE *fh = NULL;
-	int ntot = 0;
-	double total = 0.0;
-	char tmp[32];
-	double *vec1;
-	double *vec2;
-	int n_line;
-	double cc;
-
-	reflections = find_intersections(crystal_get_image(cr), cr, pmodel);
-	crystal_set_reflections(cr, reflections);
-
-	nref = num_reflections(reflections);
-	if ( nref < 10 ) {
-		ERROR("Too few reflections found.  Failing test by default.\n");
-		return 0.0;
-	}
-
-	vals[0] = malloc(nref*sizeof(long double));
-	vals[1] = malloc(nref*sizeof(long double));
-	vals[2] = malloc(nref*sizeof(long double));
-	if ( (vals[0] == NULL) || (vals[1] == NULL) || (vals[2] == NULL) ) {
-		ERROR("Couldn't allocate memory.\n");
-		return 0.0;
-	}
-
-	valid = malloc(nref*sizeof(int));
-	if ( valid == NULL ) {
-		ERROR("Couldn't allocate memory.\n");
-		return 0.0;
-	}
-	for ( i=0; i<nref; i++ ) valid[i] = 1;
-
-	scan_partialities(reflections, reflections, valid, vals, 1);
-
-	calc_either_side(cr, incr_val, valid, vals, refine, pmodel);
-
-	if ( plot ) {
-		snprintf(tmp, 32, "gradient-test-%s.dat", file);
-		fh = fopen(tmp, "w");
-	}
-
-	vec1 = malloc(nref*sizeof(double));
-	vec2 = malloc(nref*sizeof(double));
-	if ( (vec1 == NULL) || (vec2 == NULL) ) {
-		ERROR("Couldn't allocate memory.\n");
-		return 0.0;
-	}
-
-	n_invalid = 0;  n_good = 0;
-	n_nan = 0;  n_small = 0;  n_bad = 0;  n_line = 0;
-	i = 0;
-	for ( refl = first_refl(reflections, &iter);
-	      refl != NULL;
-	      refl = next_refl(refl, iter) )
-	{
-
-		long double grad1, grad2, grad;
-		double cgrad;
-		signed int h, k, l;
-
-		get_indices(refl, &h, &k, &l);
-
-		if ( !valid[i] ) {
-			n_invalid++;
-			i++;
-		} else {
-
-			double r1, r2, p;
-			int cl, ch;
-
-			grad1 = (vals[1][i] - vals[0][i]) / incr_val;
-			grad2 = (vals[2][i] - vals[1][i]) / incr_val;
-			grad = (grad1 + grad2) / 2.0;
-			i++;
-
-			cgrad = l_gradient(cr, refine, refl, pmodel);
-
-			get_partial(refl, &r1, &r2, &p, &cl, &ch);
-
-			if ( isnan(cgrad) ) {
-				n_nan++;
-				continue;
-			}
-
-			if ( plot ) {
-				fprintf(fh, "%e %Le\n", cgrad, grad);
-			}
-
-			vec1[n_line] = cgrad;
-			vec2[n_line] = grad;
-			n_line++;
-
-			if ( (fabs(cgrad) < 5e-8) && (fabs(grad) < 5e-8) ) {
-				n_small++;
-				continue;
-			}
-
-			total += fabs(cgrad - grad);
-			ntot++;
-
-			if ( !within_tolerance(grad, cgrad, 5.0)
-			  || !within_tolerance(cgrad, grad, 5.0) )
-			{
-
-				if ( !quiet ) {
-					STATUS("!- %s %3i %3i %3i"
-					       " %10.2Le %10.2e ratio = %5.2Lf"
-					       " %10.2e %10.2e\n",
-					       str, h, k, l, grad, cgrad,
-					       cgrad/grad, r1, r2);
-				}
-				n_bad++;
-
-			} else {
-
-				//STATUS("OK %s %3i %3i %3i"
-				//       " %10.2Le %10.2e ratio = %5.2Lf"
-				//       " %10.2e %10.2e\n",
-				//       str, h, k, l, grad, cgrad, cgrad/grad,
-				//       r1, r2);
-
-				n_good++;
-
-			}
-
-		}
-
-	}
-
-	STATUS("%3s: %3i within 5%%, %3i outside, %3i nan, %3i invalid, "
-	       "%3i small. ", str, n_good, n_bad, n_nan, n_invalid, n_small);
-
-	if ( plot ) {
-		fclose(fh);
-	}
-
-	cc = gsl_stats_correlation(vec1, 1, vec2, 1, n_line);
-	STATUS("CC = %+f\n", cc);
-	return cc;
-}
-
-
-int main(int argc, char *argv[])
-{
-	struct image image;
-	const double incr_frac = 1.0/1000000.0;
-	double incr_val;
-	double ax, ay, az;
-	double bx, by, bz;
-	double cx, cy, cz;
-	UnitCell *cell;
-	Crystal *cr;
-	struct quaternion orientation;
-	int i;
-	int fail = 0;
-	int quiet = 0;
-	int plot = 0;
-	int c;
-	gsl_rng *rng;
-
-	const struct option longopts[] = {
-		{"quiet",       0, &quiet,        1},
-		{"plot",        0, &plot,         1},
-		{0, 0, NULL, 0}
-	};
-
-	while ((c = getopt_long(argc, argv, "", longopts, NULL)) != -1) {
-		switch (c) {
-
-			case 0 :
-			break;
-
-			case '?' :
-			break;
-
-			default :
-			ERROR("Unhandled option '%c'\n", c);
-			break;
-
-		}
-
-	}
-
-	image.width = 1024;
-	image.height = 1024;
-	image.det = simple_geometry(&image);
-	image.det->panels[0].res = 13333.3;
-	image.det->panels[0].clen = 80e-3;
-	image.det->panels[0].coffset = 0.0;
-
-	image.lambda = ph_en_to_lambda(eV_to_J(8000.0));
-	image.div = 1e-3;
-	image.bw = 0.01;
-	image.filename = malloc(256);
-
-	cr = crystal_new();
-	if ( cr == NULL ) {
-		ERROR("Failed to allocate crystal.\n");
-		return 1;
-	}
-	crystal_set_mosaicity(cr, 0.0);
-	crystal_set_profile_radius(cr, 0.005e9);
-	crystal_set_image(cr, &image);
-
-	cell = cell_new_from_parameters(10.0e-9, 10.0e-9, 10.0e-9,
-	                                deg2rad(90.0),
-	                                deg2rad(90.0),
-	                                deg2rad(90.0));
-
-	rng = gsl_rng_alloc(gsl_rng_mt19937);
-
-	for ( i=0; i<3; i++ ) {
-
-		UnitCell *rot;
-		double val;
-		PartialityModel pmodel;
-
-		if ( i == 0 ) {
-			pmodel = PMODEL_SPHERE;
-			STATUS("Testing flat sphere model:\n");
-		} else if ( i == 1 ) {
-			pmodel = PMODEL_GAUSSIAN;
-			STATUS("Testing Gaussian model:\n");
-		} else if ( i == 2 ) {
-			pmodel = PMODEL_THIN;
-			STATUS("No need to test thin Ewald sphere model.\n");
-			continue;
-		} else {
-			ERROR("WTF?\n");
-			return 1;
-		}
-
-		orientation = random_quaternion(rng);
-		rot = cell_rotate(cell, orientation);
-		crystal_set_cell(cr, rot);
-
-		cell_get_reciprocal(rot,
-			            &ax, &ay, &az, &bx, &by,
-			            &bz, &cx, &cy, &cz);
-
-		incr_val = incr_frac * image.div;
-		val =  test_gradients(cr, incr_val, REF_DIV, "div", "div",
-		                      pmodel, quiet, plot);
-		if ( val < 0.99 ) fail = 1;
-
-	}
-
-	gsl_rng_free(rng);
-
-	return fail;
-}
diff --git a/tests/pr_p_gradient_check.c b/tests/pr_p_gradient_check.c
index 82ed2b02..b224f9b4 100644
--- a/tests/pr_p_gradient_check.c
+++ b/tests/pr_p_gradient_check.c
@@ -71,21 +71,6 @@ static void scan_partialities(RefList *reflections, RefList *compare,
 		}
 
 		get_partial(refl2, &r1, &r2, &p, &clamp_low, &clamp_high);
-		if ( clamp_low && clamp_high && (pmodel != PMODEL_SCSPHERE) ) {
-			if ( !within_tolerance(p, 1.0, 0.001) ) {
-
-				signed int h, k, l;
-
-				get_indices(refl, &h, &k, &l);
-
-				ERROR("%3i %3i %3i - double clamped but"
-				      " partiality not close to 1.0 (%5.2f)\n",
-				      h, k, l, p);
-
-			}
-			valid[i] = 0;
-		}
-
 		vals[idx][i] = p;
 		i++;
 	}
@@ -451,26 +436,18 @@ int main(int argc, char *argv[])
 
 	rng = gsl_rng_alloc(gsl_rng_mt19937);
 
-	for ( i=0; i<4; i++ ) {
+	for ( i=0; i<2; i++ ) {
 
 		UnitCell *rot;
 		double val;
 		PartialityModel pmodel;
 
 		if ( i == 0 ) {
-			pmodel = PMODEL_SPHERE;
-			STATUS("Testing flat sphere model:\n");
-		} else if ( i == 1 ) {
-			pmodel = PMODEL_GAUSSIAN;
-			/* FIXME: Gradients for Gaussian model are not good */
-			STATUS("NOT testing Gaussian model.\n");
-			continue;
-		} else if ( i == 2 ) {
-			pmodel = PMODEL_THIN;
-			STATUS("Testing Thin Ewald Sphere model:\n");
-		} else if ( i == 3 ) {
 			pmodel = PMODEL_SCSPHERE;
 			STATUS("Testing SCSphere model:\n");
+		} else if ( i == 1 ) {
+			pmodel = PMODEL_SCGAUSSIAN;
+			STATUS("Testing SCGaussian model.\n");
 		} else {
 			ERROR("WTF?\n");
 			return 1;
diff --git a/tests/pr_pl_gradient_check.c b/tests/pr_pl_gradient_check.c
deleted file mode 100644
index cb68e0e3..00000000
--- a/tests/pr_pl_gradient_check.c
+++ /dev/null
@@ -1,540 +0,0 @@
-/*
- * pr_p_gradient_check.c
- *
- * Check partiality gradients for post refinement
- *
- * Copyright © 2012-2014 Deutsches Elektronen-Synchrotron DESY,
- *                       a research centre of the Helmholtz Association.
- *
- * Authors:
- *   2012-2014 Thomas White <taw@physics.org>
- *
- * This file is part of CrystFEL.
- *
- * CrystFEL is free software: you can redistribute it and/or modify
- * it under the terms of the GNU General Public License as published by
- * the Free Software Foundation, either version 3 of the License, or
- * (at your option) any later version.
- *
- * CrystFEL is distributed in the hope that it will be useful,
- * but WITHOUT ANY WARRANTY; without even the implied warranty of
- * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
- * GNU General Public License for more details.
- *
- * You should have received a copy of the GNU General Public License
- * along with CrystFEL.  If not, see <http://www.gnu.org/licenses/>.
- *
- */
-
-#ifdef HAVE_CONFIG_H
-#include <config.h>
-#endif
-
-
-#include <stdlib.h>
-#include <stdio.h>
-#include <gsl/gsl_statistics.h>
-#include <getopt.h>
-
-#include <image.h>
-#include <cell.h>
-#include <cell-utils.h>
-#include <geometry.h>
-#include <reflist.h>
-#include "../src/post-refinement.h"
-
-
-static void scan_partialities(RefList *reflections, RefList *compare,
-                              int *valid, long double *vals[3], int idx)
-{
-	int i;
-	Reflection *refl;
-	RefListIterator *iter;
-
-	i = 0;
-	for ( refl = first_refl(reflections, &iter);
-	      refl != NULL;
-	      refl = next_refl(refl, iter) )
-	{
-		signed int h, k, l;
-		Reflection *refl2;
-		double r1, r2, p;
-		int clamp_low, clamp_high;
-
-		get_indices(refl, &h, &k, &l);
-		refl2 = find_refl(compare, h, k, l);
-		if ( refl2 == NULL ) {
-			valid[i] = 0;
-			i++;
-			continue;
-		}
-
-		get_partial(refl2, &r1, &r2, &p, &clamp_low, &clamp_high);
-		if ( clamp_low && clamp_high ) {
-			if ( !within_tolerance(p, 1.0, 0.001) ) {
-
-				signed int h, k, l;
-
-				get_indices(refl, &h, &k, &l);
-
-				ERROR("%3i %3i %3i - double clamped but"
-				      " partiality not close to 1.0 (%5.2f)\n",
-				      h, k, l, p);
-
-			}
-			valid[i] = 0;
-		}
-
-		vals[idx][i] = p * get_lorentz(refl2);
-		i++;
-	}
-}
-
-
-static UnitCell *new_shifted_cell(UnitCell *input, int k, double shift)
-{
-	UnitCell *cell;
-	double asx, asy, asz;
-	double bsx, bsy, bsz;
-	double csx, csy, csz;
-
-	cell = cell_new();
-	cell_get_reciprocal(input, &asx, &asy, &asz, &bsx, &bsy, &bsz,
-	                    &csx, &csy, &csz);
-	switch ( k )
-	{
-		case REF_ASX :  asx += shift;  break;
-		case REF_ASY :  asy += shift;  break;
-		case REF_ASZ :  asz += shift;  break;
-		case REF_BSX :  bsx += shift;  break;
-		case REF_BSY :  bsy += shift;  break;
-		case REF_BSZ :  bsz += shift;  break;
-		case REF_CSX :  csx += shift;  break;
-		case REF_CSY :  csy += shift;  break;
-		case REF_CSZ :  csz += shift;  break;
-	}
-	cell_set_reciprocal(cell, asx, asy, asz, bsx, bsy, bsz, csx, csy, csz);
-
-	return cell;
-}
-
-
-static void shift_parameter(struct image *image, int k, double shift)
-{
-	switch ( k )
-	{
-		case REF_DIV : image->div += shift;  break;
-	}
-}
-
-
-static Crystal *new_shifted_crystal(Crystal *cr, int refine, double incr_val)
-{
-	Crystal *cr_new;
-	double r;
-	UnitCell *cell;
-
-	cr_new = crystal_copy(cr);
-	if ( cr_new == NULL ) {
-		ERROR("Failed to allocate crystal.\n");
-		return NULL;
-	}
-
-	crystal_set_image(cr_new, crystal_get_image(cr));
-	r = crystal_get_profile_radius(cr_new);
-
-	switch ( refine ) {
-
-		case REF_ASX :
-		case REF_ASY :
-		case REF_ASZ :
-		case REF_BSX :
-		case REF_BSY :
-		case REF_BSZ :
-		case REF_CSX :
-		case REF_CSY :
-		case REF_CSZ :
-		cell = new_shifted_cell(crystal_get_cell(cr), refine,
-		                        incr_val);
-		crystal_set_cell(cr_new, cell);
-		break;
-
-		case REF_R :
-		cell = cell_new_from_cell(crystal_get_cell(cr));
-		crystal_set_cell(cr_new, cell);
-		crystal_set_profile_radius(cr_new, r + incr_val);
-		break;
-
-		default :
-		ERROR("Can't shift %i\n", refine);
-		break;
-
-	}
-
-	return cr_new;
-}
-
-
-static void calc_either_side(Crystal *cr, double incr_val,
-                             int *valid, long double *vals[3], int refine,
-                             PartialityModel pmodel)
-{
-	RefList *compare;
-	struct image *image = crystal_get_image(cr);
-
-	if ( (refine != REF_DIV) ) {
-
-		Crystal *cr_new;
-
-		/* Crystal properties */
-		cr_new = new_shifted_crystal(cr, refine, -incr_val);
-		compare = find_intersections(image, cr_new, pmodel);
-		scan_partialities(crystal_get_reflections(cr), compare, valid,
-		                  vals, 0);
-		cell_free(crystal_get_cell(cr_new));
-		crystal_free(cr_new);
-		reflist_free(compare);
-
-		cr_new = new_shifted_crystal(cr, refine, +incr_val);
-		compare = find_intersections(image, cr_new, pmodel);
-		scan_partialities(crystal_get_reflections(cr), compare, valid,
-		                  vals, 2);
-		cell_free(crystal_get_cell(cr_new));
-		crystal_free(cr_new);
-		reflist_free(compare);
-
-	} else {
-
-		struct image im_moved;
-
-		/* "Image" properties */
-		im_moved = *image;
-		shift_parameter(&im_moved, refine, -incr_val);
-		compare = find_intersections(&im_moved, cr, pmodel);
-		scan_partialities(crystal_get_reflections(cr), compare,
-		                  valid, vals, 0);
-		reflist_free(compare);
-
-		im_moved = *image;
-		shift_parameter(&im_moved, refine, +incr_val);
-		compare = find_intersections(&im_moved, cr, pmodel);
-		scan_partialities(crystal_get_reflections(cr), compare,
-		                  valid, vals, 2);
-		reflist_free(compare);
-
-	}
-}
-
-
-static double test_gradients(Crystal *cr, double incr_val, int refine,
-                             const char *str, const char *file,
-                             PartialityModel pmodel, int quiet, int plot)
-{
-	Reflection *refl;
-	RefListIterator *iter;
-	long double *vals[3];
-	int i;
-	int *valid;
-	int nref;
-	int n_good, n_invalid, n_small, n_nan, n_bad;
-	RefList *reflections;
-	FILE *fh = NULL;
-	int ntot = 0;
-	double total = 0.0;
-	char tmp[32];
-	double *vec1;
-	double *vec2;
-	int n_line;
-	double cc;
-
-	reflections = find_intersections(crystal_get_image(cr), cr, pmodel);
-	crystal_set_reflections(cr, reflections);
-
-	nref = num_reflections(reflections);
-	if ( nref < 10 ) {
-		ERROR("Too few reflections found.  Failing test by default.\n");
-		return 0.0;
-	}
-
-	vals[0] = malloc(nref*sizeof(long double));
-	vals[1] = malloc(nref*sizeof(long double));
-	vals[2] = malloc(nref*sizeof(long double));
-	if ( (vals[0] == NULL) || (vals[1] == NULL) || (vals[2] == NULL) ) {
-		ERROR("Couldn't allocate memory.\n");
-		return 0.0;
-	}
-
-	valid = malloc(nref*sizeof(int));
-	if ( valid == NULL ) {
-		ERROR("Couldn't allocate memory.\n");
-		return 0.0;
-	}
-	for ( i=0; i<nref; i++ ) valid[i] = 1;
-
-	scan_partialities(reflections, reflections, valid, vals, 1);
-
-	calc_either_side(cr, incr_val, valid, vals, refine, pmodel);
-
-	if ( plot ) {
-		snprintf(tmp, 32, "gradient-test-%s.dat", file);
-		fh = fopen(tmp, "w");
-	}
-
-	vec1 = malloc(nref*sizeof(double));
-	vec2 = malloc(nref*sizeof(double));
-	if ( (vec1 == NULL) || (vec2 == NULL) ) {
-		ERROR("Couldn't allocate memory.\n");
-		return 0.0;
-	}
-
-	n_invalid = 0;  n_good = 0;
-	n_nan = 0;  n_small = 0;  n_bad = 0;  n_line = 0;
-	i = 0;
-	for ( refl = first_refl(reflections, &iter);
-	      refl != NULL;
-	      refl = next_refl(refl, iter) )
-	{
-
-		long double grad1, grad2, grad;
-		double cgrad;
-		signed int h, k, l;
-
-		get_indices(refl, &h, &k, &l);
-
-		if ( !valid[i] ) {
-			n_invalid++;
-			i++;
-		} else {
-
-			double r1, r2, p, lor;
-			int cl, ch;
-
-			grad1 = (vals[1][i] - vals[0][i]) / incr_val;
-			grad2 = (vals[2][i] - vals[1][i]) / incr_val;
-			grad = (grad1 + grad2) / 2.0;
-			i++;
-
-			get_partial(refl, &r1, &r2, &p, &cl, &ch);
-			lor = get_lorentz(refl);
-
-			cgrad  = p_gradient(cr, refine, refl, pmodel) * lor;
-			cgrad += l_gradient(cr, refine, refl, pmodel) * p;
-
-			if ( isnan(cgrad) ) {
-				n_nan++;
-				continue;
-			}
-
-			if ( plot ) {
-				fprintf(fh, "%e %Le\n", cgrad, grad);
-			}
-
-			vec1[n_line] = cgrad;
-			vec2[n_line] = grad;
-			n_line++;
-
-			if ( (fabs(cgrad) < 5e-8) && (fabs(grad) < 5e-8) ) {
-				n_small++;
-				continue;
-			}
-
-			total += fabs(cgrad - grad);
-			ntot++;
-
-			if ( !within_tolerance(grad, cgrad, 5.0)
-			  || !within_tolerance(cgrad, grad, 5.0) )
-			{
-
-				if ( !quiet ) {
-					STATUS("!- %s %3i %3i %3i"
-					       " %10.2Le %10.2e ratio = %5.2Lf"
-					       " %10.2e %10.2e\n",
-					       str, h, k, l, grad, cgrad,
-					       cgrad/grad, r1, r2);
-				}
-				n_bad++;
-
-			} else {
-
-				//STATUS("OK %s %3i %3i %3i"
-				//       " %10.2Le %10.2e ratio = %5.2Lf"
-				//       " %10.2e %10.2e\n",
-				//       str, h, k, l, grad, cgrad, cgrad/grad,
-				//       r1, r2);
-
-				n_good++;
-
-			}
-
-		}
-
-	}
-
-	STATUS("%3s: %3i within 5%%, %3i outside, %3i nan, %3i invalid, "
-	       "%3i small. ", str, n_good, n_bad, n_nan, n_invalid, n_small);
-
-	if ( plot ) {
-		fclose(fh);
-	}
-
-	cc = gsl_stats_correlation(vec1, 1, vec2, 1, n_line);
-	STATUS("CC = %+f\n", cc);
-	return cc;
-}
-
-
-int main(int argc, char *argv[])
-{
-	struct image image;
-	const double incr_frac = 1.0/1000000.0;
-	double incr_val;
-	double ax, ay, az;
-	double bx, by, bz;
-	double cx, cy, cz;
-	UnitCell *cell;
-	Crystal *cr;
-	struct quaternion orientation;
-	int i;
-	int fail = 0;
-	int quiet = 0;
-	int plot = 0;
-	int c;
-	gsl_rng *rng;
-
-	const struct option longopts[] = {
-		{"quiet",       0, &quiet,        1},
-		{"plot",        0, &plot,         1},
-		{0, 0, NULL, 0}
-	};
-
-	while ((c = getopt_long(argc, argv, "", longopts, NULL)) != -1) {
-		switch (c) {
-
-			case 0 :
-			break;
-
-			case '?' :
-			break;
-
-			default :
-			ERROR("Unhandled option '%c'\n", c);
-			break;
-
-		}
-
-	}
-
-	image.width = 1024;
-	image.height = 1024;
-	image.det = simple_geometry(&image);
-	image.det->panels[0].res = 13333.3;
-	image.det->panels[0].clen = 80e-3;
-	image.det->panels[0].coffset = 0.0;
-
-	image.lambda = ph_en_to_lambda(eV_to_J(8000.0));
-	image.div = 1e-3;
-	image.bw = 0.01;
-	image.filename = malloc(256);
-
-	cr = crystal_new();
-	if ( cr == NULL ) {
-		ERROR("Failed to allocate crystal.\n");
-		return 1;
-	}
-	crystal_set_mosaicity(cr, 0.0);
-	crystal_set_profile_radius(cr, 0.005e9);
-	crystal_set_image(cr, &image);
-
-	cell = cell_new_from_parameters(10.0e-9, 10.0e-9, 10.0e-9,
-	                                deg2rad(90.0),
-	                                deg2rad(90.0),
-	                                deg2rad(90.0));
-
-	rng = gsl_rng_alloc(gsl_rng_mt19937);
-
-	for ( i=0; i<3; i++ ) {
-
-		UnitCell *rot;
-		double val;
-		PartialityModel pmodel;
-
-		if ( i == 0 ) {
-			pmodel = PMODEL_SPHERE;
-			STATUS("Testing flat sphere model:\n");
-		} else if ( i == 1 ) {
-			pmodel = PMODEL_GAUSSIAN;
-			/* FIXME: Gradients for Gaussian model are not good */
-			STATUS("NOT testing Gaussian model.\n");
-			continue;
-		} else if ( i == 2 ) {
-			pmodel = PMODEL_THIN;
-			STATUS("No need to test thin Ewald sphere model.\n");
-			continue;
-		} else {
-			ERROR("WTF?\n");
-			return 1;
-		}
-
-		orientation = random_quaternion(rng);
-		rot = cell_rotate(cell, orientation);
-		crystal_set_cell(cr, rot);
-
-		cell_get_reciprocal(rot,
-			            &ax, &ay, &az, &bx, &by,
-			            &bz, &cx, &cy, &cz);
-
-		incr_val = incr_frac * image.div;
-		val =  test_gradients(cr, incr_val, REF_DIV, "div", "div",
-		                      pmodel, quiet, plot);
-		if ( val < 0.99 ) fail = 1;
-
-		incr_val = incr_frac * crystal_get_profile_radius(cr);
-		val = test_gradients(cr, incr_val, REF_R, "R", "R", pmodel,
-		                     quiet, plot);
-		if ( val < 0.99 ) fail = 1;
-
-		incr_val = incr_frac * ax;
-		val = test_gradients(cr, incr_val, REF_ASX, "ax*", "x", pmodel,
-		                     quiet, plot);
-		if ( val < 0.99 ) fail = 1;
-		incr_val = incr_frac * bx;
-		val = test_gradients(cr, incr_val, REF_BSX, "bx*", "x", pmodel,
-		                     quiet, plot);
-		if ( val < 0.99 ) fail = 1;
-		incr_val = incr_frac * cx;
-		val = test_gradients(cr, incr_val, REF_CSX, "cx*", "x", pmodel,
-		                     quiet, plot);
-		if ( val < 0.99 ) fail = 1;
-
-		incr_val = incr_frac * ay;
-		val = test_gradients(cr, incr_val, REF_ASY, "ay*", "y", pmodel,
-		                     quiet, plot);
-		if ( val < 0.99 ) fail = 1;
-		incr_val = incr_frac * by;
-		val = test_gradients(cr, incr_val, REF_BSY, "by*", "y", pmodel,
-		                     quiet, plot);
-		if ( val < 0.99 ) fail = 1;
-		incr_val = incr_frac * cy;
-		val = test_gradients(cr, incr_val, REF_CSY, "cy*", "y", pmodel,
-		                     quiet, plot);
-		if ( val < 0.99 ) fail = 1;
-
-		incr_val = incr_frac * az;
-		val = test_gradients(cr, incr_val, REF_ASZ, "az*", "z", pmodel,
-		                     quiet, plot);
-		if ( val < 0.99 ) fail = 1;
-		incr_val = incr_frac * bz;
-		val = test_gradients(cr, incr_val, REF_BSZ, "bz*", "z", pmodel,
-		                     quiet, plot);
-		if ( val < 0.99 ) fail = 1;
-		incr_val = incr_frac * cz;
-		val = test_gradients(cr, incr_val, REF_CSZ, "cz*", "z", pmodel,
-		                     quiet, plot);
-		if ( val < 0.99 ) fail = 1;
-
-	}
-
-	gsl_rng_free(rng);
-
-	return fail;
-}
diff --git a/tests/prof2d_check.c b/tests/prof2d_check.c
index 8563c7d9..c7ebd8b2 100644
--- a/tests/prof2d_check.c
+++ b/tests/prof2d_check.c
@@ -138,7 +138,7 @@ int main(int argc, char *argv[])
 	image.n_crystals = 1;
 	image.crystals = &cr;
 
-	list = find_intersections(&image, cr, PMODEL_SPHERE);
+	list = find_intersections(&image, cr, PMODEL_SCSPHERE);
 	crystal_set_reflections(cr, list);
 
 	for ( fs=0; fs<w; fs++ ) {