Move a load more stuff into libcrystfel

20 files changed, 0 insertions, 5833 deletions

diff --git a/src/cl-utils.c b/src/cl-utils.c
deleted file mode 100644
index 65a09363..00000000
--- a/src/cl-utils.c
+++ /dev/null
@@ -1,200 +0,0 @@
-/*
- * cl-utils.c
- *
- * OpenCL utility functions
- *
- * (c) 2006-2011 Thomas White <taw@physics.org>
- *
- * Part of CrystFEL - crystallography with a FEL
- *
- */
-
-#ifdef HAVE_CONFIG_H
-#include <config.h>
-#endif
-
-#include <stdlib.h>
-#include <stdio.h>
-#include <string.h>
-
-#ifdef HAVE_CL_CL_H
-#include <CL/cl.h>
-#else
-#include <cl.h>
-#endif
-
-#include "utils.h"
-
-
-const char *clError(cl_int err)
-{
-	switch ( err ) {
-	case CL_SUCCESS : return "no error";
-	case CL_INVALID_PLATFORM : return "invalid platform";
-	case CL_INVALID_KERNEL : return "invalid kernel";
-	case CL_INVALID_ARG_INDEX : return "invalid argument index";
-	case CL_INVALID_ARG_VALUE : return "invalid argument value";
-	case CL_INVALID_MEM_OBJECT : return "invalid memory object";
-	case CL_INVALID_SAMPLER : return "invalid sampler";
-	case CL_INVALID_ARG_SIZE : return "invalid argument size";
-	case CL_INVALID_COMMAND_QUEUE  : return "invalid command queue";
-	case CL_INVALID_CONTEXT : return "invalid context";
-	case CL_INVALID_VALUE : return "invalid value";
-	case CL_INVALID_EVENT_WAIT_LIST : return "invalid wait list";
-	case CL_MAP_FAILURE : return "map failure";
-	case CL_MEM_OBJECT_ALLOCATION_FAILURE : return "object allocation failure";
-	case CL_OUT_OF_HOST_MEMORY : return "out of host memory";
-	case CL_OUT_OF_RESOURCES : return "out of resources";
-	case CL_INVALID_KERNEL_NAME : return "invalid kernel name";
-	case CL_INVALID_KERNEL_ARGS : return "invalid kernel arguments";
-	case CL_INVALID_WORK_GROUP_SIZE : return "invalid work group size";
-	case CL_IMAGE_FORMAT_NOT_SUPPORTED : return "image format not supported";
-	case CL_INVALID_WORK_DIMENSION : return "invalid work dimension";
-	default :
-		return "unknown error";
-	}
-}
-
-
-static char *get_device_string(cl_device_id dev, cl_device_info info)
-{
-	int r;
-	size_t size;
-	char *val;
-
-	r = clGetDeviceInfo(dev, info, 0, NULL, &size);
-	if ( r != CL_SUCCESS ) {
-		ERROR("Couldn't get device vendor size: %s\n",
-		      clError(r));
-		return NULL;
-	}
-	val = malloc(size);
-	r = clGetDeviceInfo(dev, info, size, val, NULL);
-	if ( r != CL_SUCCESS ) {
-		ERROR("Couldn't get dev vendor: %s\n", clError(r));
-		return NULL;
-	}
-
-	return val;
-}
-
-
-cl_device_id get_cl_dev(cl_context ctx, int n)
-{
-	cl_device_id *dev;
-	cl_int r;
-	size_t size;
-	int i, num_devs;
-
-	/* Get the required size of the array */
-	r = clGetContextInfo(ctx, CL_CONTEXT_DEVICES, 0, NULL, &size);
-	if ( r != CL_SUCCESS ) {
-		ERROR("Couldn't get array size for devices: %s\n", clError(r));
-		return 0;
-	}
-
-	dev = malloc(size);
-	r = clGetContextInfo(ctx, CL_CONTEXT_DEVICES, size, dev, NULL);
-	if ( r != CL_SUCCESS ) {
-		ERROR("Couldn't get device: %s\n", clError(r));
-		return 0;
-	}
-	num_devs = size/sizeof(cl_device_id);
-
-	if ( n >= num_devs ) {
-		ERROR("Device ID out of range\n");
-		return 0;
-	}
-
-	if ( n < 0 ) {
-
-		STATUS("Available devices:\n");
-		for ( i=0; i<num_devs; i++ ) {
-
-			char *vendor;
-			char *name;
-
-			vendor = get_device_string(dev[i], CL_DEVICE_VENDOR);
-			name = get_device_string(dev[i], CL_DEVICE_NAME);
-
-			STATUS("Device %i: %s %s\n", i, vendor, name);
-
-		}
-		n = 0;
-
-		STATUS("Using device 0.  Use --gpu-dev to choose another.\n");
-
-	} else {
-
-		char *vendor;
-		char *name;
-
-		vendor = get_device_string(dev[n], CL_DEVICE_VENDOR);
-		name = get_device_string(dev[n], CL_DEVICE_NAME);
-
-		STATUS("Using device %i: %s %s\n", n, vendor, name);
-
-	}
-
-	return dev[n];
-}
-
-
-static void show_build_log(cl_program prog, cl_device_id dev)
-{
-	cl_int r;
-	char log[4096];
-	size_t s;
-
-	r = clGetProgramBuildInfo(prog, dev, CL_PROGRAM_BUILD_LOG, 4096, log,
-	                          &s);
-
-	STATUS("%s\n", log);
-}
-
-
-cl_program load_program(const char *filename, cl_context ctx,
-                        cl_device_id dev, cl_int *err, const char *extra_cflags)
-{
-	FILE *fh;
-	cl_program prog;
-	char *source;
-	size_t len;
-	cl_int r;
-	char cflags[1024] = "";
-
-	fh = fopen(filename, "r");
-	if ( fh == NULL ) {
-		ERROR("Couldn't open '%s'\n", filename);
-		*err = CL_INVALID_PROGRAM;
-		return 0;
-	}
-	source = malloc(16384);
-	len = fread(source, 1, 16383, fh);
-	fclose(fh);
-	source[len] = '\0';
-
-	prog = clCreateProgramWithSource(ctx, 1, (const char **)&source,
-	                                 NULL, err);
-	if ( *err != CL_SUCCESS ) {
-		ERROR("Couldn't load source\n");
-		return 0;
-	}
-
-	snprintf(cflags, 1023, "-Werror ");
-	strncat(cflags, "-I"DATADIR"/crystfel/ ", 1023-strlen(cflags));
-	strncat(cflags, "-cl-no-signed-zeros ", 1023-strlen(cflags));
-	strncat(cflags, extra_cflags, 1023-strlen(cflags));
-
-	r = clBuildProgram(prog, 0, NULL, cflags, NULL, NULL);
-	if ( r != CL_SUCCESS ) {
-		ERROR("Couldn't build program '%s'\n", filename);
-		show_build_log(prog, dev);
-		*err = r;
-		return 0;
-	}
-
-	free(source);
-	*err = CL_SUCCESS;
-	return prog;
-}
diff --git a/src/cl-utils.h b/src/cl-utils.h
deleted file mode 100644
index 21a7ecd2..00000000
--- a/src/cl-utils.h
+++ /dev/null
@@ -1,27 +0,0 @@
-/*
- * cl-utils.h
- *
- * OpenCL utility functions
- *
- * (c) 2006-2010 Thomas White <taw@physics.org>
- *
- * Part of CrystFEL - crystallography with a FEL
- *
- */
-
-#ifndef CLUTILS_H
-#define CLUTILS_H
-
-#ifdef HAVE_CONFIG_H
-#include <config.h>
-#endif
-
-
-extern const char *clError(cl_int err);
-extern cl_device_id get_cl_dev(cl_context ctx, int n);
-extern cl_program load_program(const char *filename, cl_context ctx,
-                        cl_device_id dev, cl_int *err,
-                        const char *extra_cflags);
-
-
-#endif	/* CLUTILS_H */
diff --git a/src/diffraction-gpu.c b/src/diffraction-gpu.c
deleted file mode 100644
index 605b1514..00000000
--- a/src/diffraction-gpu.c
+++ /dev/null
@@ -1,529 +0,0 @@
-/*
- * diffraction-gpu.c
- *
- * Calculate diffraction patterns by Fourier methods (GPU version)
- *
- * (c) 2006-2011 Thomas White <taw@physics.org>
- *
- * Part of CrystFEL - crystallography with a FEL
- *
- */
-
-#ifdef HAVE_CONFIG_H
-#include <config.h>
-#endif
-
-#include <stdlib.h>
-#include <math.h>
-#include <stdio.h>
-#include <string.h>
-#include <complex.h>
-
-#ifdef HAVE_CL_CL_H
-#include <CL/cl.h>
-#else
-#include <cl.h>
-#endif
-
-#include "image.h"
-#include "utils.h"
-#include "cell.h"
-#include "diffraction.h"
-#include "cl-utils.h"
-#include "beam-parameters.h"
-
-
-#define SAMPLING (4)
-#define BWSAMPLING (10)
-#define DIVSAMPLING (1)
-#define SINC_LUT_ELEMENTS (4096)
-
-
-struct gpu_context
-{
-	cl_context ctx;
-	cl_command_queue cq;
-	cl_program prog;
-	cl_kernel kern;
-	cl_mem intensities;
-	cl_mem flags;
-
-	/* Array of sinc LUTs */
-	cl_mem *sinc_luts;
-	cl_float **sinc_lut_ptrs;
-	int max_sinc_lut;  /* Number of LUTs, i.e. one greater than the maximum
-	                    * index.  This equals the highest allowable "n". */
-};
-
-
-static void check_sinc_lut(struct gpu_context *gctx, int n)
-{
-	cl_int err;
-	cl_image_format fmt;
-	int i;
-
-	if ( n > gctx->max_sinc_lut ) {
-
-		gctx->sinc_luts = realloc(gctx->sinc_luts,
-		                          n*sizeof(*gctx->sinc_luts));
-		gctx->sinc_lut_ptrs = realloc(gctx->sinc_lut_ptrs,
-		                              n*sizeof(*gctx->sinc_lut_ptrs));
-
-		for ( i=gctx->max_sinc_lut; i<n; i++ ) {
-			gctx->sinc_lut_ptrs[i] = NULL;
-		}
-
-		gctx->max_sinc_lut = n;
-	}
-
-	if ( gctx->sinc_lut_ptrs[n-1] != NULL ) return;
-
-	/* Create a new sinc LUT */
-	gctx->sinc_lut_ptrs[n-1] = malloc(SINC_LUT_ELEMENTS*sizeof(cl_float));
-	gctx->sinc_lut_ptrs[n-1][0] = n;
-	if ( n == 1 ) {
-		for ( i=1; i<SINC_LUT_ELEMENTS; i++ ) {
-			gctx->sinc_lut_ptrs[n-1][i] = 1.0;
-		}
-	} else {
-		for ( i=1; i<SINC_LUT_ELEMENTS; i++ ) {
-			double x, val;
-			x = (double)i/SINC_LUT_ELEMENTS;
-			val = fabs(sin(M_PI*n*x)/sin(M_PI*x));
-			gctx->sinc_lut_ptrs[n-1][i] = val;
-		}
-	}
-
-	fmt.image_channel_order = CL_INTENSITY;
-	fmt.image_channel_data_type = CL_FLOAT;
-
-	gctx->sinc_luts[n-1] = clCreateImage2D(gctx->ctx,
-	                                CL_MEM_READ_ONLY | CL_MEM_COPY_HOST_PTR,
-	                                &fmt, SINC_LUT_ELEMENTS, 1, 0,
-	                                gctx->sinc_lut_ptrs[n-1], &err);
-	if ( err != CL_SUCCESS ) {
-		ERROR("Couldn't create LUT for %i\n", n);
-		return;
-	}
-}
-
-
-static int set_arg_float(struct gpu_context *gctx, int idx, float val)
-{
-	cl_int err;
-	err = clSetKernelArg(gctx->kern, idx, sizeof(cl_float), &val);
-	if ( err != CL_SUCCESS ) {
-		ERROR("Couldn't set kernel argument %i: %s\n",
-		      idx, clError(err));
-		return 1;
-	}
-
-	return 0;
-}
-
-
-static int set_arg_int(struct gpu_context *gctx, int idx, int val)
-{
-	cl_int err;
-
-	err = clSetKernelArg(gctx->kern, idx, sizeof(cl_int), &val);
-	if ( err != CL_SUCCESS ) {
-		ERROR("Couldn't set kernel argument %i: %s\n",
-		      idx, clError(err));
-		return 1;
-	}
-
-	return 0;
-}
-
-
-static int set_arg_mem(struct gpu_context *gctx, int idx, cl_mem val)
-{
-	cl_int err;
-
-	err = clSetKernelArg(gctx->kern, idx, sizeof(cl_mem), &val);
-	if ( err != CL_SUCCESS ) {
-		ERROR("Couldn't set kernel argument %i: %s\n",
-		      idx, clError(err));
-		return 1;
-	}
-
-	return 0;
-}
-
-
-void get_diffraction_gpu(struct gpu_context *gctx, struct image *image,
-                         int na, int nb, int nc, UnitCell *ucell)
-{
-	cl_int err;
-	double ax, ay, az;
-	double bx, by, bz;
-	double cx, cy, cz;
-	float klow, khigh;
-	int i;
-	cl_float16 cell;
-	cl_int4 ncells;
-	const int sampling = SAMPLING;
-	cl_float bwstep;
-	int n_inf = 0;
-	int n_neg = 0;
-	cl_float divxlow, divxstep;
-	cl_float divylow, divystep;
-	int n_nan = 0;
-	int sprod;
-
-	if ( gctx == NULL ) {
-		ERROR("GPU setup failed.\n");
-		return;
-	}
-
-	cell_get_cartesian(ucell, &ax, &ay, &az, &bx, &by, &bz, &cx, &cy, &cz);
-	cell.s[0] = ax;  cell.s[1] = ay;  cell.s[2] = az;
-	cell.s[3] = bx;  cell.s[4] = by;  cell.s[5] = bz;
-	cell.s[6] = cx;  cell.s[7] = cy;  cell.s[8] = cz;
-
-	/* Calculate wavelength */
-	klow = 1.0/(image->lambda*(1.0 + image->beam->bandwidth/2.0));
-	khigh = 1.0/(image->lambda*(1.0 - image->beam->bandwidth/2.0));
-	bwstep = (khigh-klow) / BWSAMPLING;
-
-	/* Calculate divergence stuff */
-	divxlow = -image->beam->divergence/2.0;
-	divylow = -image->beam->divergence/2.0;
-	divxstep = image->beam->divergence / DIVSAMPLING;
-	divystep = image->beam->divergence / DIVSAMPLING;
-
-	ncells.s[0] = na;
-	ncells.s[1] = nb;
-	ncells.s[2] = nc;
-	ncells.s[3] = 0;  /* unused */
-
-	/* Ensure all required LUTs are available */
-	check_sinc_lut(gctx, na);
-	check_sinc_lut(gctx, nb);
-	check_sinc_lut(gctx, nc);
-
-	if ( set_arg_float(gctx, 2, klow) ) return;
-	if ( set_arg_mem(gctx, 9, gctx->intensities) ) return;
-	if ( set_arg_int(gctx, 12, sampling) ) return;
-	if ( set_arg_float(gctx, 14, bwstep) ) return;
-	if ( set_arg_mem(gctx, 15, gctx->sinc_luts[na-1]) ) return;
-	if ( set_arg_mem(gctx, 16, gctx->sinc_luts[nb-1]) ) return;
-	if ( set_arg_mem(gctx, 17, gctx->sinc_luts[nc-1]) ) return;
-	if ( set_arg_mem(gctx, 18, gctx->flags) ) return;
-	if ( set_arg_float(gctx, 23, divxlow) ) return;
-	if ( set_arg_float(gctx, 24, divxstep) ) return;
-	if ( set_arg_int(gctx, 25, DIVSAMPLING) ) return;
-	if ( set_arg_float(gctx, 26, divylow) ) return;
-	if ( set_arg_float(gctx, 27, divystep) ) return;
-	if ( set_arg_int(gctx, 28, DIVSAMPLING) ) return;
-
-	/* Unit cell */
-	err = clSetKernelArg(gctx->kern, 8, sizeof(cl_float16), &cell);
-	if ( err != CL_SUCCESS ) {
-		ERROR("Couldn't set unit cell: %s\n", clError(err));
-		return;
-	}
-
-	/* Local memory for reduction */
-	sprod = BWSAMPLING*SAMPLING*SAMPLING*DIVSAMPLING*DIVSAMPLING;
-	err = clSetKernelArg(gctx->kern, 13, sprod*sizeof(cl_float), NULL);
-	if ( err != CL_SUCCESS ) {
-		ERROR("Couldn't set local memory: %s\n", clError(err));
-		return;
-	}
-
-	/* Allocate memory for the result */
-	image->data = calloc(image->width * image->height, sizeof(float));
-	image->twotheta = calloc(image->width * image->height, sizeof(double));
-
-	/* Iterate over panels */
-	for ( i=0; i<image->det->n_panels; i++ ) {
-
-		size_t dims[3];
-		size_t ldims[3] = {SAMPLING, SAMPLING,
-		                   BWSAMPLING * DIVSAMPLING * DIVSAMPLING};
-		struct panel *p;
-		cl_mem tt;
-		size_t tt_size;
-		cl_mem diff;
-		size_t diff_size;
-		float *diff_ptr;
-		float *tt_ptr;
-		int pan_width, pan_height;
-		int fs, ss;
-
-		p = &image->det->panels[i];
-
-		pan_width = 1 + p->max_fs - p->min_fs;
-		pan_height = 1 + p->max_ss - p->min_ss;
-
-		/* Buffer for the results of this panel */
-		diff_size = pan_width * pan_height * sizeof(cl_float);
-		diff = clCreateBuffer(gctx->ctx, CL_MEM_WRITE_ONLY,
-	                              diff_size, NULL, &err);
-		if ( err != CL_SUCCESS ) {
-			ERROR("Couldn't allocate diffraction memory\n");
-			return;
-		}
-		tt_size = pan_width * pan_height * sizeof(cl_float);
-		tt = clCreateBuffer(gctx->ctx, CL_MEM_WRITE_ONLY, tt_size,
-			            NULL, &err);
-		if ( err != CL_SUCCESS ) {
-			ERROR("Couldn't allocate twotheta memory\n");
-			return;
-		}
-
-		if ( set_arg_mem(gctx, 0, diff) ) return;
-		if ( set_arg_mem(gctx, 1, tt) ) return;
-		if ( set_arg_int(gctx, 3, pan_width) ) return;
-		if ( set_arg_float(gctx, 4, p->cnx) ) return;
-		if ( set_arg_float(gctx, 5, p->cny) ) return;
-		if ( set_arg_float(gctx, 6, p->res) ) return;
-		if ( set_arg_float(gctx, 7, p->clen) ) return;
-		if ( set_arg_int(gctx, 10, p->min_fs) ) return;
-		if ( set_arg_int(gctx, 11, p->min_ss) ) return;
-		if ( set_arg_float(gctx, 19, p->fsx) ) return;
-		if ( set_arg_float(gctx, 20, p->fsy) ) return;
-		if ( set_arg_float(gctx, 21, p->ssx) ) return;
-		if ( set_arg_float(gctx, 22, p->ssy) ) return;
-
-		dims[0] = pan_width * SAMPLING;
-		dims[1] = pan_height * SAMPLING;
-		dims[2] = BWSAMPLING * DIVSAMPLING * DIVSAMPLING;
-
-		err = clEnqueueNDRangeKernel(gctx->cq, gctx->kern, 3, NULL,
-		                             dims, ldims, 0, NULL, NULL);
-		if ( err != CL_SUCCESS ) {
-			ERROR("Couldn't enqueue diffraction kernel: %s\n",
-			      clError(err));
-			return;
-		}
-
-		clFinish(gctx->cq);
-
-		diff_ptr = clEnqueueMapBuffer(gctx->cq, diff, CL_TRUE,
-		                              CL_MAP_READ, 0, diff_size,
-		                              0, NULL, NULL, &err);
-		if ( err != CL_SUCCESS ) {
-			ERROR("Couldn't map diffraction buffer: %s\n",
-			      clError(err));
-			return;
-		}
-		tt_ptr = clEnqueueMapBuffer(gctx->cq, tt, CL_TRUE, CL_MAP_READ,
-		                            0, tt_size, 0, NULL, NULL, &err);
-		if ( err != CL_SUCCESS ) {
-			ERROR("Couldn't map tt buffer\n");
-			return;
-		}
-
-		for ( fs=0; fs<pan_width; fs++ ) {
-		for ( ss=0; ss<pan_height; ss++ ) {
-
-			float val, tt;
-			int tfs, tss;
-
-			val = diff_ptr[fs + pan_width*ss];
-			if ( isinf(val) ) n_inf++;
-			if ( val < 0.0 ) n_neg++;
-			if ( isnan(val) ) n_nan++;
-			tt = tt_ptr[fs + pan_width*ss];
-
-			tfs = p->min_fs + fs;
-			tss = p->min_ss + ss;
-			image->data[tfs + image->width*tss] = val;
-			image->twotheta[tfs + image->width*tss] = tt;
-
-		}
-		}
-
-		clEnqueueUnmapMemObject(gctx->cq, diff, diff_ptr,
-		                        0, NULL, NULL);
-		clEnqueueUnmapMemObject(gctx->cq, tt, tt_ptr,
-		                        0, NULL, NULL);
-
-		clReleaseMemObject(diff);
-		clReleaseMemObject(tt);
-
-	}
-
-
-	if ( n_neg + n_inf + n_nan ) {
-		ERROR("WARNING: The GPU calculation produced %i negative"
-		      " values, %i infinities and %i NaNs.\n",
-		      n_neg, n_inf, n_nan);
-	}
-
-}
-
-
-/* Setup the OpenCL stuff, create buffers, load the structure factor table */
-struct gpu_context *setup_gpu(int no_sfac,
-                              const double *intensities, unsigned char *flags,
-                              const char *sym, int dev_num)
-{
-	struct gpu_context *gctx;
-	cl_uint nplat;
-	cl_platform_id platforms[8];
-	cl_context_properties prop[3];
-	cl_int err;
-	cl_device_id dev;
-	size_t intensities_size;
-	float *intensities_ptr;
-	size_t flags_size;
-	float *flags_ptr;
-	size_t maxwgsize;
-	int i;
-	char cflags[512] = "";
-
-	STATUS("Setting up GPU...\n");
-
-	err = clGetPlatformIDs(8, platforms, &nplat);
-	if ( err != CL_SUCCESS ) {
-		ERROR("Couldn't get platform IDs: %i\n", err);
-		return NULL;
-	}
-	if ( nplat == 0 ) {
-		ERROR("Couldn't find at least one platform!\n");
-		return NULL;
-	}
-	prop[0] = CL_CONTEXT_PLATFORM;
-	prop[1] = (cl_context_properties)platforms[0];
-	prop[2] = 0;
-
-	gctx = malloc(sizeof(*gctx));
-	gctx->ctx = clCreateContextFromType(prop, CL_DEVICE_TYPE_GPU,
-	                                    NULL, NULL, &err);
-	if ( err != CL_SUCCESS ) {
-		ERROR("Couldn't create OpenCL context: %i\n", err);
-		free(gctx);
-		return NULL;
-	}
-
-	dev = get_cl_dev(gctx->ctx, dev_num);
-
-	gctx->cq = clCreateCommandQueue(gctx->ctx, dev, 0, &err);
-	if ( err != CL_SUCCESS ) {
-		ERROR("Couldn't create OpenCL command queue\n");
-		free(gctx);
-		return NULL;
-	}
-
-	/* Create a single-precision version of the scattering factors */
-	intensities_size = IDIM*IDIM*IDIM*sizeof(cl_float);
-	intensities_ptr = malloc(intensities_size);
-	if ( intensities != NULL ) {
-		for ( i=0; i<IDIM*IDIM*IDIM; i++ ) {
-			intensities_ptr[i] = intensities[i];
-		}
-	} else {
-		for ( i=0; i<IDIM*IDIM*IDIM; i++ ) {
-			intensities_ptr[i] = 1e5;
-		}
-		strncat(cflags, "-DFLAT_INTENSITIES ", 511-strlen(cflags));
-	}
-	gctx->intensities = clCreateBuffer(gctx->ctx,
-	                             CL_MEM_READ_ONLY | CL_MEM_COPY_HOST_PTR,
-	                             intensities_size, intensities_ptr, &err);
-	if ( err != CL_SUCCESS ) {
-		ERROR("Couldn't allocate intensities memory\n");
-		free(gctx);
-		return NULL;
-	}
-	free(intensities_ptr);
-
-	if ( sym != NULL ) {
-		if ( strcmp(sym, "1") == 0 ) {
-			strncat(cflags, "-DPG1 ", 511-strlen(cflags));
-		} else if ( strcmp(sym, "-1") == 0 ) {
-			strncat(cflags, "-DPG1BAR ", 511-strlen(cflags));
-		} else if ( strcmp(sym, "6/mmm") == 0 ) {
-			strncat(cflags, "-DPG6MMM ", 511-strlen(cflags));
-		} else if ( strcmp(sym, "6") == 0 ) {
-			strncat(cflags, "-DPG6 ", 511-strlen(cflags));
-		} else if ( strcmp(sym, "6/m") == 0 ) {
-			strncat(cflags, "-DPG6M ", 511-strlen(cflags));
-		} else {
-			ERROR("Sorry!  Point group '%s' is not currently"
-			      " supported on the GPU."
-			      " I'm using '1' instead.\n", sym);
-			strncat(cflags, "-DPG1 ", 511-strlen(cflags));
-		}
-	} else {
-		if ( intensities != NULL ) {
-			ERROR("You gave me an intensities file but no point"
-			      " group.  I'm assuming '1'.\n");
-			strncat(cflags, "-DPG1 ", 511-strlen(cflags));
-		}
-	}
-
-	/* Create a flag array */
-	flags_size = IDIM*IDIM*IDIM*sizeof(cl_float);
-	flags_ptr = malloc(flags_size);
-	if ( flags != NULL ) {
-		for ( i=0; i<IDIM*IDIM*IDIM; i++ ) {
-			flags_ptr[i] = flags[i];
-		}
-	} else {
-		for ( i=0; i<IDIM*IDIM*IDIM; i++ ) {
-			flags_ptr[i] = 1.0;
-		}
-	}
-	gctx->flags = clCreateBuffer(gctx->ctx,
-	                             CL_MEM_READ_ONLY | CL_MEM_COPY_HOST_PTR,
-	                             flags_size, flags_ptr, &err);
-	if ( err != CL_SUCCESS ) {
-		ERROR("Couldn't allocate flag buffer\n");
-		free(gctx);
-		return NULL;
-	}
-	free(flags_ptr);
-
-	gctx->prog = load_program(DATADIR"/crystfel/diffraction.cl", gctx->ctx,
-	                          dev, &err, cflags);
-	if ( err != CL_SUCCESS ) {
-		free(gctx);
-		return NULL;
-	}
-
-	gctx->kern = clCreateKernel(gctx->prog, "diffraction", &err);
-	if ( err != CL_SUCCESS ) {
-		ERROR("Couldn't create kernel\n");
-		free(gctx);
-		return NULL;
-	}
-
-	gctx->max_sinc_lut = 0;
-	gctx->sinc_lut_ptrs = NULL;
-	gctx->sinc_luts = NULL;
-
-	clGetDeviceInfo(dev, CL_DEVICE_MAX_WORK_GROUP_SIZE,
-	                sizeof(size_t), &maxwgsize, NULL);
-	STATUS("Maximum work group size = %lli\n", (long long int)maxwgsize);
-
-	return gctx;
-}
-
-
-void cleanup_gpu(struct gpu_context *gctx)
-{
-	int i;
-
-	clReleaseProgram(gctx->prog);
-	clReleaseMemObject(gctx->intensities);
-
-	/* Release LUTs */
-	for ( i=1; i<=gctx->max_sinc_lut; i++ ) {
-		if ( gctx->sinc_lut_ptrs[i-1] != NULL ) {
-			clReleaseMemObject(gctx->sinc_luts[i-1]);
-			free(gctx->sinc_lut_ptrs[i-1]);
-		}
-	}
-
-	free(gctx->sinc_luts);
-	free(gctx->sinc_lut_ptrs);
-
-	clReleaseCommandQueue(gctx->cq);
-	clReleaseContext(gctx->ctx);
-	free(gctx);
-}
diff --git a/src/diffraction-gpu.h b/src/diffraction-gpu.h
deleted file mode 100644
index a3bde4e1..00000000
--- a/src/diffraction-gpu.h
+++ /dev/null
@@ -1,57 +0,0 @@
-/*
- * diffraction-gpu.h
- *
- * Calculate diffraction patterns by Fourier methods (GPU version)
- *
- * (c) 2006-2011 Thomas White <taw@physics.org>
- *
- * Part of CrystFEL - crystallography with a FEL
- *
- */
-
-#ifdef HAVE_CONFIG_H
-#include <config.h>
-#endif
-
-#ifndef DIFFRACTION_GPU_H
-#define DIFFRACTION_GPU_H
-
-#include "image.h"
-#include "cell.h"
-
-struct gpu_context;
-
-#if HAVE_OPENCL
-
-extern void get_diffraction_gpu(struct gpu_context *gctx, struct image *image,
-                                int na, int nb, int nc, UnitCell *ucell);
-extern struct gpu_context *setup_gpu(int no_sfac,
-                                     const double *intensities,
-                                     const unsigned char *flags,
-                                     const char *sym, int dev_num);
-extern void cleanup_gpu(struct gpu_context *gctx);
-
-#else
-
-static void get_diffraction_gpu(struct gpu_context *gctx, struct image *image,
-                                int na, int nb, int nc, UnitCell *ucell)
-{
-	/* Do nothing */
-	ERROR("This copy of CrystFEL was not compiled with OpenCL support.\n");
-}
-
-static struct gpu_context *setup_gpu(int no_sfac,
-                                     const double *intensities,
-                                     const unsigned char *flags,
-                                     const char *sym, int dev_num)
-{
-	return NULL;
-}
-
-static void cleanup_gpu(struct gpu_context *gctx)
-{
-}
-
-#endif
-
-#endif	/* DIFFRACTION_GPU_H */
diff --git a/src/diffraction.c b/src/diffraction.c
deleted file mode 100644
index 9532a6ce..00000000
--- a/src/diffraction.c
+++ /dev/null
@@ -1,463 +0,0 @@
-/*
- * diffraction.c
- *
- * Calculate diffraction patterns by Fourier methods
- *
- * (c) 2006-2011 Thomas White <taw@physics.org>
- *
- * Part of CrystFEL - crystallography with a FEL
- *
- */
-
-
-#include <stdlib.h>
-#include <math.h>
-#include <stdio.h>
-#include <string.h>
-#include <complex.h>
-#include <assert.h>
-#include <fenv.h>
-
-#include "image.h"
-#include "utils.h"
-#include "cell.h"
-#include "diffraction.h"
-#include "beam-parameters.h"
-#include "symmetry.h"
-
-
-#define SAMPLING (4)
-#define BWSAMPLING (10)
-#define DIVSAMPLING (1)
-#define SINC_LUT_ELEMENTS (4096)
-
-
-static double *get_sinc_lut(int n)
-{
-	int i;
-	double *lut;
-
-	lut = malloc(SINC_LUT_ELEMENTS*sizeof(double));
-	lut[0] = n;
-	if ( n == 1 ) {
-		for ( i=1; i<SINC_LUT_ELEMENTS; i++ ) {
-			lut[i] = 1.0;
-		}
-	} else {
-		for ( i=1; i<SINC_LUT_ELEMENTS; i++ ) {
-			double x, val;
-			x = (double)i/SINC_LUT_ELEMENTS;
-			val = fabs(sin(M_PI*n*x)/sin(M_PI*x));
-			lut[i] = val;
-		}
-	}
-
-	return lut;
-}
-
-
-static double interpolate_lut(double *lut, double val)
-{
-	double i, pos, f;
-	unsigned int low, high;
-
-	pos = SINC_LUT_ELEMENTS * modf(fabs(val), &i);
-	low = (int)pos;  /* Discard fractional part */
-	high = low + 1;
-	f = modf(pos, &i);  /* Fraction */
-	if ( high == SINC_LUT_ELEMENTS ) high = 0;
-
-	return (1.0-f)*lut[low] + f*lut[high];
-}
-
-
-static double lattice_factor(struct rvec q, double ax, double ay, double az,
-                                            double bx, double by, double bz,
-                                            double cx, double cy, double cz,
-                                            double *lut_a, double *lut_b,
-                                            double *lut_c)
-{
-	struct rvec Udotq;
-	double f1, f2, f3;
-
-	Udotq.u = ax*q.u + ay*q.v + az*q.w;
-	Udotq.v = bx*q.u + by*q.v + bz*q.w;
-	Udotq.w = cx*q.u + cy*q.v + cz*q.w;
-
-	f1 = interpolate_lut(lut_a, Udotq.u);
-	f2 = interpolate_lut(lut_b, Udotq.v);
-	f3 = interpolate_lut(lut_c, Udotq.w);
-
-	return f1 * f2 * f3;
-}
-
-
-static double sym_lookup_intensity(const double *intensities,
-                                   const unsigned char *flags,
-                                   const SymOpList *sym,
-                                   signed int h, signed int k, signed int l)
-{
-	int i;
-	double ret = 0.0;
-
-	for ( i=0; i<num_equivs(sym, NULL); i++ ) {
-
-		signed int he;
-		signed int ke;
-		signed int le;
-		double f, val;
-
-		get_equiv(sym, NULL, i, h, k, l, &he, &ke, &le);
-
-		f = (double)lookup_flag(flags, he, ke, le);
-		val = lookup_intensity(intensities, he, ke, le);
-
-		ret += f*val;
-
-	}
-
-	return ret;
-}
-
-
-static double sym_lookup_phase(const double *phases,
-                               const unsigned char *flags, const SymOpList *sym,
-                               signed int h, signed int k, signed int l)
-{
-	int i;
-	double ret = 0.0;
-
-	for ( i=0; i<num_equivs(sym, NULL); i++ ) {
-
-		signed int he;
-		signed int ke;
-		signed int le;
-		double f, val;
-
-		get_equiv(sym, NULL, i, h, k, l, &he, &ke, &le);
-
-		f = (double)lookup_flag(flags, he, ke, le);
-		val = lookup_phase(phases, he, ke, le);
-
-		ret += f*val;
-
-	}
-
-	return ret;
-}
-
-
-static double interpolate_linear(const double *ref, const unsigned char *flags,
-                                 const SymOpList *sym, float hd,
-                                 signed int k, signed int l)
-{
-	signed int h;
-	double val1, val2;
-	float f;
-
-	h = (signed int)hd;
-	if ( hd < 0.0 ) h -= 1;
-	f = hd - (float)h;
-	assert(f >= 0.0);
-
-	val1 = sym_lookup_intensity(ref, flags, sym, h, k, l);
-	val2 = sym_lookup_intensity(ref, flags, sym, h+1, k, l);
-
-	val1 = val1;
-	val2 = val2;
-
-	return (1.0-f)*val1 + f*val2;
-}
-
-
-static double interpolate_bilinear(const double *ref,
-                                   const unsigned char *flags,
-                                   const SymOpList *sym,
-                                   float hd, float kd, signed int l)
-{
-	signed int k;
-	double val1, val2;
-	float f;
-
-	k = (signed int)kd;
-	if ( kd < 0.0 ) k -= 1;
-	f = kd - (float)k;
-	assert(f >= 0.0);
-
-	val1 = interpolate_linear(ref, flags, sym, hd, k, l);
-	val2 = interpolate_linear(ref, flags, sym, hd, k+1, l);
-
-	return (1.0-f)*val1 + f*val2;
-}
-
-
-static double interpolate_intensity(const double *ref,
-                                    const unsigned char *flags,
-                                    const SymOpList *sym,
-                                    float hd, float kd, float ld)
-{
-	signed int l;
-	double val1, val2;
-	float f;
-
-	l = (signed int)ld;
-	if ( ld < 0.0 ) l -= 1;
-	f = ld - (float)l;
-	assert(f >= 0.0);
-
-	val1 = interpolate_bilinear(ref, flags, sym, hd, kd, l);
-	val2 = interpolate_bilinear(ref, flags, sym, hd, kd, l+1);
-
-	return (1.0-f)*val1 + f*val2;
-}
-
-
-static double complex interpolate_phased_linear(const double *ref,
-                                                const double *phases,
-                                                const unsigned char *flags,
-                                                const SymOpList *sym,
-                                                float hd,
-                                                signed int k, signed int l)
-{
-	signed int h;
-	double val1, val2;
-	float f;
-	double ph1, ph2;
-	double re1, re2, im1, im2;
-	double re, im;
-
-	h = (signed int)hd;
-	if ( hd < 0.0 ) h -= 1;
-	f = hd - (float)h;
-	assert(f >= 0.0);
-
-	val1 = sym_lookup_intensity(ref, flags, sym, h, k, l);
-	val2 = sym_lookup_intensity(ref, flags, sym, h+1, k, l);
-	ph1 = sym_lookup_phase(phases, flags, sym, h, k, l);
-	ph2 = sym_lookup_phase(phases, flags, sym, h+1, k, l);
-
-	val1 = val1;
-	val2 = val2;
-
-	/* Calculate real and imaginary parts */
-	re1 = val1 * cos(ph1);
-	im1 = val1 * sin(ph1);
-	re2 = val2 * cos(ph2);
-	im2 = val2 * sin(ph2);
-
-	re = (1.0-f)*re1 + f*re2;
-	im = (1.0-f)*im1 + f*im2;
-
-	return re + im*I;
-}
-
-
-static double complex interpolate_phased_bilinear(const double *ref,
-                                                  const double *phases,
-                                                  const unsigned char *flags,
-                                                  const SymOpList *sym,
-                                                  float hd, float kd,
-                                                  signed int l)
-{
-	signed int k;
-	double complex val1, val2;
-	float f;
-
-	k = (signed int)kd;
-	if ( kd < 0.0 ) k -= 1;
-	f = kd - (float)k;
-	assert(f >= 0.0);
-
-	val1 = interpolate_phased_linear(ref, phases, flags, sym, hd, k, l);
-	val2 = interpolate_phased_linear(ref, phases, flags, sym, hd, k+1, l);
-
-	return (1.0-f)*val1 + f*val2;
-}
-
-
-static double interpolate_phased_intensity(const double *ref,
-                                           const double *phases,
-                                           const unsigned char *flags,
-                                           const SymOpList *sym,
-                                           float hd, float kd, float ld)
-{
-	signed int l;
-	double complex val1, val2;
-	float f;
-
-	l = (signed int)ld;
-	if ( ld < 0.0 ) l -= 1;
-	f = ld - (float)l;
-	assert(f >= 0.0);
-
-	val1 = interpolate_phased_bilinear(ref, phases, flags, sym,
-	                                   hd, kd, l);
-	val2 = interpolate_phased_bilinear(ref, phases, flags, sym,
-	                                   hd, kd, l+1);
-
-	return cabs((1.0-f)*val1 + f*val2);
-}
-
-
-/* Look up the structure factor for the nearest Bragg condition */
-static double molecule_factor(const double *intensities, const double *phases,
-                              const unsigned char *flags, struct rvec q,
-                              double ax, double ay, double az,
-                              double bx, double by, double bz,
-                              double cx, double cy, double cz,
-                              GradientMethod m, const SymOpList *sym)
-{
-	float hd, kd, ld;
-	signed int h, k, l;
-	double r;
-
-	hd = q.u * ax + q.v * ay + q.w * az;
-	kd = q.u * bx + q.v * by + q.w * bz;
-	ld = q.u * cx + q.v * cy + q.w * cz;
-
-	/* No flags -> flat intensity distribution */
-	if ( flags == NULL ) return 1.0e5;
-
-	switch ( m ) {
-	case GRADIENT_MOSAIC :
-		fesetround(1);  /* Round to nearest */
-		h = (signed int)rint(hd);
-		k = (signed int)rint(kd);
-		l = (signed int)rint(ld);
-		if ( abs(h) > INDMAX ) r = 0.0;
-		else if ( abs(k) > INDMAX ) r = 0.0;
-		else if ( abs(l) > INDMAX ) r = 0.0;
-		else r = sym_lookup_intensity(intensities, flags, sym, h, k, l);
-		break;
-	case GRADIENT_INTERPOLATE :
-		r = interpolate_intensity(intensities, flags, sym, hd, kd, ld);
-		break;
-	case GRADIENT_PHASED :
-		r = interpolate_phased_intensity(intensities, phases, flags,
-		                                 sym, hd, kd, ld);
-		break;
-	default:
-		ERROR("This gradient method not implemented yet.\n");
-		exit(1);
-	}
-
-	return r;
-}
-
-
-void get_diffraction(struct image *image, int na, int nb, int nc,
-                     const double *intensities, const double *phases,
-                     const unsigned char *flags, UnitCell *cell,
-                     GradientMethod m, const SymOpList *sym)
-{
-	unsigned int fs, ss;
-	double ax, ay, az;
-	double bx, by, bz;
-	double cx, cy, cz;
-	float klow, khigh, bwstep;
-	double *lut_a;
-	double *lut_b;
-	double *lut_c;
-	double divxlow, divylow, divxstep, divystep;
-
-	cell_get_cartesian(cell, &ax, &ay, &az, &bx, &by, &bz, &cx, &cy, &cz);
-
-	/* Allocate (and zero) the "diffraction array" */
-	image->data = calloc(image->width * image->height, sizeof(float));
-
-	/* Needed later for Lorentz calculation */
-	image->twotheta = malloc(image->width * image->height * sizeof(double));
-
-	klow = 1.0/(image->lambda*(1.0 + image->beam->bandwidth/2.0));
-	khigh = 1.0/(image->lambda*(1.0 - image->beam->bandwidth/2.0));
-	bwstep = (khigh-klow) / BWSAMPLING;
-
-	divxlow = -image->beam->divergence/2.0;
-	divylow = -image->beam->divergence/2.0;
-	divxstep = image->beam->divergence / DIVSAMPLING;
-	divystep = image->beam->divergence / DIVSAMPLING;
-
-	lut_a = get_sinc_lut(na);
-	lut_b = get_sinc_lut(nb);
-	lut_c = get_sinc_lut(nc);
-
-	for ( fs=0; fs<image->width; fs++ ) {
-	for ( ss=0; ss<image->height; ss++ ) {
-
-		int fs_step, ss_step, kstep;
-		int divxval, divyval;
-		int idx = fs + image->width*ss;
-
-		for ( fs_step=0; fs_step<SAMPLING; fs_step++ ) {
-		for ( ss_step=0; ss_step<SAMPLING; ss_step++ ) {
-		for ( kstep=0; kstep<BWSAMPLING; kstep++ ) {
-		for ( divxval=0; divxval<DIVSAMPLING; divxval++ ) {
-		for ( divyval=0; divyval<DIVSAMPLING; divyval++ ) {
-
-			double k;
-			double intensity;
-			double f_lattice, I_lattice;
-			double I_molecule;
-			struct rvec q, qn;
-			double twotheta;
-			const double dfs = (double)fs
-			                    + ((double)fs_step / SAMPLING);
-			const double dss = (double)ss
-			                    + ((double)ss_step / SAMPLING);
-
-			double xdiv = divxlow + divxstep*(double)divxval;
-			double ydiv = divylow + divystep*(double)divyval;
-
-			/* Calculate k this time round */
-			k = klow + (double)kstep * bwstep;
-
-			qn = get_q(image, dfs, dss, &twotheta, k);
-
-			/* x divergence */
-			q.u =  qn.u*cos(xdiv) +qn.w*sin(xdiv);
-			q.v =  qn.v;
-			q.w = -qn.u*sin(xdiv) +qn.w*cos(xdiv);
-
-			qn = q;
-
-			/* y divergence */
-			q.v =  qn.v*cos(ydiv) +qn.w*sin(ydiv);
-			q.w = -qn.v*sin(ydiv) +qn.w*cos(ydiv);
-
-			f_lattice = lattice_factor(q, ax, ay, az,
-			                              bx, by, bz,
-			                              cx, cy, cz,
-			                              lut_a, lut_b, lut_c);
-
-			I_molecule = molecule_factor(intensities,
-			                             phases, flags, q,
-			                             ax,ay,az,bx,by,bz,cx,cy,cz,
-			                             m, sym);
-
-			I_lattice = pow(f_lattice, 2.0);
-			intensity = I_lattice * I_molecule;
-
-			image->data[idx] += intensity;
-
-			if ( fs_step + ss_step + kstep == 0 ) {
-				image->twotheta[idx] = twotheta;
-			}
-
-		}
-		}
-		}
-		}
-		}
-
-		image->data[idx] /= (SAMPLING*SAMPLING*BWSAMPLING
-		                     *DIVSAMPLING*DIVSAMPLING);
-
-
-	}
-	progress_bar(fs, image->width-1, "Calculating diffraction");
-	}
-
-	free(lut_a);
-	free(lut_b);
-	free(lut_c);
-}
diff --git a/src/diffraction.h b/src/diffraction.h
deleted file mode 100644
index f71d3cce..00000000
--- a/src/diffraction.h
+++ /dev/null
@@ -1,34 +0,0 @@
-/*
- * diffraction.h
- *
- * Calculate diffraction patterns by Fourier methods
- *
- * (c) 2006-2011 Thomas White <taw@physics.org>
- *
- * Part of CrystFEL - crystallography with a FEL
- *
- */
-
-#ifdef HAVE_CONFIG_H
-#include <config.h>
-#endif
-
-#ifndef DIFFRACTION_H
-#define DIFFRACTION_H
-
-#include "image.h"
-#include "cell.h"
-#include "symmetry.h"
-
-typedef enum {
-	GRADIENT_MOSAIC,
-	GRADIENT_INTERPOLATE,
-	GRADIENT_PHASED
-} GradientMethod;
-
-extern void get_diffraction(struct image *image, int na, int nb, int nc,
-                            const double *intensities, const double *phases,
-                            const unsigned char *flags, UnitCell *cell,
-                            GradientMethod m, const SymOpList *sym);
-
-#endif	/* DIFFRACTION_H */
diff --git a/src/filters.c b/src/filters.c
deleted file mode 100644
index c4e409df..00000000
--- a/src/filters.c
+++ /dev/null
@@ -1,130 +0,0 @@
-/*
- * filters.c
- *
- * Image filtering
- *
- * (c) 2006-2010 Thomas White <taw@physics.org>
- *
- * Part of CrystFEL - crystallography with a FEL
- *
- */
-
-
-#ifdef HAVE_CONFIG_H
-#include <config.h>
-#endif
-
-#include <stdlib.h>
-#include <stdio.h>
-#include <math.h>
-#include <string.h>
-#include <assert.h>
-#include <gsl/gsl_statistics_int.h>
-#ifdef GSL_FUDGE
-#include <gsl/gsl_blas.h>
-#endif
-
-#include "image.h"
-
-
-static int compare_vals(const void *ap, const void *bp)
-{
-	const signed int a = *(signed int *)ap;
-	const signed int b = *(signed int *)bp;
-
-	if ( a > b ) return 1;
-	if ( a < b ) return -1;
-	return 0;
-}
-
-
-static void clean_panel(struct image *image, int sx, int sy)
-{
-	int x, y;
-	const int s = sizeof(signed int);
-
-	for ( x=0; x<512; x++ ) {
-
-		signed int vals[128];
-		double m;
-
-		for ( y=0; y<128; y++ ) {
-			vals[y] = image->data[(x+sx)+(y+sy)*image->width];
-		}
-
-		qsort(&vals[0], 128, s, compare_vals);
-
-		m = gsl_stats_int_median_from_sorted_data(vals, 1, 128);
-
-		for ( y=0; y<128; y++ ) {
-			image->data[(x+sx)+(y+sy)*image->width] -= m;
-		}
-
-	}
-}
-
-
-/* Pre-processing to make life easier */
-void filter_cm(struct image *image)
-{
-	int px, py;
-
-	if ( (image->width != 1024) || (image->height != 1024) ) return;
-
-	for ( px=0; px<2; px++ ) {
-	for ( py=0; py<8; py++ ) {
-
-		clean_panel(image, 512*px, 128*py);
-
-	}
-	}
-
-}
-
-
-void filter_noise(struct image *image, float *old)
-{
-	int x, y;
-
-	for ( x=0; x<image->width; x++ ) {
-	for ( y=0; y<image->height; y++ ) {
-
-		int dx, dy;
-		int val = image->data[x+image->width*y];
-
-		if ( old != NULL ) old[x+image->width*y] = val;
-
-		/* FIXME: This isn't really the right thing to do
-		 * at the edges. */
-		if ( (x==0) || (x==image->width-1)
-		  || (y==0) || (y==image->height-1) ) {
-			if ( val < 0 ) val = 0;
-			continue;
-		}
-
-		for ( dx=-1; dx<=+1; dx++ ) {
-		for ( dy=-1; dy<=+1; dy++ ) {
-
-			int val2;
-
-			val2 = image->data[(x+dx)+image->width*(y+dy)];
-
-			if ( val2 < 0 ) val = 0;
-
-		}
-		}
-
-		image->data[x+image->width*y] = val;
-
-	}
-	}
-}
-
-
-#ifdef GSL_FUDGE
-/* Force the linker to bring in CBLAS to make GSL happy */
-void filters_fudge_gslcblas()
-{
-        STATUS("%p\n", cblas_sgemm);
-}
-#endif
diff --git a/src/filters.h b/src/filters.h
deleted file mode 100644
index 6b35d7e8..00000000
--- a/src/filters.h
+++ /dev/null
@@ -1,25 +0,0 @@
-/*
- * peaks.h
- *
- * Image filtering
- *
- * (c) 2006-2010 Thomas White <taw@physics.org>
- *
- * Part of CrystFEL - crystallography with a FEL
- *
- */
-
-
-#ifndef FILTERS_H
-#define FILTERS_H
-
-#ifdef HAVE_CONFIG_H
-#include <config.h>
-#endif
-
-
-extern void filter_cm(struct image *image);
-extern void filter_noise(struct image *image, float *old);
-
-
-#endif	/* FILTERS_H */
diff --git a/src/geometry.c b/src/geometry.c
deleted file mode 100644
index 485abba3..00000000
--- a/src/geometry.c
+++ /dev/null
@@ -1,341 +0,0 @@
-/*
- * geometry.c
- *
- * Geometry of diffraction
- *
- * (c) 2006-2011 Thomas White <taw@physics.org>
- *
- * Part of CrystFEL - crystallography with a FEL
- *
- */
-
-#ifdef HAVE_CONFIG_H
-#include <config.h>
-#endif
-
-
-#include <stdlib.h>
-#include <assert.h>
-
-#include "utils.h"
-#include "cell.h"
-#include "image.h"
-#include "peaks.h"
-#include "beam-parameters.h"
-#include "reflist.h"
-#include "reflist-utils.h"
-#include "symmetry.h"
-
-
-static signed int locate_peak(double x, double y, double z, double k,
-                              struct detector *det, double *xdap, double *ydap)
-{
-	int i;
-	signed int found = -1;
-	const double den = k + z;
-
-	*xdap = -1;  *ydap = -1;
-
-	for ( i=0; i<det->n_panels; i++ ) {
-
-		double xd, yd;
-		double fs, ss, plx, ply;
-		struct panel *p;
-
-		p = &det->panels[i];
-
-		/* Coordinates of peak relative to central beam, in m */
-		xd = p->clen * x / den;
-		yd = p->clen * y / den;
-
-		/* Convert to pixels */
-		xd *= p->res;
-		yd *= p->res;
-
-		/* Convert to relative to the panel corner */
-		plx = xd - p->cnx;
-		ply = yd - p->cny;
-
-		fs = p->xfs*plx + p->yfs*ply;
-		ss = p->xss*plx + p->yss*ply;
-
-		fs += p->min_fs;
-		ss += p->min_ss;
-
-		/* Now, is this on this panel? */
-		if ( fs < p->min_fs ) continue;
-		if ( fs > p->max_fs ) continue;
-		if ( ss < p->min_ss ) continue;
-		if ( ss > p->max_ss ) continue;
-
-		/* If peak appears on multiple panels, reject it */
-		if ( found != -1 ) return -1;
-
-		/* Woohoo! */
-		found = i;
-		*xdap = fs;
-		*ydap = ss;
-
-	}
-
-	return found;
-}
-
-
-static double excitation_error(double xl, double yl, double zl,
-                               double ds, double k, double divergence,
-                               double tt)
-{
-	double al;
-	double r;
-	double delta;
-
-	al = M_PI_2 - asin(-zl/ds);
-
-	r = ( ds * sin(al) / sin(tt) ) - k;
-
-	delta = sqrt(2.0 * pow(ds, 2.0) * (1.0-cos(divergence)));
-	if ( divergence > 0.0 ) {
-		r += delta;
-	} else {
-		r -= delta;
-	}
-
-	return r;
-}
-
-
-static double partiality(double r1, double r2, double r)
-{
-	double q1, q2;
-	double p1, p2;
-
-	/* Calculate degrees of penetration */
-	q1 = (r1 + r)/(2.0*r);
-	q2 = (r2 + r)/(2.0*r);
-
-	/* Convert to partiality */
-	p1 = 3.0*pow(q1,2.0) - 2.0*pow(q1,3.0);
-	p2 = 3.0*pow(q2,2.0) - 2.0*pow(q2,3.0);
-
-	return p2 - p1;
-}
-
-
-static Reflection *check_reflection(struct image *image,
-                                    signed int h, signed int k, signed int l,
-                                    double asx, double asy, double asz,
-                                    double bsx, double bsy, double bsz,
-                                    double csx, double csy, double csz)
-{
-	const int output = 0;
-	double xl, yl, zl;
-	double ds, ds_sq;
-	double rlow, rhigh;     /* "Excitation error" */
-	signed int p;           /* Panel number */
-	double xda, yda;        /* Position on detector */
-	int close, inside;
-	double part;            /* Partiality */
-	int clamp_low = 0;
-	int clamp_high = 0;
-	double bandwidth = image->bw;
-	double divergence = image->div;
-	double lambda = image->lambda;
-	double klow, kcen, khigh;    /* Wavenumber */
-	Reflection *refl;
-	double tt;
-
-	/* "low" gives the largest Ewald sphere,
-	 * "high" gives the smallest Ewald sphere. */
-	klow = 1.0/(lambda - lambda*bandwidth/2.0);
-	kcen = 1.0/lambda;
-	khigh = 1.0/(lambda + lambda*bandwidth/2.0);
-
-	/* Get the coordinates of the reciprocal lattice point */
-	zl = h*asz + k*bsz + l*csz;
-	/* Throw out if it's "in front".  A tiny bit "in front" is OK. */
-	if ( zl > image->profile_radius ) return NULL;
-	xl = h*asx + k*bsx + l*csx;
-	yl = h*asy + k*bsy + l*csy;
-
-	tt = angle_between(0.0, 0.0, 1.0,  xl, yl, zl+kcen);
-	if ( tt > deg2rad(90.0) ) return NULL;
-
-	ds_sq = modulus_squared(xl, yl, zl);  /* d*^2 */
-	ds = sqrt(ds_sq);
-
-	/* Calculate excitation errors */
-	rlow = excitation_error(xl, yl, zl, ds, klow, -divergence/2.0, tt);
-	rhigh = excitation_error(xl, yl, zl, ds, khigh, +divergence/2.0, tt);
-
-	/* Is the reciprocal lattice point close to either extreme of
-	 * the sphere, maybe just outside the "Ewald volume"? */
-	close = (fabs(rlow) < image->profile_radius)
-	     || (fabs(rhigh) < image->profile_radius);
-
-	/* Is the reciprocal lattice point somewhere between the
-	 * extremes of the sphere, i.e. inside the "Ewald volume"? */
-	inside = signbit(rlow) ^ signbit(rhigh);
-
-	/* Can't be both inside and close */
-	if ( inside ) close = 0;
-
-	/* Neither?  Skip it. */
-	if ( !(close || inside) ) return NULL;
-
-	/* If the "lower" Ewald sphere is a long way away, use the
-	 * position at which the Ewald sphere would just touch the
-	 * reflection. */
-	if ( rlow < -image->profile_radius ) {
-		rlow = -image->profile_radius;
-		clamp_low = -1;
-	}
-	if ( rlow > +image->profile_radius ) {
-		rlow = +image->profile_radius;
-		clamp_low = +1;
-	}
-	/* Likewise the "higher" Ewald sphere */
-	if ( rhigh < -image->profile_radius ) {
-		rhigh = -image->profile_radius;
-		clamp_high = -1;
-	}
-	if ( rhigh > +image->profile_radius ) {
-		rhigh = +image->profile_radius;
-		clamp_high = +1;
-	}
-	assert(clamp_low <= clamp_high);
-	/* The six possible combinations of clamp_{low,high} (including
-	 * zero) correspond to the six situations in Table 3 of Rossmann
-	 * et al. (1979). */
-
-	/* Calculate partiality */
-	part = partiality(rlow, rhigh, image->profile_radius);
-
-	/* Locate peak on detector. */
-	p = locate_peak(xl, yl, zl, kcen, image->det, &xda, &yda);
-	if ( p == -1 ) return NULL;
-
-	/* Add peak to list */
-	refl = reflection_new(h, k, l);
-	set_detector_pos(refl, 0.0, xda, yda);
-	set_partial(refl, rlow, rhigh, part, clamp_low, clamp_high);
-	set_symmetric_indices(refl, h, k, l);
-	set_redundancy(refl, 1);
-
-	if ( output ) {
-		printf("%3i %3i %3i %6f (at %5.2f,%5.2f) %5.2f\n",
-		       h, k, l, 0.0, xda, yda, part);
-	}
-
-	return refl;
-}
-
-
-RefList *find_intersections(struct image *image, UnitCell *cell)
-{
-	double asx, asy, asz;
-	double bsx, bsy, bsz;
-	double csx, csy, csz;
-	RefList *reflections;
-	int hmax, kmax, lmax;
-	double mres;
-	signed int h, k, l;
-
-	reflections = reflist_new();
-
-	/* Cell angle check from Foadi and Evans (2011) */
-	if ( !cell_is_sensible(cell) ) {
-		ERROR("Invalid unit cell parameters given to"
-		      " find_intersections()\n");
-		cell_print(cell);
-		return NULL;
-	}
-
-	cell_get_reciprocal(cell, &asx, &asy, &asz,
-	                          &bsx, &bsy, &bsz,
-	                          &csx, &csy, &csz);
-
-	/* We add a horrific 20% fudge factor because bandwidth, divergence
-	 * and so on mean reflections appear beyond the largest q */
-	mres = 1.2 * largest_q(image);
-
-	hmax = mres / modulus(asx, asy, asz);
-	kmax = mres / modulus(bsx, bsy, bsz);
-	lmax = mres / modulus(csx, csy, csz);
-
-	if ( (hmax >= 256) || (kmax >= 256) || (lmax >= 256) ) {
-		ERROR("Unit cell is stupidly large.\n");
-		cell_print(cell);
-		if ( hmax >= 256 ) hmax = 255;
-		if ( kmax >= 256 ) kmax = 255;
-		if ( lmax >= 256 ) lmax = 255;
-	}
-
-	for ( h=-hmax; h<=hmax; h++ ) {
-	for ( k=-kmax; k<=kmax; k++ ) {
-	for ( l=-lmax; l<=lmax; l++ ) {
-
-		Reflection *refl;
-
-		refl = check_reflection(image, h, k, l,
-		                        asx,asy,asz,bsx,bsy,bsz,csx,csy,csz);
-
-		if ( refl != NULL ) {
-			refl = add_refl_to_list(refl, reflections);
-		}
-
-	}
-	}
-	}
-
-	return reflections;
-}
-
-
-/* Calculate partialities and apply them to the image's reflections */
-void update_partialities(struct image *image)
-{
-	Reflection *refl;
-	RefListIterator *iter;
-	RefList *predicted;
-	double asx, asy, asz;
-	double bsx, bsy, bsz;
-	double csx, csy, csz;
-
-	cell_get_reciprocal(image->indexed_cell, &asx, &asy, &asz,
-	                    &bsx, &bsy, &bsz, &csx, &csy, &csz);
-
-	/* Scratch list to give check_reflection() something to add to */
-	predicted = reflist_new();
-
-	for ( refl = first_refl(image->reflections, &iter);
-	      refl != NULL;
-	      refl = next_refl(refl, iter) )
-	{
-		Reflection *vals;
-		double r1, r2, p, x, y;
-		signed int h, k, l;
-		int clamp1, clamp2;
-
-		get_symmetric_indices(refl, &h, &k, &l);
-
-		vals = check_reflection(image, h, k, l,
-		                        asx,asy,asz,bsx,bsy,bsz,csx,csy,csz);
-
-		if ( vals == NULL ) {
-			set_redundancy(refl, 0);
-			continue;
-		}
-		set_redundancy(refl, 1);
-
-		/* Transfer partiality stuff */
-		get_partial(vals, &r1, &r2, &p, &clamp1, &clamp2);
-		set_partial(refl, r1, r2, p, clamp1, clamp2);
-
-		/* Transfer detector location */
-		get_detector_pos(vals, &x, &y);
-		set_detector_pos(refl, 0.0, x, y);
-	}
-
-	reflist_free(predicted);
-}
diff --git a/src/geometry.h b/src/geometry.h
deleted file mode 100644
index ddf04b80..00000000
--- a/src/geometry.h
+++ /dev/null
@@ -1,26 +0,0 @@
-/*
- * geometry.h
- *
- * Geometry of diffraction
- *
- * (c) 2006-2010 Thomas White <taw@physics.org>
- *
- * Part of CrystFEL - crystallography with a FEL
- *
- */
-
-#ifndef GEOMETRY_H
-#define GEOMETRY_H
-
-
-#ifdef HAVE_CONFIG_H
-#include <config.h>
-#endif
-
-#include "reflist.h"
-
-extern RefList *find_intersections(struct image *image, UnitCell *cell);
-
-extern void update_partialities(struct image *image);
-
-#endif	/* GEOMETRY_H */
diff --git a/src/peaks.c b/src/peaks.c
deleted file mode 100644
index f1a58d23..00000000
--- a/src/peaks.c
+++ /dev/null
@@ -1,594 +0,0 @@
-/*
- * peaks.c
- *
- * Peak search and other image analysis
- *
- * (c) 2006-2011 Thomas White <taw@physics.org>
- *	    2011 Andrew Martin
- *
- * Part of CrystFEL - crystallography with a FEL
- *
- */
-
-
-#ifdef HAVE_CONFIG_H
-#include <config.h>
-#endif
-
-#include <stdlib.h>
-#include <stdio.h>
-#include <math.h>
-#include <string.h>
-#include <assert.h>
-#include <gsl/gsl_statistics_int.h>
-#include <pthread.h>
-#include <fenv.h>
-
-#include "image.h"
-#include "utils.h"
-#include "index.h"
-#include "peaks.h"
-#include "detector.h"
-#include "filters.h"
-#include "diffraction.h"
-#include "reflist-utils.h"
-#include "beam-parameters.h"
-
-
-/* How close a peak must be to an indexed position to be considered "close"
- * for the purposes of double hit detection and sanity checking. */
-#define PEAK_CLOSE (30.0)
-
-/* How close a peak must be to an indexed position to be considered "close"
- * for the purposes of integration. */
-#define PEAK_REALLY_CLOSE (10.0)
-
-/* Degree of polarisation of X-ray beam */
-#define POL (1.0)
-
-static int cull_peaks_in_panel(struct image *image, struct panel *p)
-{
-	int i, n;
-	int nelim = 0;
-
-	n = image_feature_count(image->features);
-
-	for ( i=0; i<n; i++ ) {
-
-		struct imagefeature *f;
-		int j, ncol;
-
-		f = image_get_feature(image->features, i);
-		if ( f == NULL ) continue;
-
-		if ( f->fs < p->min_fs ) continue;
-		if ( f->fs > p->max_fs ) continue;
-		if ( f->ss < p->min_ss ) continue;
-		if ( f->ss > p->max_ss ) continue;
-
-		/* How many peaks are in the same column? */
-		ncol = 0;
-		for ( j=0; j<n; j++ ) {
-
-			struct imagefeature *g;
-
-			if ( i==j ) continue;
-
-			g = image_get_feature(image->features, j);
-			if ( g == NULL ) continue;
-
-			if ( p->badrow == 'f' ) {
-				if ( fabs(f->ss - g->ss) < 2.0 ) ncol++;
-			} else if ( p->badrow == 's' ) {
-				if ( fabs(f->fs - g->fs) < 2.0 ) ncol++;
-			} /* else do nothing */
-
-		}
-
-		/* More than three? */
-		if ( ncol <= 3 ) continue;
-
-		/* Yes?  Delete them all... */
-		nelim = 0;
-		for ( j=0; j<n; j++ ) {
-			struct imagefeature *g;
-			g = image_get_feature(image->features, j);
-			if ( g == NULL ) continue;
-			if ( p->badrow == 'f' ) {
-				if ( fabs(f->ss - g->ss) < 2.0 ) {
-					image_remove_feature(image->features,
-					                     j);
-					nelim++;
-				}
-			} else if ( p->badrow == 's' ) {
-				if ( fabs(f->fs - g->ss) < 2.0 ) {
-					image_remove_feature(image->features,
-					                     j);
-					nelim++;
-				}
-			} else {
-				ERROR("Invalid badrow direction.\n");
-				abort();
-			}
-
-		}
-
-	}
-
-	return nelim;
-}
-
-
-/* Post-processing of the peak list to remove noise */
-static int cull_peaks(struct image *image)
-{
-	int nelim = 0;
-	struct panel *p;
-	int i;
-
-	for ( i=0; i<image->det->n_panels; i++ ) {
-		p = &image->det->panels[i];
-		if ( p->badrow != '-' ) {
-			nelim += cull_peaks_in_panel(image, p);
-		}
-	}
-
-	return nelim;
-}
-
-
-/* Returns non-zero if peak has been vetoed.
- * i.e. don't use result if return value is not zero. */
-int integrate_peak(struct image *image, int cfs, int css,
-                   double *pfs, double *pss, double *intensity,
-                   double *pbg, double *pmax, double *sigma,
-                   int do_polar, int centroid, int bgsub)
-{
-	signed int fs, ss;
-	double lim, out_lim, mid_lim;
-	double lim_sq, out_lim_sq, mid_lim_sq;
-	double total = 0.0;
-	double fsct = 0.0;
-	double ssct = 0.0;
-	double noise = 0.0;
-	int noise_counts = 0;
-	double max = 0.0;
-	struct panel *p = NULL;
-        int pixel_counts = 0;
-        double noise_mean = 0.0;
-	double noise_meansq = 0.0;
-	struct beam_params *beam;
-	double aduph;
-
-	beam = image->beam;
-	if ( beam != NULL ) {
-		aduph = image->beam->adu_per_photon;
-	} else {
-		aduph = 1.0;
-	}
-
-	p = find_panel(image->det, cfs, css);
-	if ( p == NULL ) return 1;
-	if ( p->no_index ) return 1;
-
-	lim = p->integr_radius;
-	mid_lim = 3.0 + lim;
-	out_lim = 6.0 + lim;
-	lim_sq = pow(lim, 2.0);
-	mid_lim_sq = pow(mid_lim, 2.0);
-	out_lim_sq = pow(out_lim, 2.0);
-
-	for ( fs=-out_lim; fs<+out_lim; fs++ ) {
-	for ( ss=-out_lim; ss<+out_lim; ss++ ) {
-
-		double val;
-		double tt = 0.0;
-		double phi, pa, pb, pol;
-		uint16_t flags;
-		struct panel *p2;
-		int idx;
-
-		/* Outer mask radius */
-		if ( fs*fs + ss*ss > out_lim_sq ) continue;
-
-		if ( ((fs+cfs)>=image->width) || ((fs+cfs)<0) ) continue;
-		if ( ((ss+css)>=image->height) || ((ss+css)<0) ) continue;
-
-		/* Strayed off one panel? */
-		p2 = find_panel(image->det, fs+cfs, ss+css);
-		if ( p2 != p ) return 1;
-
-		idx = fs+cfs+image->width*(ss+css);
-
-		/* Veto this peak if we tried to integrate in a bad region */
-		if ( image->flags != NULL ) {
-
-			flags = image->flags[idx];
-
-			/* It must have all the "good" bits to be valid */
-			if ( !((flags & image->det->mask_good)
-			                   == image->det->mask_good) ) return 1;
-
-			/* If it has any of the "bad" bits, reject */
-			if ( flags & image->det->mask_bad ) return 1;
-
-		}
-
-		val = image->data[idx];
-
-		if ( do_polar ) {
-
-			tt = get_tt(image, fs+cfs, ss+css);
-
-			phi = atan2(ss+css, fs+cfs);
-			pa = pow(sin(phi)*sin(tt), 2.0);
-			pb = pow(cos(tt), 2.0);
-			pol = 1.0 - 2.0*POL*(1-pa) + POL*(1.0+pb);
-
-			val /= pol;
-
-		}
-
-		if ( val > max ) max = val;
-
-		/* If outside inner mask, estimate noise from this region */
-		if ( fs*fs + ss*ss > mid_lim_sq ) {
-
-			/* Noise
-			 * noise and noise_meansq are both in photons (^2) */
-			noise += val / image->beam->adu_per_photon;
-			noise_counts++;
-			noise_meansq += pow(val, 2.0);
-
-		} else if ( fs*fs + ss*ss < lim_sq ) {
-
-			/* Peak */
-			pixel_counts++;
-			total += val;
-			fsct += val*(cfs+fs);
-			ssct += val*(css+ss);
-
-		}
-
-	}
-	}
-
-	noise_mean = noise / noise_counts;  /* photons */
-
-	/* The centroid is excitingly undefined if there is no intensity */
-	centroid = 0;
-	if ( centroid && (total != 0) ) {
-		*pfs = ((double)fsct / total) + 0.5;
-		*pss = ((double)ssct / total) + 0.5;
-	} else {
-		*pfs = (double)cfs + 0.5;
-		*pss = (double)css + 0.5;
-	}
-	if ( bgsub ) {
-		*intensity = total - aduph * pixel_counts*noise_mean; /* ADU */
-	} else {
-		*intensity = total;  /* ADU */
-	}
-
-	if ( in_bad_region(image->det, *pfs, *pss) ) return 1;
-
-	if ( sigma != NULL ) {
-
-		/* First term is standard deviation of background per pixel
-		 * sqrt(pixel_counts) - increase of error for integrated value
-		 * sqrt(2) - increase of error for background subtraction  */
-		*sigma = sqrt(noise_meansq/noise_counts-(noise_mean*noise_mean))
-		          * sqrt(2.0*pixel_counts) * aduph;
-
-	}
-
-	if ( pbg != NULL ) {
-		*pbg = aduph * (noise / noise_counts);
-	}
-	if ( pmax != NULL ) {
-		*pmax = max;
-	}
-
-	return 0;
-}
-
-
-static void search_peaks_in_panel(struct image *image, float threshold,
-                                  float min_gradient, float min_snr,
-                                  struct panel *p)
-{
-	int fs, ss, stride;
-	float *data;
-	double d;
-	int idx;
-	double f_fs = 0.0;
-	double f_ss = 0.0;
-	double intensity = 0.0;
-	double sigma = 0.0;
-	double pbg = 0.0;
-	double pmax = 0.0;
-	int nrej_dis = 0;
-	int nrej_pro = 0;
-	int nrej_fra = 0;
-	int nrej_bad = 0;
-	int nrej_snr = 0;
-	int nacc = 0;
-	int ncull;
-	const int pws = p->peak_sep/2;
-
-	data = image->data;
-	stride = image->width;
-
-	for ( fs = p->min_fs+1; fs <= p->max_fs-1; fs++ ) {
-	for ( ss = p->min_ss+1; ss <= p->max_ss-1; ss++ ) {
-
-		double dx1, dx2, dy1, dy2;
-		double dxs, dys;
-		double grad;
-		int mask_fs, mask_ss;
-		int s_fs, s_ss;
-		double max;
-		unsigned int did_something;
-		int r;
-
-		/* Overall threshold */
-		if ( data[fs+stride*ss] < threshold ) continue;
-
-		/* Get gradients */
-		dx1 = data[fs+stride*ss] - data[(fs+1)+stride*ss];
-		dx2 = data[(fs-1)+stride*ss] - data[fs+stride*ss];
-		dy1 = data[fs+stride*ss] - data[(fs+1)+stride*(ss+1)];
-		dy2 = data[fs+stride*(ss-1)] - data[fs+stride*ss];
-
-		/* Average gradient measurements from both sides */
-		dxs = ((dx1*dx1) + (dx2*dx2)) / 2;
-		dys = ((dy1*dy1) + (dy2*dy2)) / 2;
-
-		/* Calculate overall gradient */
-		grad = dxs + dys;
-
-		if ( grad < min_gradient ) continue;
-
-		mask_fs = fs;
-		mask_ss = ss;
-
-		do {
-
-			max = data[mask_fs+stride*mask_ss];
-			did_something = 0;
-
-			for ( s_ss=biggest(mask_ss-pws/2,
-			                   p->min_ss);
-			      s_ss<=smallest(mask_ss+pws/2,
-			                     p->max_ss);
-			      s_ss++ ) {
-			for ( s_fs=biggest(mask_fs-pws/2,
-			                   p->min_fs);
-			      s_fs<=smallest(mask_fs+pws/2,
-			                     p->max_fs);
-			      s_fs++ ) {
-
-				if ( data[s_fs+stride*s_ss] > max ) {
-					max = data[s_fs+stride*s_ss];
-					mask_fs = s_fs;
-					mask_ss = s_ss;
-					did_something = 1;
-				}
-
-			}
-			}
-
-			/* Abort if drifted too far from the foot point */
-			if ( distance(mask_fs, mask_ss, fs, ss) >
-			     p->peak_sep/2.0 )
-			{
-				break;
-			}
-
-		} while ( did_something );
-
-		/* Too far from foot point? */
-		if ( distance(mask_fs, mask_ss, fs, ss) > p->peak_sep/2.0 ) {
-			nrej_dis++;
-			continue;
-		}
-
-		/* Should be enforced by bounds used above.  Muppet check. */
-		assert(mask_fs <= p->max_fs);
-		assert(mask_ss <= p->max_ss);
-		assert(mask_fs >= p->min_fs);
-		assert(mask_ss >= p->min_ss);
-
-		/* Centroid peak and get better coordinates.
-		 * Don't bother doing polarisation/SA correction, because the
-		 * intensity of this peak is only an estimate at this stage. */
-		r = integrate_peak(image, mask_fs, mask_ss,
-		                   &f_fs, &f_ss, &intensity,
-		                   &pbg, &pmax, &sigma, 0, 1, 1);
-
-		if ( r ) {
-			/* Bad region - don't detect peak */
-			nrej_bad++;
-			continue;
-		}
-
-		/* It is possible for the centroid to fall outside the image */
-		if ( (f_fs < p->min_fs) || (f_fs > p->max_fs)
-		  || (f_ss < p->min_ss) || (f_ss > p->max_ss) ) {
-			nrej_fra++;
-			continue;
-		}
-
-		if (intensity/sigma < min_snr) {
-			nrej_snr++;
-			continue;
-		}
-
-		/* Check for a nearby feature */
-		image_feature_closest(image->features, f_fs, f_ss, &d, &idx);
-		if ( d < p->peak_sep/2.0 ) {
-			nrej_pro++;
-			continue;
-		}
-
-		/* Add using "better" coordinates */
-		image_add_feature(image->features, f_fs, f_ss, image, intensity,
-		                  NULL);
-		nacc++;
-
-	}
-	}
-
-	if ( image->det != NULL ) {
-		ncull = cull_peaks(image);
-		nacc -= ncull;
-	} else {
-		STATUS("Not culling peaks because I don't have a "
-		       "detector geometry file.\n");
-		ncull = 0;
-	}
-
-//	STATUS("%i accepted, %i box, %i proximity, %i outside panel, "
-//	       "%i in bad regions, %i with SNR < %g, %i badrow culled.\n",
-//	       nacc, nrej_dis, nrej_pro, nrej_fra, nrej_bad,
-//	       nrej_snr, min_snr, ncull);
-}
-
-
-void search_peaks(struct image *image, float threshold, float min_gradient,
-                  float min_snr)
-{
-	int i;
-
-	if ( image->features != NULL ) {
-		image_feature_list_free(image->features);
-	}
-	image->features = image_feature_list_new();
-
-	for ( i=0; i<image->det->n_panels; i++ ) {
-
-		struct panel *p = &image->det->panels[i];
-
-		if ( p->no_index ) continue;
-		search_peaks_in_panel(image, threshold, min_gradient, min_snr, p);
-
-	}
-}
-
-
-int peak_sanity_check(struct image *image)
-{
-
-	int i;
-	int n_feat = 0;
-	int n_sane = 0;
-	double ax, ay, az;
-	double bx, by, bz;
-	double cx, cy, cz;
-	double min_dist = 0.25;
-
-	/* Round towards nearest */
-	fesetround(1);
-
-	/* Cell basis vectors for this image */
-	cell_get_cartesian(image->indexed_cell, &ax, &ay, &az,
-	                                        &bx, &by, &bz,
-	                                        &cx, &cy, &cz);
-
-	/* Loop over peaks, checking proximity to nearest reflection */
-	for ( i=0; i<image_feature_count(image->features); i++ ) {
-
-		struct imagefeature *f;
-		struct rvec q;
-		double h,k,l,hd,kd,ld;
-
-		/* Assume all image "features" are genuine peaks */
-		f = image_get_feature(image->features, i);
-		if ( f == NULL ) continue;
-		n_feat++;
-
-		/* Reciprocal space position of found peak */
-		q = get_q(image, f->fs, f->ss, NULL, 1.0/image->lambda);
-
-		/* Decimal and fractional Miller indices of nearest
-		 * reciprocal lattice point */
-		hd = q.u * ax + q.v * ay + q.w * az;
-		kd = q.u * bx + q.v * by + q.w * bz;
-		ld = q.u * cx + q.v * cy + q.w * cz;
-		h = lrint(hd);
-		k = lrint(kd);
-		l = lrint(ld);
-
-		/* Check distance */
-		if ( (fabs(h - hd) < min_dist) && (fabs(k - kd) < min_dist)
-		  && (fabs(l - ld) < min_dist) )
-		{
-			n_sane++;
-			continue;
-		}
-
-	}
-
-	/* return 0 means fail test, return 1 means pass test */
-	// printf("%d out of %d peaks are \"sane\"\n",n_sane,n_feat);
-	if ( (float)n_sane / (float)n_feat < 0.5 ) return 0;
-
-	return 1;
-}
-
-
-/* Integrate the list of predicted reflections in "image" */
-void integrate_reflections(struct image *image, int polar, int use_closer,
-                           int bgsub)
-{
-	Reflection *refl;
-	RefListIterator *iter;
-
-	for ( refl = first_refl(image->reflections, &iter);
-	      refl != NULL;
-	      refl = next_refl(refl, iter) ) {
-
-		double fs, ss, intensity;
-		double d;
-		int idx;
-		double bg, max;
-		double sigma;
-		double pfs, pss;
-		int r;
-
-		get_detector_pos(refl, &pfs, &pss);
-
-		/* Is there a really close feature which was detected? */
-		if ( use_closer ) {
-
-			struct imagefeature *f;
-
-			if ( image->features != NULL ) {
-				f = image_feature_closest(image->features,
-					                  pfs, pss, &d, &idx);
-			} else {
-				f = NULL;
-			}
-			if ( (f != NULL) && (d < PEAK_REALLY_CLOSE) ) {
-
-				pfs = f->fs;
-				pss = f->ss;
-
-			}
-		}
-
-		r = integrate_peak(image, pfs, pss, &fs, &ss,
-		                   &intensity, &bg, &max, &sigma, polar, 0,
-		                   bgsub);
-
-		/* Record intensity and set redundancy to 1 on success */
-		if ( r == 0 ) {
-			set_int(refl, intensity);
-			set_esd_intensity(refl, sigma);
-			set_redundancy(refl, 1);
-		} else {
-			set_redundancy(refl, 0);
-		}
-
-	}
-}
diff --git a/src/peaks.h b/src/peaks.h
deleted file mode 100644
index 9d475ea9..00000000
--- a/src/peaks.h
+++ /dev/null
@@ -1,38 +0,0 @@
-/*
- * peaks.h
- *
- * Peak search and other image analysis
- *
- * (c) 2006-2010 Thomas White <taw@physics.org>
- *
- * Part of CrystFEL - crystallography with a FEL
- *
- */
-
-
-#ifndef PEAKS_H
-#define PEAKS_H
-
-#ifdef HAVE_CONFIG_H
-#include <config.h>
-#endif
-
-#include <pthread.h>
-
-#include "reflist.h"
-
-extern void search_peaks(struct image *image, float threshold,
-                         float min_gradient, float min_snr);
-
-extern void integrate_reflections(struct image *image,
-                                  int polar, int use_closer, int bgsub);
-
-extern int peak_sanity_check(struct image * image);
-
-/* Exported so it can be poked by integration_check */
-extern int integrate_peak(struct image *image, int cfs, int css,
-                          double *pfs, double *pss, double *intensity,
-                          double *pbg, double *pmax, double *sigma,
-                          int do_polar, int centroid, int bgsub);
-
-#endif	/* PEAKS_H */
diff --git a/src/reflist-utils.c b/src/reflist-utils.c
deleted file mode 100644
index b64e9979..00000000
--- a/src/reflist-utils.c
+++ /dev/null
@@ -1,502 +0,0 @@
-/*
- * reflist-utils.c
- *
- * Utilities to complement the core reflist.c
- *
- * (c) 2006-2011 Thomas White <taw@physics.org>
- *
- * Part of CrystFEL - crystallography with a FEL
- *
- */
-
-
-#include <stdio.h>
-#include <assert.h>
-
-
-#include "reflist.h"
-#include "cell.h"
-#include "utils.h"
-#include "reflist-utils.h"
-#include "symmetry.h"
-
-
-/**
- * SECTION:reflist-utils
- * @short_description: Reflection list utilities
- * @title: RefList utilities
- * @section_id:
- * @see_also:
- * @include: "reflist-utils.h"
- * @Image:
- *
- * There are some utility functions associated with the core %RefList.
- **/
-
-
-double *intensities_from_list(RefList *list)
-{
-	Reflection *refl;
-	RefListIterator *iter;
-	double *out = new_list_intensity();
-
-	for ( refl = first_refl(list, &iter);
-	      refl != NULL;
-	      refl = next_refl(refl, iter) ) {
-
-		signed int h, k, l;
-		double intensity = get_intensity(refl);
-
-		get_indices(refl, &h, &k, &l);
-
-		set_intensity(out, h, k, l, intensity);
-
-	}
-
-	return out;
-}
-
-
-double *phases_from_list(RefList *list)
-{
-	Reflection *refl;
-	RefListIterator *iter;
-	double *out = new_list_phase();
-
-	for ( refl = first_refl(list, &iter);
-	      refl != NULL;
-	      refl = next_refl(refl, iter) ) {
-
-		signed int h, k, l;
-		double phase = get_phase(refl, NULL);
-
-		get_indices(refl, &h, &k, &l);
-
-		set_phase(out, h, k, l, phase);
-
-	}
-
-	return out;
-
-}
-
-
-unsigned char *flags_from_list(RefList *list)
-{
-	Reflection *refl;
-	RefListIterator *iter;
-	unsigned char *out = new_list_flag();
-
-	for ( refl = first_refl(list, &iter);
-	      refl != NULL;
-	      refl = next_refl(refl, iter) ) {
-
-		signed int h, k, l;
-
-		get_indices(refl, &h, &k, &l);
-
-		set_flag(out, h, k, l, 1);
-
-	}
-
-	return out;
-
-}
-
-
-int check_list_symmetry(RefList *list, const SymOpList *sym)
-{
-	Reflection *refl;
-	RefListIterator *iter;
-	SymOpMask *mask;
-
-	mask = new_symopmask(sym);
-	if ( mask == NULL ) {
-		ERROR("Couldn't create mask for list symmetry check.\n");
-		return 1;
-	}
-
-	for ( refl = first_refl(list, &iter);
-	      refl != NULL;
-	      refl = next_refl(refl, iter) ) {
-
-		int j;
-		int found = 0;
-		signed int h, k, l;
-		int n;
-
-		get_indices(refl, &h, &k, &l);
-
-		special_position(sym, mask, h, k, l);
-		n = num_equivs(sym, mask);
-
-		for ( j=0; j<n; j++ ) {
-
-			signed int he, ke, le;
-			Reflection *f;
-
-			get_equiv(sym, mask, j, h, k, l, &he, &ke, &le);
-
-			f = find_refl(list, he, ke, le);
-			if ( f != NULL ) found++;
-
-		}
-
-		assert(found != 0);  /* That'd just be silly */
-		if ( found > 1 ) {
-
-			STATUS("Found %i %i %i: %i times:\n", h, k, l, found);
-
-			for ( j=0; j<n; j++ ) {
-
-				signed int he, ke, le;
-				Reflection *f;
-
-				get_equiv(sym, mask, j, h, k, l, &he, &ke, &le);
-
-				f = find_refl(list, he, ke, le);
-				if ( f != NULL ) {
-					STATUS("%3i %3i %3i\n", he, ke, le);
-				}
-
-			}
-			free_symopmask(mask);
-
-			return 1;  /* Symmetry is wrong! */
-		}
-
-	}
-
-	free_symopmask(mask);
-
-	return 0;
-}
-
-
-int find_equiv_in_list(RefList *list, signed int h, signed int k,
-                       signed int l, const SymOpList *sym, signed int *hu,
-                       signed int *ku, signed int *lu)
-{
-	int i;
-	int found = 0;
-
-	for ( i=0; i<num_equivs(sym, NULL); i++ ) {
-
-		signed int he, ke, le;
-		Reflection *f;
-		get_equiv(sym, NULL, i, h, k, l, &he, &ke, &le);
-		f = find_refl(list, he, ke, le);
-
-		/* There must only be one equivalent.  If there are more, it
-		 * indicates that the user lied about the input symmetry.
-		 * This situation should have been checked for earlier by
-		 * calling check_symmetry() with 'items' and 'mero'. */
-
-		if ( (f != NULL) && !found ) {
-			*hu = he;  *ku = ke;  *lu = le;
-			return 1;
-		}
-
-	}
-
-	return 0;
-}
-
-
-/**
- * write_reflections_to_file:
- * @fh: File handle to write to
- * @list: The reflection list to write
- * @cell: Unit cell to use for generating 1/d values, or NULL.
- *
- * This function writes the contents of @list to @fh, using @cell to generate
- * 1/d values to ease later processing.  If @cell is NULL, 1/d values will not
- * be included ('-' will be written in their place).
- *
- * Reflections which have a redundancy of zero will not be written.
- *
- * The resulting list can be read back with read_reflections_from_file().
- **/
-void write_reflections_to_file(FILE *fh, RefList *list, UnitCell *cell)
-{
-	Reflection *refl;
-	RefListIterator *iter;
-
-	fprintf(fh, "  h   k   l          I    phase   sigma(I) "
-		     " 1/d(nm^-1)  counts  fs/px  ss/px\n");
-
-	for ( refl = first_refl(list, &iter);
-	      refl != NULL;
-	      refl = next_refl(refl, iter) ) {
-
-		signed int h, k, l;
-		double intensity, esd_i, s, ph;
-		int red;
-		double fs, ss;
-		char res[16];
-		char phs[16];
-		int have_phase;
-
-		get_indices(refl, &h, &k, &l);
-		get_detector_pos(refl, &fs, &ss);
-		intensity = get_intensity(refl);
-		esd_i = get_esd_intensity(refl);
-		red = get_redundancy(refl);
-		ph = get_phase(refl, &have_phase);
-
-		/* Reflections with redundancy = 0 are not written */
-		if ( red == 0 ) continue;
-
-		if ( cell != NULL ) {
-			s = 2.0 * resolution(cell, h, k, l);
-			snprintf(res, 16, "%10.2f", s/1e9);
-		} else {
-			strcpy(res, "         -");
-		}
-
-		if ( have_phase ) {
-			snprintf(phs, 16, "%8.2f", rad2deg(ph));
-		} else {
-			strncpy(phs, "       -", 15);
-		}
-
-		fprintf(fh,
-		       "%3i %3i %3i %10.2f %s %10.2f  %s %7i %6.1f %6.1f\n",
-		       h, k, l, intensity, phs, esd_i, res, red,
-		       fs, ss);
-
-	}
-}
-
-
-/**
- * write_reflist:
- * @filename: Filename
- * @list: The reflection list to write
- * @cell: Unit cell to use for generating 1/d values, or NULL.
- *
- * This function writes the contents of @list to @file, using @cell to generate
- * 1/d values to ease later processing.  If @cell is NULL, 1/d values will not
- * be included ('-' will be written in their place).
- *
- * Reflections which have a redundancy of zero will not be written.
- *
- * The resulting list can be read back with read_reflections_from_file() or
- * read_reflections().
- *
- * This is a convenience function which simply opens @filename and then calls
- * write_reflections_to_file.
- *
- * Returns: zero on success, non-zero on failure.
- **/
-int write_reflist(const char *filename, RefList *list, UnitCell *cell)
-{
-	FILE *fh;
-
-	if ( filename == NULL ) {
-		fh = stdout;
-	} else {
-		fh = fopen(filename, "w");
-	}
-
-	if ( fh == NULL ) {
-		ERROR("Couldn't open output file '%s'.\n", filename);
-		return 1;
-	}
-
-	write_reflections_to_file(fh, list, cell);
-	fprintf(fh, REFLECTION_END_MARKER"\n");
-
-	fclose(fh);
-
-	return 0;
-}
-
-
-RefList *read_reflections_from_file(FILE *fh)
-{
-	char *rval = NULL;
-	int first = 1;
-	RefList *out;
-
-	out = reflist_new();
-
-	do {
-
-		char line[1024];
-		signed int h, k, l;
-		float intensity, sigma, fs, ss;
-		char phs[1024];
-		char ress[1024];
-		int cts;
-		int r;
-		Reflection *refl;
-
-		rval = fgets(line, 1023, fh);
-		if ( rval == NULL ) continue;
-		chomp(line);
-
-		if ( strcmp(line, REFLECTION_END_MARKER) == 0 ) return out;
-
-		r = sscanf(line, "%i %i %i %f %s %f %s %i %f %f",
-		           &h, &k, &l, &intensity, phs, &sigma, ress, &cts,
-		           &fs, &ss);
-		if ( (r != 10) && (!first) ) {
-			reflist_free(out);
-			return NULL;
-		}
-
-		first = 0;
-		if ( r == 10 ) {
-
-			double ph;
-			char *v;
-
-			refl = add_refl(out, h, k, l);
-			set_int(refl, intensity);
-			set_detector_pos(refl, 0.0, fs, ss);
-			set_esd_intensity(refl, sigma);
-			set_redundancy(refl, cts);
-
-			ph = strtod(phs, &v);
-			if ( v != NULL ) set_ph(refl, deg2rad(ph));
-
-			/* The 1/d value is actually ignored. */
-
-		}
-
-	} while ( rval != NULL );
-
-	/* Got read error of some kind before finding PEAK_LIST_END_MARKER */
-	return NULL;
-}
-
-
-RefList *read_reflections(const char *filename)
-{
-	FILE *fh;
-	RefList *out;
-
-	if ( filename == NULL ) {
-		fh = stdout;
-	} else {
-		fh = fopen(filename, "r");
-	}
-
-	if ( fh == NULL ) {
-		ERROR("Couldn't open input file '%s'.\n", filename);
-		return NULL;
-	}
-
-	out = read_reflections_from_file(fh);
-
-	fclose(fh);
-
-	return out;
-}
-
-
-/**
- * asymmetric_indices:
- * @in: A %RefList
- * @sym: A %SymOpList
- *
- * This function creates a newly allocated copy of @in, but indexed using the
- * asymmetric indices according to @sym instead of the original indices.  The
- * original indices are stored and can be retrieved using
- * get_symmetric_indices() if required.
- *
- * Returns: the new %RefList, or NULL on failure.
- **/
-RefList *asymmetric_indices(RefList *in, const SymOpList *sym)
-{
-	Reflection *refl;
-	RefListIterator *iter;
-	RefList *new;
-
-	new = reflist_new();
-	if ( new == NULL ) return NULL;
-
-	for ( refl = first_refl(in, &iter);
-	      refl != NULL;
-	      refl = next_refl(refl, iter) ) {
-
-		signed int h, k, l;
-		signed int ha, ka, la;
-		Reflection *cr;
-
-		get_indices(refl, &h, &k, &l);
-
-		get_asymm(sym, h, k, l, &ha, &ka, &la);
-
-		cr = add_refl(new, ha, ka, la);
-		assert(cr != NULL);
-
-		copy_data(cr, refl);
-		set_symmetric_indices(cr, h, k, l);
-
-	}
-
-	return new;
-}
-
-
-/**
- * resolution_limits:
- * @list: A %RefList
- * @cell: A %UnitCell
- * @rmin: Place to store the minimum 1/d value
- * @rmax: Place to store the maximum 1/d value
- *
- * This function calculates the minimum and maximum values of 1/d, where
- * 2dsin(theta) = wavelength.  The answers are in m^-1.
- **/
-void resolution_limits(RefList *list, UnitCell *cell,
-                       double *rmin, double *rmax)
-{
-	Reflection *refl;
-	RefListIterator *iter;
-
-	*rmin = INFINITY;
-	*rmax = 0.0;
-
-	for ( refl = first_refl(list, &iter);
-	      refl != NULL;
-	      refl = next_refl(refl, iter) )
-	{
-		double r;
-		signed int h, k, l;
-
-		get_indices(refl, &h, &k, &l);
-		r = 2.0 * resolution(cell, h, k, l);
-
-		if ( r > *rmax ) *rmax = r;
-		if ( r < *rmin ) *rmin = r;
-	}
-}
-
-
-/**
- * max_intensity:
- * @list: A %RefList
- *
- * Returns: The maximum intensity in @list.
- **/
-double max_intensity(RefList *list)
-{
-	Reflection *refl;
-	RefListIterator *iter;
-	double max;
-
-	max = -INFINITY;
-
-	for ( refl = first_refl(list, &iter);
-	      refl != NULL;
-	      refl = next_refl(refl, iter) )
-	{
-		double val = get_intensity(refl);
-		if ( val > max ) max = val;
-	}
-
-	return max;
-}
diff --git a/src/reflist-utils.h b/src/reflist-utils.h
deleted file mode 100644
index d14e5f8e..00000000
--- a/src/reflist-utils.h
+++ /dev/null
@@ -1,52 +0,0 @@
-/*
- * reflist-utils.h
- *
- * Utilities to complement the core reflist.c
- *
- * (c) 2006-2011 Thomas White <taw@physics.org>
- *
- * Part of CrystFEL - crystallography with a FEL
- *
- */
-
-#ifdef HAVE_CONFIG_H
-#include <config.h>
-#endif
-
-#ifndef REFLIST_UTILS_H
-#define REFLIST_UTILS_H
-
-
-#include "reflist.h"
-#include "cell.h"
-#include "symmetry.h"
-
-
-#define REFLECTION_END_MARKER "End of reflections"
-
-
-extern void write_reflections_to_file(FILE *fh, RefList *list, UnitCell *cell);
-
-extern int write_reflist(const char *filename, RefList *list, UnitCell *cell);
-
-extern RefList *read_reflections_from_file(FILE *fh);
-
-extern RefList *read_reflections(const char *filename);
-
-extern double *intensities_from_list(RefList *list);
-extern double *phases_from_list(RefList *list);
-extern unsigned char *flags_from_list(RefList *list);
-
-extern int check_list_symmetry(RefList *list, const SymOpList *sym);
-extern int find_equiv_in_list(RefList *list, signed int h, signed int k,
-                              signed int l, const SymOpList *sym, signed int *hu,
-                              signed int *ku, signed int *lu);
-
-extern RefList *asymmetric_indices(RefList *in, const SymOpList *sym);
-
-extern void resolution_limits(RefList *list, UnitCell *cell,
-                              double *rmin, double *rmax);
-
-extern double max_intensity(RefList *list);
-
-#endif	/* REFLIST_UTILS_H */
diff --git a/src/statistics.c b/src/statistics.c
deleted file mode 100644
index 7990672d..00000000
--- a/src/statistics.c
+++ /dev/null
@@ -1,668 +0,0 @@
-/*
- * statistics.c
- *
- * Structure-factor statistics
- *
- * (c) 2006-2010 Thomas White <taw@physics.org>
- *
- * Part of CrystFEL - crystallography with a FEL
- *
- */
-
-
-#ifdef HAVE_CONFIG_H
-#include <config.h>
-#endif
-
-#include <math.h>
-#include <stdlib.h>
-#include <gsl/gsl_errno.h>
-#include <gsl/gsl_min.h>
-#include <gsl/gsl_statistics.h>
-
-#include "statistics.h"
-#include "utils.h"
-
-/**
- * SECTION:statistics
- * @short_description: Intensity statistics and R-factors
- * @title: Statistics
- * @section_id:
- * @see_also:
- * @include: "statistics.h"
- * @Image:
- *
- * These functions are for calculating various figures of merit.
- */
-
-
-struct r_params {
-	RefList *list1;
-	RefList *list2;
-	int fom;              /* Which FoM to use (see the enum just below) */
-};
-
-enum {
-	R_1_ZERO,
-	R_1_IGNORE,
-	R_2,
-	R_1_I,
-	R_DIFF_ZERO,
-	R_DIFF_IGNORE,
-	R_DIFF_INTENSITY,
-};
-
-
-/* Return the least squares optimal scaling factor when comparing intensities.
- * list1,list2 are the two intensity lists to compare.
- */
-double stat_scale_intensity(RefList *list1, RefList *list2)
-{
-	double top = 0.0;
-	double bot = 0.0;
-	Reflection *refl1;
-	RefListIterator *iter;
-
-	for ( refl1 = first_refl(list1, &iter);
-	      refl1 != NULL;
-	      refl1 = next_refl(refl1, iter) )
-	{
-		double i1, i2;
-		signed int h, k, l;
-		Reflection *refl2;
-
-		get_indices(refl1, &h, &k, &l);
-		refl2 = find_refl(list2, h, k, l);
-		if ( refl2 == NULL ) continue;  /* No common reflection */
-
-		i1 = get_intensity(refl1);
-		i2 = get_intensity(refl2);
-
-		top += i1 * i2;
-		bot += i2 * i2;
-
-	}
-
-	return top/bot;
-}
-
-
-/* Return the least squares optimal scaling factor when comparing the square
- * roots of the intensities (i.e. one approximation to the structure factor
- * moduli).
- * list1,list2 are the two intensity lists to compare (they contain intensities,
- * not square rooted intensities).
- */
-static double stat_scale_sqrti(RefList *list1, RefList *list2)
-{
-	double top = 0.0;
-	double bot = 0.0;
-	Reflection *refl1;
-	RefListIterator *iter;
-
-	for ( refl1 = first_refl(list1, &iter);
-	      refl1 != NULL;
-	      refl1 = next_refl(refl1, iter) )
-	{
-		double i1, i2;
-		double f1, f2;
-		signed int h, k, l;
-		Reflection *refl2;
-
-		get_indices(refl1, &h, &k, &l);
-		refl2 = find_refl(list2, h, k, l);
-		if ( refl2 == NULL ) continue;  /* No common reflection */
-
-		i1 = get_intensity(refl1);
-		i2 = get_intensity(refl2);
-
-		if ( i1 < 0.0 ) continue;
-		f1 = sqrt(i1);
-
-		if ( i2 < 0.0 ) continue;
-		f2 = sqrt(i2);
-
-		top += f1 * f2;
-		bot += f2 * f2;
-
-	}
-
-	return top/bot;
-}
-
-
-static double internal_r1_ignorenegs(RefList *list1, RefList *list2,
-                                     double scale)
-{
-	double top = 0.0;
-	double bot = 0.0;
-	Reflection *refl1;
-	RefListIterator *iter;
-
-	for ( refl1 = first_refl(list1, &iter);
-	      refl1 != NULL;
-	      refl1 = next_refl(refl1, iter) )
-	{
-		double i1, i2;
-		double f1, f2;
-		signed int h, k, l;
-		Reflection *refl2;
-
-		get_indices(refl1, &h, &k, &l);
-		refl2 = find_refl(list2, h, k, l);
-		if ( refl2 == NULL ) continue;  /* No common reflection */
-
-		i1 = get_intensity(refl1);
-		i2 = get_intensity(refl2);
-
-		if ( i1 < 0.0 ) continue;
-		f1 = sqrt(i1);
-
-		if ( i2 < 0.0 ) continue;
-		f2 = sqrt(i2);
-		f2 *= scale;
-
-		top += fabs(f1 - f2);
-		bot += f1;
-
-	}
-
-	return top/bot;
-}
-
-
-static double internal_r1_negstozero(RefList *list1, RefList *list2,
-                                     double scale)
-{
-	double top = 0.0;
-	double bot = 0.0;
-	Reflection *refl1;
-	RefListIterator *iter;
-
-	for ( refl1 = first_refl(list1, &iter);
-	      refl1 != NULL;
-	      refl1 = next_refl(refl1, iter) )
-	{
-		double i1, i2;
-		double f1, f2;
-		signed int h, k, l;
-		Reflection *refl2;
-
-		get_indices(refl1, &h, &k, &l);
-		refl2 = find_refl(list2, h, k, l);
-		if ( refl2 == NULL ) continue;  /* No common reflection */
-
-		i1 = get_intensity(refl1);
-		i2 = get_intensity(refl2);
-
-		f1 = i1 > 0.0 ? sqrt(i1) : 0.0;
-
-		f2 = i2 > 0.0 ? sqrt(i2) : 0.0;
-		f2 *= scale;
-
-		top += fabs(f1 - f2);
-		bot += f1;
-
-	}
-
-	return top/bot;
-}
-
-
-static double internal_r2(RefList *list1, RefList *list2, double scale)
-{
-	double top = 0.0;
-	double bot = 0.0;
-	Reflection *refl1;
-	RefListIterator *iter;
-
-	for ( refl1 = first_refl(list1, &iter);
-	      refl1 != NULL;
-	      refl1 = next_refl(refl1, iter) )
-	{
-		double i1, i2;
-		signed int h, k, l;
-		Reflection *refl2;
-
-		get_indices(refl1, &h, &k, &l);
-		refl2 = find_refl(list2, h, k, l);
-		if ( refl2 == NULL ) continue;  /* No common reflection */
-
-		i1 = get_intensity(refl1);
-		i2 = get_intensity(refl2);
-
-		i2 *= scale;
-
-		top += pow(i1 - i2, 2.0);
-		bot += pow(i1, 2.0);
-
-	}
-
-	return sqrt(top/bot);
-}
-
-
-static double internal_r_i(RefList *list1, RefList *list2, double scale)
-{
-	double top = 0.0;
-	double bot = 0.0;
-	Reflection *refl1;
-	RefListIterator *iter;
-
-	for ( refl1 = first_refl(list1, &iter);
-	      refl1 != NULL;
-	      refl1 = next_refl(refl1, iter) )
-	{
-		double i1, i2;
-		signed int h, k, l;
-		Reflection *refl2;
-
-		get_indices(refl1, &h, &k, &l);
-		refl2 = find_refl(list2, h, k, l);
-		if ( refl2 == NULL ) continue;  /* No common reflection */
-
-		i1 = get_intensity(refl1);
-		i2 = get_intensity(refl2);
-		i2 *= scale;
-
-		top += fabs(i1-i2);
-		bot += fabs(i1);
-
-	}
-
-	return top/bot;
-}
-
-
-static double internal_rdiff_intensity(RefList *list1, RefList *list2,
-                                       double scale)
-{
-	double top = 0.0;
-	double bot = 0.0;
-	Reflection *refl1;
-	RefListIterator *iter;
-
-	for ( refl1 = first_refl(list1, &iter);
-	      refl1 != NULL;
-	      refl1 = next_refl(refl1, iter) )
-	{
-		double i1, i2;
-		signed int h, k, l;
-		Reflection *refl2;
-
-		get_indices(refl1, &h, &k, &l);
-		refl2 = find_refl(list2, h, k, l);
-		if ( refl2 == NULL ) continue;  /* No common reflection */
-
-		i1 = get_intensity(refl1);
-		i2 = get_intensity(refl2);
-		i2 *= scale;
-
-		top += fabs(i1 - i2);
-		bot += i1 + i2;
-
-	}
-
-	return 2.0*top/bot;
-}
-
-
-static double internal_rdiff_negstozero(RefList *list1, RefList *list2,
-                                        double scale)
-{
-	double top = 0.0;
-	double bot = 0.0;
-	Reflection *refl1;
-	RefListIterator *iter;
-
-	for ( refl1 = first_refl(list1, &iter);
-	      refl1 != NULL;
-	      refl1 = next_refl(refl1, iter) )
-	{
-		double i1, i2;
-		double f1, f2;
-		signed int h, k, l;
-		Reflection *refl2;
-
-		get_indices(refl1, &h, &k, &l);
-		refl2 = find_refl(list2, h, k, l);
-		if ( refl2 == NULL ) continue;  /* No common reflection */
-
-		i1 = get_intensity(refl1);
-		i2 = get_intensity(refl2);
-
-		f1 = i1 > 0.0 ? sqrt(i1) : 0.0;
-
-		f2 = i2 > 0.0 ? sqrt(i2) : 0.0;
-		f2 *= scale;
-
-		top += fabs(f1 - f2);
-		bot += f1 + f2;
-
-	}
-
-	return 2.0*top/bot;
-}
-
-
-static double internal_rdiff_ignorenegs(RefList *list1, RefList *list2,
-                                        double scale)
-{
-	double top = 0.0;
-	double bot = 0.0;
-	Reflection *refl1;
-	RefListIterator *iter;
-
-	for ( refl1 = first_refl(list1, &iter);
-	      refl1 != NULL;
-	      refl1 = next_refl(refl1, iter) )
-	{
-		double i1, i2;
-		double f1, f2;
-		signed int h, k, l;
-		Reflection *refl2;
-
-		get_indices(refl1, &h, &k, &l);
-		refl2 = find_refl(list2, h, k, l);
-		if ( refl2 == NULL ) continue;  /* No common reflection */
-
-		i1 = get_intensity(refl1);
-		i2 = get_intensity(refl2);
-
-		if ( i1 < 0.0 ) continue;
-		f1 = sqrt(i1);
-
-		if ( i2 < 0.0 ) continue;
-		f2 = sqrt(i2);
-		f2 *= scale;
-
-		top += fabs(f1 - f2);
-		bot += f1 + f2;
-
-	}
-
-	return 2.0*top/bot;
-}
-
-
-static double calc_r(double scale, void *params)
-{
-	struct r_params *rp = params;
-
-	switch ( rp->fom ) {
-	case R_1_ZERO :
-		return internal_r1_negstozero(rp->list1, rp->list2, scale);
-	case R_1_IGNORE :
-		return internal_r1_ignorenegs(rp->list1, rp->list2, scale);
-	case R_2 :
-		return internal_r2(rp->list1, rp->list2, scale);
-
-	case R_1_I :
-		return internal_r_i(rp->list1, rp->list2, scale);
-
-	case R_DIFF_ZERO :
-		return internal_rdiff_negstozero(rp->list1, rp->list2,scale);
-	case R_DIFF_IGNORE :
-		return internal_rdiff_ignorenegs(rp->list1, rp->list2, scale);
-	case R_DIFF_INTENSITY :
-		return internal_rdiff_intensity(rp->list1, rp->list2, scale);
-	}
-
-	ERROR("No such FoM!\n");
-	abort();
-}
-
-
-static double r_minimised(RefList *list1, RefList *list2, double *scalep, int fom,
-                          int u)
-{
-	gsl_function F;
-	gsl_min_fminimizer *s;
-	int status;
-	double scale = 1.0;
-	struct r_params rp;
-	int iter = 0;
-
-	rp.list1 = list1;
-	rp.list2 = list2;
-	rp.fom = fom;
-
-	if ( u ) {
-
-		scale = 1.0;
-
-	} else {
-
-		F.function = &calc_r;
-		F.params = &rp;
-
-		s = gsl_min_fminimizer_alloc(gsl_min_fminimizer_brent);
-
-		/* Initial guess */
-		switch ( fom ) {
-		case R_1_ZERO :
-		case R_1_IGNORE :
-		case R_DIFF_ZERO :
-		case R_DIFF_IGNORE :
-			scale = stat_scale_sqrti(list1, list2);
-			break;
-		case R_2 :
-		case R_1_I :
-		case R_DIFF_INTENSITY :
-			scale = stat_scale_intensity(list1, list2);
-			break;
-		}
-		//STATUS("Initial scale factor estimate: %5.2e\n", scale);
-
-		/* Probably within an order of magnitude either side */
-		gsl_min_fminimizer_set(s, &F, scale, scale/10.0, scale*10.0);
-
-		do {
-
-			double lo, up;
-
-			/* Iterate */
-			if ( gsl_min_fminimizer_iterate(s) ) {
-				ERROR("Failed to find scale factor.\n");
-				return NAN;
-			}
-
-			/* Get the current estimate */
-			scale = gsl_min_fminimizer_x_minimum(s);
-			lo = gsl_min_fminimizer_x_lower(s);
-			up = gsl_min_fminimizer_x_upper(s);
-
-			/* Check for convergence */
-			status = gsl_min_test_interval(lo, up, 0.001, 0.0);
-
-			iter++;
-
-		} while ( status == GSL_CONTINUE );
-
-		if ( status != GSL_SUCCESS ) {
-			ERROR("Scale factor minimisation failed.\n");
-		}
-
-		gsl_min_fminimizer_free(s);
-
-	}
-
-	//STATUS("Final scale factor: %5.2e\n", scale);
-	*scalep = scale;
-	return calc_r(scale, &rp);
-}
-
-
-double stat_r1_ignore(RefList *list1, RefList *list2, double *scalep, int u)
-{
-	return r_minimised(list1, list2, scalep, R_1_IGNORE, u);
-}
-
-
-double stat_r1_zero(RefList *list1, RefList *list2, double *scalep, int u)
-{
-	return r_minimised(list1, list2, scalep, R_1_ZERO, u);
-}
-
-
-double stat_r2(RefList *list1, RefList *list2, double *scalep, int u)
-{
-	return r_minimised(list1, list2, scalep, R_2, u);
-}
-
-
-double stat_r1_i(RefList *list1, RefList *list2, double *scalep, int u)
-{
-	return r_minimised(list1, list2, scalep, R_1_I, u);
-}
-
-
-double stat_rdiff_zero(RefList *list1, RefList *list2, double *scalep, int u)
-{
-	return r_minimised(list1, list2, scalep, R_DIFF_ZERO, u);
-}
-
-
-double stat_rdiff_ignore(RefList *list1, RefList *list2, double *scalep, int u)
-{
-	return r_minimised(list1, list2, scalep, R_DIFF_IGNORE, u);
-}
-
-
-double stat_rdiff_intensity(RefList *list1, RefList *list2, double *scalep, int u)
-{
-	return r_minimised(list1, list2,  scalep, R_DIFF_INTENSITY, u);
-}
-
-
-double stat_pearson_i(RefList *list1, RefList *list2)
-{
-	double *vec1, *vec2;
-	int ni = num_reflections(list1);
-	double val;
-	int nacc = 0;
-	Reflection *refl1;
-	RefListIterator *iter;
-
-	vec1 = malloc(ni*sizeof(double));
-	vec2 = malloc(ni*sizeof(double));
-
-	for ( refl1 = first_refl(list1, &iter);
-	      refl1 != NULL;
-	      refl1 = next_refl(refl1, iter) )
-	{
-		double i1, i2;
-		signed int h, k, l;
-		Reflection *refl2;
-
-		get_indices(refl1, &h, &k, &l);
-		refl2 = find_refl(list2, h, k, l);
-		if ( refl2 == NULL ) continue;  /* No common reflection */
-
-		i1 = get_intensity(refl1);
-		i2 = get_intensity(refl2);
-
-		vec1[nacc] = i1;
-		vec2[nacc] = i2;
-		nacc++;
-	}
-
-	val = gsl_stats_correlation(vec1, 1, vec2, 1, nacc);
-	free(vec1);
-	free(vec2);
-
-	return val;
-}
-
-
-double stat_pearson_f_ignore(RefList *list1, RefList *list2)
-{
-	double *vec1, *vec2;
-	int ni = num_reflections(list1);
-	double val;
-	int nacc = 0;
-	Reflection *refl1;
-	RefListIterator *iter;
-
-	vec1 = malloc(ni*sizeof(double));
-	vec2 = malloc(ni*sizeof(double));
-
-	for ( refl1 = first_refl(list1, &iter);
-	      refl1 != NULL;
-	      refl1 = next_refl(refl1, iter) )
-	{
-		double i1, i2;
-		double f1, f2;
-		signed int h, k, l;
-		Reflection *refl2;
-
-		get_indices(refl1, &h, &k, &l);
-		refl2 = find_refl(list2, h, k, l);
-		if ( refl2 == NULL ) continue;  /* No common reflection */
-
-		i1 = get_intensity(refl1);
-		i2 = get_intensity(refl2);
-
-		if ( i1 < 0.0 ) continue;
-		if ( i2 < 0.0 ) continue;
-
-		f1 = sqrt(i1);
-		f2 = sqrt(i2);
-
-		vec1[nacc] = f1;
-		vec2[nacc] = f2;
-		nacc++;
-
-	}
-
-	val = gsl_stats_correlation(vec1, 1, vec2, 1, nacc);
-	free(vec1);
-	free(vec2);
-
-	return val;
-}
-
-
-double stat_pearson_f_zero(RefList *list1, RefList *list2)
-{
-	double *vec1, *vec2;
-	int ni = num_reflections(list1);
-	double val;
-	int nacc = 0;
-	Reflection *refl1;
-	RefListIterator *iter;
-
-	vec1 = malloc(ni*sizeof(double));
-	vec2 = malloc(ni*sizeof(double));
-
-	for ( refl1 = first_refl(list1, &iter);
-	      refl1 != NULL;
-	      refl1 = next_refl(refl1, iter) )
-	{
-		double i1, i2;
-		double f1, f2;
-		signed int h, k, l;
-		Reflection *refl2;
-
-		get_indices(refl1, &h, &k, &l);
-		refl2 = find_refl(list2, h, k, l);
-		if ( refl2 == NULL ) continue;  /* No common reflection */
-
-		i1 = get_intensity(refl1);
-		i2 = get_intensity(refl2);
-
-		f1 = i1 > 0.0 ? sqrt(i1) : 0.0;
-		f2 = i2 > 0.0 ? sqrt(i2) : 0.0;
-
-		vec1[nacc] = f1;
-		vec2[nacc] = f2;
-		nacc++;
-
-	}
-
-	val = gsl_stats_correlation(vec1, 1, vec2, 1, nacc);
-	free(vec1);
-	free(vec2);
-
-	return val;
-}
diff --git a/src/statistics.h b/src/statistics.h
deleted file mode 100644
index 8fa78ea6..00000000
--- a/src/statistics.h
+++ /dev/null
@@ -1,45 +0,0 @@
-/*
- * statistics.h
- *
- * Structure-factor statistics
- *
- * (c) 2006-2010 Thomas White <taw@physics.org>
- *
- * Part of CrystFEL - crystallography with a FEL
- *
- */
-
-#ifdef HAVE_CONFIG_H
-#include <config.h>
-#endif
-
-#ifndef STATISTICS_H
-#define STATISTICS_H
-
-
-#include "reflist.h"
-
-extern double stat_scale_intensity(RefList *list1, RefList *list2);
-
-extern double stat_r1_zero(RefList *list1, RefList *list2,
-                           double *scalep, int u);
-extern double stat_r1_ignore(RefList *list1, RefList *list2,
-                             double *scalep, int u);
-
-extern double stat_r2(RefList *list1, RefList *list2, double *scalep, int u);
-
-extern double stat_r1_i(RefList *list1, RefList *list2, double *scalep, int u);
-
-extern double stat_rdiff_zero(RefList *list1, RefList *list2,
-                              double *scalep, int u);
-extern double stat_rdiff_ignore(RefList *list1, RefList *list2,
-                                double *scalep, int u);
-extern double stat_rdiff_intensity(RefList *list1, RefList *list2,
-                                   double *scalep, int u);
-
-extern double stat_pearson_i(RefList *list1, RefList *list2);
-extern double stat_pearson_f_zero(RefList *list1, RefList *list2);
-extern double stat_pearson_f_ignore(RefList *list1, RefList *list2);
-
-
-#endif	/* STATISTICS_H */
diff --git a/src/stream.c b/src/stream.c
deleted file mode 100644
index a7cdc2d9..00000000
--- a/src/stream.c
+++ /dev/null
@@ -1,487 +0,0 @@
-/*
- * stream.c
- *
- * Stream tools
- *
- * (c) 2006-2011 Thomas White <taw@physics.org>
- * (c) 2011 Rick Kirian <rkirian@asu.edu>
- *
- * Part of CrystFEL - crystallography with a FEL
- *
- */
-
-#ifdef HAVE_CONFIG_H
-#include <config.h>
-#endif
-
-
-#include <stdlib.h>
-#include <stdio.h>
-#include <string.h>
-
-#include "cell.h"
-#include "utils.h"
-#include "image.h"
-#include "stream.h"
-#include "reflist.h"
-#include "reflist-utils.h"
-
-
-#define CHUNK_START_MARKER "----- Begin chunk -----"
-#define CHUNK_END_MARKER "----- End chunk -----"
-#define PEAK_LIST_START_MARKER "Peaks from peak search"
-#define PEAK_LIST_END_MARKER "End of peak list"
-#define REFLECTION_START_MARKER "Reflections measured after indexing"
-/* REFLECTION_END_MARKER is over in reflist-utils.h because it is also
- * used to terminate a standalone list of reflections */
-
-static void exclusive(const char *a, const char *b)
-{
-	ERROR("The stream options '%s' and '%s' are mutually exclusive.\n",
-	      a, b);
-}
-
-
-int parse_stream_flags(const char *a)
-{
-	int n, i;
-	int ret = STREAM_NONE;
-	char **flags;
-
-	n = assplode(a, ",", &flags, ASSPLODE_NONE);
-
-	for ( i=0; i<n; i++ ) {
-
-		if ( strcmp(flags[i], "integrated") == 0) {
-
-			ret |= STREAM_INTEGRATED;
-
-		} else if ( strcmp(flags[i], "peaks") == 0) {
-			if ( ret & STREAM_PEAKS_IF_INDEXED ) {
-				exclusive("peaks", "peaksifindexed");
-				return -1;
-			}
-			if ( ret & STREAM_PEAKS_IF_NOT_INDEXED ) {
-				exclusive("peaks", "peaksifnotindexed");
-				return -1;
-			}
-			ret |= STREAM_PEAKS;
-
-		} else if ( strcmp(flags[i], "peaksifindexed") == 0) {
-			if ( ret & STREAM_PEAKS ) {
-				exclusive("peaks", "peaksifindexed");
-				return -1;
-			}
-			if ( ret & STREAM_PEAKS_IF_NOT_INDEXED ) {
-				exclusive("peaksifnotindexed",
-				          "peaksifindexed");
-				return -1;
-			}
-			ret |= STREAM_PEAKS_IF_INDEXED;
-
-		} else if ( strcmp(flags[i], "peaksifnotindexed") == 0) {
-			if ( ret & STREAM_PEAKS ) {
-				exclusive("peaks", "peaksifnotindexed");
-				return -1;
-			}
-			if ( ret & STREAM_PEAKS_IF_INDEXED ) {
-				exclusive("peaksifnotindexed",
-				          "peaksifindexed");
-				return -1;
-			}
-			ret |= STREAM_PEAKS_IF_NOT_INDEXED;
-
-		} else {
-			ERROR("Unrecognised stream flag '%s'\n", flags[i]);
-			return -1;
-		}
-
-		free(flags[i]);
-
-	}
-	free(flags);
-
-	return ret;
-}
-
-
-int count_patterns(FILE *fh)
-{
-	char *rval;
-
-	int n_total_patterns = 0;
-	do {
-		char line[1024];
-
-		rval = fgets(line, 1023, fh);
-		if ( rval == NULL ) continue;
-		chomp(line);
-		if ( strcmp(line, CHUNK_END_MARKER) == 0 ) n_total_patterns++;
-
-	} while ( rval != NULL );
-
-	if ( ferror(fh) ) {
-		ERROR("Read error while counting patterns.\n");
-		return 0;
-	}
-
-	return n_total_patterns;
-}
-
-
-static int read_peaks(FILE *fh, struct image *image)
-{
-	char *rval = NULL;
-	int first = 1;
-
-	image->features = image_feature_list_new();
-
-	do {
-
-		char line[1024];
-		float x, y, d, intensity;
-		int r;
-
-		rval = fgets(line, 1023, fh);
-		if ( rval == NULL ) continue;
-		chomp(line);
-
-		if ( strcmp(line, PEAK_LIST_END_MARKER) == 0 ) return 0;
-
-		r = sscanf(line, "%f %f %f %f", &x, &y, &d, &intensity);
-		if ( (r != 4) && (!first) ) {
-			ERROR("Failed to parse peak list line.\n");
-			ERROR("The failed line was: '%s'\n", line);
-			return 1;
-		}
-
-		first = 0;
-		if ( r == 4 ) {
-			image_add_feature(image->features, x, y,
-			                  image, intensity, NULL);
-		}
-
-	} while ( rval != NULL );
-
-	/* Got read error of some kind before finding PEAK_LIST_END_MARKER */
-	return 1;
-}
-
-
-static void write_peaks(struct image *image, FILE *ofh)
-{
-	int i;
-
-	fprintf(ofh, PEAK_LIST_START_MARKER"\n");
-	fprintf(ofh, "  fs/px   ss/px  (1/d)/nm^-1   Intensity\n");
-
-	for ( i=0; i<image_feature_count(image->features); i++ ) {
-
-		struct imagefeature *f;
-		struct rvec r;
-		double q;
-
-		f = image_get_feature(image->features, i);
-		if ( f == NULL ) continue;
-
-		r = get_q(image, f->fs, f->ss, NULL, 1.0/image->lambda);
-		q = modulus(r.u, r.v, r.w);
-
-		fprintf(ofh, "%7.2f %7.2f   %10.2f  %10.2f\n",
-		       f->fs, f->ss, q/1.0e9, f->intensity);
-
-	}
-
-	fprintf(ofh, PEAK_LIST_END_MARKER"\n");
-}
-
-
-void write_chunk(FILE *ofh, struct image *i, struct hdfile *hdfile, int f)
-{
-	double asx, asy, asz;
-	double bsx, bsy, bsz;
-	double csx, csy, csz;
-	double a, b, c, al, be, ga;
-
-	fprintf(ofh, CHUNK_START_MARKER"\n");
-
-	fprintf(ofh, "Image filename: %s\n", i->filename);
-
-	if ( i->indexed_cell != NULL ) {
-
-		cell_get_parameters(i->indexed_cell, &a, &b, &c,
-		                                         &al, &be, &ga);
-		fprintf(ofh, "Cell parameters %7.5f %7.5f %7.5f nm,"
-			     " %7.5f %7.5f %7.5f deg\n",
-			     a*1.0e9, b*1.0e9, c*1.0e9,
-			     rad2deg(al), rad2deg(be), rad2deg(ga));
-
-		cell_get_reciprocal(i->indexed_cell, &asx, &asy, &asz,
-			                  &bsx, &bsy, &bsz,
-			                  &csx, &csy, &csz);
-		fprintf(ofh, "astar = %+9.7f %+9.7f %+9.7f nm^-1\n",
-			asx/1e9, asy/1e9, asz/1e9);
-		fprintf(ofh, "bstar = %+9.7f %+9.7f %+9.7f nm^-1\n",
-			bsx/1e9, bsy/1e9, bsz/1e9);
-		fprintf(ofh, "cstar = %+9.7f %+9.7f %+9.7f nm^-1\n",
-		        csx/1e9, csy/1e9, csz/1e9);
-
-	} else {
-
-		fprintf(ofh, "No unit cell from indexing.\n");
-
-	}
-
-	if ( i->i0_available ) {
-		fprintf(ofh, "I0 = %7.5f (arbitrary units)\n", i->i0);
-	} else {
-		fprintf(ofh, "I0 = invalid\n");
-	}
-
-	fprintf(ofh, "photon_energy_eV = %f\n",
-	        J_to_eV(ph_lambda_to_en(i->lambda)));
-
-	if ( i->det != NULL ) {
-
-		int j;
-
-		for ( j=0; j<i->det->n_panels; j++ ) {
-			fprintf(ofh, "camera_length_%s = %f\n",
-			        i->det->panels[j].name, i->det->panels[j].clen);
-		}
-
-	}
-
-	copy_hdf5_fields(hdfile, i->copyme, ofh);
-
-	if ( (f & STREAM_PEAKS)
-	  || ((f & STREAM_PEAKS_IF_INDEXED) && (i->indexed_cell != NULL))
-	  || ((f & STREAM_PEAKS_IF_NOT_INDEXED) && (i->indexed_cell == NULL)) )
-	{
-		fprintf(ofh, "\n");
-		write_peaks(i, ofh);
-	}
-
-	if ( f & STREAM_INTEGRATED ) {
-
-		fprintf(ofh, "\n");
-
-		if ( i->reflections != NULL ) {
-
-			fprintf(ofh, REFLECTION_START_MARKER"\n");
-			write_reflections_to_file(ofh, i->reflections,
-			                               i->indexed_cell);
-			fprintf(ofh, REFLECTION_END_MARKER"\n");
-
-		} else {
-
-			fprintf(ofh, "No integrated reflections.\n");
-
-		}
-	}
-
-	fprintf(ofh, CHUNK_END_MARKER"\n\n");
-}
-
-
-static int find_start_of_chunk(FILE *fh)
-{
-	char *rval = NULL;
-	char line[1024];
-
-	do {
-
-		rval = fgets(line, 1023, fh);
-
-		/* Trouble? */
-		if ( rval == NULL ) return 1;
-
-		chomp(line);
-
-	} while ( strcmp(line, CHUNK_START_MARKER) != 0 );
-
-	return 0;
-}
-
-
-/* Read the next chunk from a stream and fill in 'image' */
-int read_chunk(FILE *fh, struct image *image)
-{
-	char line[1024];
-	char *rval = NULL;
-	struct rvec as, bs, cs;
-	int have_as = 0;
-	int have_bs = 0;
-	int have_cs = 0;
-	int have_filename = 0;
-	int have_cell = 0;
-	int have_ev = 0;
-
-	if ( find_start_of_chunk(fh) ) return 1;
-
-	image->i0_available = 0;
-	image->i0 = 1.0;
-	image->lambda = -1.0;
-	image->features = NULL;
-	image->reflections = NULL;
-	image->indexed_cell = NULL;
-
-	do {
-
-		float u, v, w;
-
-		rval = fgets(line, 1023, fh);
-
-		/* Trouble? */
-		if ( rval == NULL ) break;
-
-		chomp(line);
-
-		if ( strncmp(line, "Image filename: ", 16) == 0 ) {
-			image->filename = strdup(line+16);
-			have_filename = 1;
-		}
-
-		if ( strncmp(line, "camera_length_", 14) == 0 ) {
-			if ( image->det != NULL ) {
-
-				int k;
-				char name[1024];
-				struct panel *p;
-
-				for ( k=0; k<strlen(line)-14; k++ ) {
-					char ch = line[k+14];
-					name[k] = ch;
-					if ( (ch == ' ') || (ch == '=') ) {
-						name[k] = '\0';
-						break;
-					}
-				}
-
-				p = find_panel_by_name(image->det, name);
-				if ( p == NULL ) {
-					ERROR("No panel '%s'\n", name);
-				} else {
-					p->clen = atof(line+14+k+3);
-				}
-
-			}
-		}
-
-		if ( strncmp(line, "I0 = ", 5) == 0 ) {
-			image->i0 = atof(line+5);
-			image->i0_available = 1;
-		}
-
-		if ( sscanf(line, "astar = %f %f %f", &u, &v, &w) == 3 ) {
-			as.u = u*1e9;  as.v = v*1e9;  as.w = w*1e9;
-			have_as = 1;
-		}
-
-		if ( sscanf(line, "bstar = %f %f %f", &u, &v, &w) == 3 ) {
-			bs.u = u*1e9;  bs.v = v*1e9;  bs.w = w*1e9;
-			have_bs = 1;
-		}
-
-		if ( sscanf(line, "cstar = %f %f %f", &u, &v, &w) == 3 ) {
-			cs.u = u*1e9;  cs.v = v*1e9;  cs.w = w*1e9;
-			have_cs = 1;
-		}
-
-		if ( have_as && have_bs && have_cs ) {
-			if ( image->indexed_cell != NULL ) {
-				ERROR("Duplicate cell found in stream!\n");
-				cell_free(image->indexed_cell);
-			}
-			image->indexed_cell = cell_new_from_reciprocal_axes(as,
-			                                                    bs,
-			                                                    cs);
-			have_cell = 1;
-			have_as = 0;  have_bs = 0;  have_cs = 0;
-		}
-
-		if ( strncmp(line, "photon_energy_eV = ", 19) == 0 ) {
-			image->lambda = ph_en_to_lambda(eV_to_J(atof(line+19)));
-			have_ev = 1;
-		}
-
-		if ( strcmp(line, PEAK_LIST_START_MARKER) == 0 ) {
-			if ( read_peaks(fh, image) ) {
-				ERROR("Failed while reading peaks\n");
-				return 1;
-			}
-		}
-
-		if ( strcmp(line, REFLECTION_START_MARKER) == 0 ) {
-			image->reflections = read_reflections_from_file(fh);
-			if ( image->reflections == NULL ) {
-				ERROR("Failed while reading reflections\n");
-				return 1;
-			}
-		}
-
-		if ( strcmp(line, CHUNK_END_MARKER) == 0 ) {
-			if ( have_filename && have_ev ) return 0;
-			ERROR("Incomplete chunk found in input file.\n");
-			return 1;
-		}
-
-	} while ( 1 );
-
-	return 1;  /* Either error or EOF, don't care because we will complain
-	            * on the terminal if it was an error. */
-}
-
-
-void write_stream_header(FILE *ofh, int argc, char *argv[])
-{
-	int i;
-
-	fprintf(ofh, "CrystFEL stream format 2.0\n");
-	fprintf(ofh, "Command line:");
-	for ( i=0; i<argc; i++ ) {
-		fprintf(ofh, " %s", argv[i]);
-	}
-	fprintf(ofh, "\n");
-	fflush(ofh);
-}
-
-
-int skip_some_files(FILE *fh, int n)
-{
-	char *rval = NULL;
-	int n_patterns = 0;
-
-	do {
-
-		char line[1024];
-
-		if ( n_patterns == n ) return 0;
-
-		rval = fgets(line, 1023, fh);
-		if ( rval == NULL ) continue;
-		if ( strcmp(line, CHUNK_END_MARKER) == 0 ) n_patterns++;
-
-	} while ( rval != NULL );
-
-	return 1;
-}
-
-int is_stream(const char *filename) {
-	FILE *fh;
-	char line[1024];
-	char *rval = NULL;
-	fh = fopen(filename, "r");
-	rval = fgets(line, 1023, fh);
-	fclose(fh);
-	if ( rval == NULL ) {
-		return -1;
-	}
-	if ( strncmp(line, "CrystFEL stream format 2.0", 26) == 0 ) {
-		return 1;
-	}
-	else {
-		return 0;
-	}
-	return -1;
-}
diff --git a/src/stream.h b/src/stream.h
deleted file mode 100644
index ba218fb9..00000000
--- a/src/stream.h
+++ /dev/null
@@ -1,49 +0,0 @@
-/*
- * stream.h
- *
- * Stream tools
- *
- * (c) 2006-2011 Thomas White <taw@physics.org>
- *
- * Part of CrystFEL - crystallography with a FEL
- *
- */
-
-#ifndef STREAM_H
-#define STREAM_H
-
-#ifdef HAVE_CONFIG_H
-#include <config.h>
-#endif
-
-
-struct image;
-struct hdfile;
-
-/* Possible options dictating what goes into the output stream */
-enum
-{
-	STREAM_NONE                 = 0,
-	STREAM_INTEGRATED           = 1<<0,
-	STREAM_PEAKS                = 1<<2,
-	STREAM_PEAKS_IF_INDEXED     = 1<<3,
-	STREAM_PEAKS_IF_NOT_INDEXED = 1<<4,
-};
-
-
-extern int count_patterns(FILE *fh);
-
-extern void write_stream_header(FILE *ofh, int argc, char *argv[]);
-
-extern void write_chunk(FILE *ofh, struct image *image, struct hdfile *hdfile,
-                        int flags);
-
-extern int parse_stream_flags(const char *a);
-
-extern int read_chunk(FILE *fh, struct image *image);
-
-extern int skip_some_files(FILE *fh, int n);
-
-extern int is_stream(const char *filename);
-
-#endif	/* STREAM_H */
diff --git a/src/symmetry.c b/src/symmetry.c
deleted file mode 100644
index f0b24146..00000000
--- a/src/symmetry.c
+++ /dev/null
@@ -1,1503 +0,0 @@
-/*
- * symmetry.c
- *
- * Symmetry
- *
- * (c) 2006-2010 Thomas White <taw@physics.org>
- *
- * Part of CrystFEL - crystallography with a FEL
- *
- */
-
-
-#ifdef HAVE_CONFIG_H
-#include <config.h>
-#endif
-
-#include <stdlib.h>
-#include <stdio.h>
-#include <string.h>
-#include <math.h>
-#include <assert.h>
-
-#include "symmetry.h"
-#include "utils.h"
-
-
-/**
- * SECTION:symmetry
- * @short_description: Point symmetry handling
- * @title: Symmetry
- * @section_id:
- * @see_also:
- * @include: "symmetry.h"
- * @Image:
- *
- * Routines to handle point symmetry.
- */
-
-
-struct sym_op
-{
-	signed int *h;
-	signed int *k;
-	signed int *l;  /* Contributions to h, k and l from h, k, i and l */
-	int order;
-};
-
-
-struct _symoplist
-{
-	struct sym_op *ops;
-	int n_ops;
-	int max_ops;
-	char *name;
-	int *divisors;
-	int num_equivs;
-};
-
-
-struct _symopmask
-{
-	const SymOpList *list;
-	int *mask;
-};
-
-
-
-static void alloc_ops(SymOpList *ops)
-{
-	ops->ops = realloc(ops->ops, ops->max_ops*sizeof(struct sym_op));
-	ops->divisors = realloc(ops->divisors, ops->max_ops*sizeof(int));
-}
-
-
-/**
- * new_symopmask:
- * @list: A %SymOpList
- *
- * Returns: a new %SymOpMask, which you can use when filtering out special
- * reflections.
- **/
-SymOpMask *new_symopmask(const SymOpList *list)
-{
-	SymOpMask *m;
-	int i;
-
-	m = malloc(sizeof(struct _symopmask));
-	if ( m == NULL ) return NULL;
-
-	m->list = list;
-	m->mask = malloc(sizeof(int)*list->n_ops);
-	if ( m->mask == NULL ) {
-		free(m);
-		return NULL;
-	}
-
-	for ( i=0; i<list->n_ops; i++ ) {
-		m->mask[i] = 1;
-	}
-
-	return m;
-}
-
-
-/* Creates a new SymOpList */
-static SymOpList *new_symoplist()
-{
-	SymOpList *new;
-	new = malloc(sizeof(SymOpList));
-	if ( new == NULL ) return NULL;
-	new->max_ops = 16;
-	new->n_ops = 0;
-	new->ops = NULL;
-	new->divisors = NULL;
-	new->name = NULL;
-	new->num_equivs = 1;
-	alloc_ops(new);
-	return new;
-}
-
-
-/**
- * free_symoplist:
- * @ops: A %SymOpList to free
- *
- * Frees a %SymOpList and all associated resources.
- **/
-void free_symoplist(SymOpList *ops)
-{
-	int i;
-
-	if ( ops == NULL ) return;
-	for ( i=0; i<ops->n_ops; i++ ) {
-		free(ops->ops[i].h);
-		free(ops->ops[i].k);
-		free(ops->ops[i].l);
-	}
-	if ( ops->ops != NULL ) free(ops->ops);
-	if ( ops->name != NULL ) free(ops->name);
-	free(ops);
-}
-
-/**
- * free_symopmask:
- * @m: A %SymOpMask to free
- *
- * Frees a %SymOpMask and all associated resources.
- **/
-void free_symopmask(SymOpMask *m)
-{
-	if ( m == NULL ) return;
-	free(m->mask);
-	free(m);
-}
-
-
-/* This returns the number of operations in "ops".  This might be different
- * to num_equivs() if the point group is being constructed. */
-static int num_ops(const SymOpList *ops)
-{
-	return ops->n_ops;
-}
-
-
-/* Add a operation to a SymOpList */
-static void add_symop(SymOpList *ops,
-                      signed int *h, signed int *k, signed int *l,
-                      int order)
-{
-	int n;
-
-	if ( ops->n_ops == ops->max_ops ) {
-		/* Pretty sure this never happens, but still... */
-		ops->max_ops += 16;
-		alloc_ops(ops);
-	}
-
-	n = ops->n_ops;
-	ops->ops[n].h = h;
-	ops->ops[n].k = k;
-	ops->ops[n].l = l;
-	ops->ops[n].order = order;
-	ops->n_ops++;
-}
-
-
-/* Add a operation to a SymOpList */
-static void add_copied_op(SymOpList *ops, struct sym_op *copyme)
-{
-	int n;
-	signed int *h, *k, *l;
-
-	if ( ops->n_ops == ops->max_ops ) {
-		ops->max_ops += 16;
-		alloc_ops(ops);
-	}
-
-	n = ops->n_ops;
-
-	h = malloc(3*sizeof(signed int));
-	k = malloc(3*sizeof(signed int));
-	l = malloc(3*sizeof(signed int));
-
-	memcpy(h, copyme->h, 3*sizeof(signed int));
-	memcpy(k, copyme->k, 3*sizeof(signed int));
-	memcpy(l, copyme->l, 3*sizeof(signed int));
-
-	ops->ops[n].h = h;
-	ops->ops[n].k = k;
-	ops->ops[n].l = l;
-	ops->ops[n].order = copyme->order;
-
-	ops->n_ops++;
-}
-
-
-/**
- * num_equivs:
- * @ops: A %SymOpList
- * @m: A %SymOpMask, which has been shown to special_position()
- *
- * Returns: the number of equivalent reflections for a general reflection
- * in point group "ops", which were not flagged by your call to
- * special_position().
- **/
-int num_equivs(const SymOpList *ops, const SymOpMask *m)
-{
-	int n = num_ops(ops);
-	int i;
-	int c;
-
-	if ( m == NULL ) return n;
-
-	c = 0;
-	for ( i=0; i<n; i++ ) {
-		if ( m->mask[i] ) c++;
-	}
-
-	return c;
-}
-
-
-static signed int *v(signed int h, signed int k, signed int i, signed int l)
-{
-	signed int *vec = malloc(3*sizeof(signed int));
-	/* Convert back to 3-index form now */
-	vec[0] = h-i;  vec[1] = k-i;  vec[2] = l;
-	return vec;
-}
-
-
-static void combine_ops(signed int *h1, signed int *k1, signed int *l1,
-                        signed int *h2, signed int *k2, signed int *l2,
-                        signed int *hnew, signed int *knew, signed int *lnew)
-{
-	/* Yay matrices */
-	hnew[0] = h1[0]*h2[0] + h1[1]*k2[0] + h1[2]*l2[0];
-	hnew[1] = h1[0]*h2[1] + h1[1]*k2[1] + h1[2]*l2[1];
-	hnew[2] = h1[0]*h2[2] + h1[1]*k2[2] + h1[2]*l2[2];
-
-	knew[0] = k1[0]*h2[0] + k1[1]*k2[0] + k1[2]*l2[0];
-	knew[1] = k1[0]*h2[1] + k1[1]*k2[1] + k1[2]*l2[1];
-	knew[2] = k1[0]*h2[2] + k1[1]*k2[2] + k1[2]*l2[2];
-
-	lnew[0] = l1[0]*h2[0] + l1[1]*k2[0] + l1[2]*l2[0];
-	lnew[1] = l1[0]*h2[1] + l1[1]*k2[1] + l1[2]*l2[1];
-	lnew[2] = l1[0]*h2[2] + l1[1]*k2[2] + l1[2]*l2[2];
-}
-
-
-static void combine_and_add_symop(struct sym_op *opi, int oi,
-                                  struct sym_op *opj,
-                                  SymOpList *s)
-{
-	int i;
-	signed int *h, *k, *l;
-
-	h = malloc(3*sizeof(signed int));
-	k = malloc(3*sizeof(signed int));
-	l = malloc(3*sizeof(signed int));
-	assert(h != NULL);
-	assert(k != NULL);
-	assert(l != NULL);
-
-	memcpy(h, opj->h, 3*sizeof(signed int));
-	memcpy(k, opj->k, 3*sizeof(signed int));
-	memcpy(l, opj->l, 3*sizeof(signed int));
-
-	for ( i=0; i<oi; i++ ) {
-
-		signed int hfs[3], kfs[3], lfs[3];
-
-		combine_ops(h, k, l, opi->h, opi->k, opi->l, hfs, kfs, lfs);
-
-		memcpy(h, hfs, 3*sizeof(signed int));
-		memcpy(k, kfs, 3*sizeof(signed int));
-		memcpy(l, lfs, 3*sizeof(signed int));
-
-	}
-
-//	STATUS("Creating %3i %3i %3i\n", h[0], h[1], h[2]);
-//	STATUS("         %3i %3i %3i\n", k[0], k[1], k[2]);
-//	STATUS("         %3i %3i %3i\n", l[0], l[1], l[2]);
-
-	add_symop(s, h, k, l, 1);
-}
-
-
-/* Fill in the other operations for a point group starting from its
- * generators */
-static SymOpList *expand_ops(SymOpList *s)
-{
-	int n, i;
-	SymOpList *e;
-
-	e = new_symoplist();
-	if ( e == NULL ) return NULL;
-	e->name = strdup(symmetry_name(s));
-
-	add_symop(e, v(1,0,0,0), v(0,1,0,0), v(0,0,0,1), 1);  /* I */
-
-	n = num_ops(s);
-	for ( i=0; i<n; i++ ) {
-
-		int j, nj;
-		struct sym_op *opi = &s->ops[i];
-
-		/* Apply op 'i' to all the current ops in the list */
-		nj = num_ops(e);
-		for ( j=0; j<nj; j++ ) {
-
-			int oi;
-
-			for ( oi=0; oi<opi->order-1; oi++ ) {
-				combine_and_add_symop(opi, oi+1, &e->ops[j], e);
-			}
-
-		}
-
-	}
-
-	free_symoplist(s);
-
-	return e;
-}
-
-
-/********************************* Triclinic **********************************/
-
-static SymOpList *make_1bar()
-{
-	SymOpList *new = new_symoplist();
-	add_symop(new, v(-1,0,0,0), v(0,-1,0,0), v(0,0,0,-1), 2);  /* -I */
-	new->name = strdup("-1");
-	return expand_ops(new);
-}
-
-
-static SymOpList *make_1()
-{
-	SymOpList *new = new_symoplist();
-	new->name = strdup("1");
-	return expand_ops(new);
-}
-
-
-/********************************* Monoclinic *********************************/
-
-static SymOpList *make_2m()
-{
-	SymOpList *new = new_symoplist();
-	add_symop(new, v(-1,0,0,0), v(0,-1,0,0), v(0,0,0,1), 2); /* 2 // l */
-	add_symop(new, v(1,0,0,0), v(0,1,0,0), v(0,0,0,-1), 2);  /* m -| l */
-	new->name = strdup("2/m");
-	return expand_ops(new);
-}
-
-
-static SymOpList *make_2_uab()
-{
-	SymOpList *new = new_symoplist();
-	add_symop(new, v(-1,0,0,0), v(0,1,0,0), v(0,0,0,-1), 2); /* 2 // k */
-	new->name = strdup("2_uab");
-	return expand_ops(new);
-}
-
-
-static SymOpList *make_2()
-{
-	SymOpList *new = new_symoplist();
-	add_symop(new, v(-1,0,0,0), v(0,-1,0,0), v(0,0,0,1), 2); /* 2 // l */
-	new->name = strdup("2");
-	return expand_ops(new);
-}
-
-
-static SymOpList *make_m()
-{
-	SymOpList *new = new_symoplist();
-	add_symop(new, v(1,0,0,0), v(0,1,0,0), v(0,0,0,-1), 2);  /* m -| l */
-	new->name = strdup("m");
-	return expand_ops(new);
-}
-
-
-/******************************** Orthorhombic ********************************/
-
-static SymOpList *make_mmm()
-{
-	SymOpList *new = new_symoplist();
-	add_symop(new, v(-1,0,0,0), v(0,-1,0,0), v(0,0,0,1), 2); /* 2 // l */
-	add_symop(new, v(-1,0,0,0), v(0,1,0,0), v(0,0,0,-1), 2); /* 2 // k */
-	add_symop(new, v(1,0,0,0), v(0,-1,0,0), v(0,0,0,1), 2);  /* m -| k */
-	new->name = strdup("mmm");
-	return expand_ops(new);
-}
-
-
-static SymOpList *make_222()
-{
-	SymOpList *new = new_symoplist();
-	add_symop(new, v(-1,0,0,0), v(0,-1,0,0), v(0,0,0,1), 2); /* 2 // l */
-	add_symop(new, v(-1,0,0,0), v(0,1,0,0), v(0,0,0,-1), 2); /* 2 // k */
-	new->name = strdup("222");
-	return expand_ops(new);
-}
-
-
-static SymOpList *make_mm2()
-{
-	SymOpList *new = new_symoplist();
-	add_symop(new, v(-1,0,0,0), v(0,-1,0,0), v(0,0,0,1), 2); /* 2 // l */
-	add_symop(new, v(1,0,0,0), v(0,-1,0,0), v(0,0,0,1), 2);  /* m -| k */
-	new->name = strdup("mm2");
-	return expand_ops(new);
-}
-
-
-/********************************* Tetragonal *********************************/
-
-static SymOpList *make_4m()
-{
-	SymOpList *new = new_symoplist();
-	add_symop(new, v(0,-1,0,0), v(1,0,0,0), v(0,0,0,1), 4); /* 4 // l */
-	add_symop(new, v(1,0,0,0), v(0,1,0,0), v(0,0,0,-1), 2); /* m -| l */
-	new->name = strdup("4/m");
-	return expand_ops(new);
-}
-
-
-static SymOpList *make_4()
-{
-	SymOpList *new = new_symoplist();
-	add_symop(new, v(0,-1,0,0), v(1,0,0,0), v(0,0,0,1), 4); /* 4 // l */
-	new->name = strdup("4");
-	return expand_ops(new);
-}
-
-
-static SymOpList *make_4mm()
-{
-	SymOpList *new = new_symoplist();
-	add_symop(new, v(0,-1,0,0), v(1,0,0,0), v(0,0,0,1), 4); /* 4 // l */
-	add_symop(new, v(-1,0,0,0), v(0,1,0,0), v(0,0,0,1), 2); /* m -| l */
-	new->name = strdup("4mm");
-	return expand_ops(new);
-}
-
-
-static SymOpList *make_422()
-{
-	SymOpList *new = new_symoplist();
-	add_symop(new, v(0,-1,0,0), v(1,0,0,0), v(0,0,0,1), 4);  /* 4 // l */
-	add_symop(new, v(-1,0,0,0), v(0,1,0,0), v(0,0,0,-1), 2); /* 2 // k */
-	new->name = strdup("422");
-	return expand_ops(new);
-}
-
-
-static SymOpList *make_4bar()
-{
-	SymOpList *new = new_symoplist();
-	add_symop(new, v(0,1,0,0), v(-1,0,0,0), v(0,0,0,-1), 4); /* -4 // l */
-	new->name = strdup("-4");
-	return expand_ops(new);
-}
-
-
-static SymOpList *make_4bar2m()
-{
-	SymOpList *new = new_symoplist();
-	add_symop(new, v(0,1,0,0), v(-1,0,0,0), v(0,0,0,-1), 4); /* -4 // l */
-	add_symop(new, v(-1,0,0,0), v(0,1,0,0), v(0,0,0,-1), 2); /* 2 // k */
-	new->name = strdup("-42m");
-	return expand_ops(new);
-}
-
-
-static SymOpList *make_4barm2()
-{
-	SymOpList *new = new_symoplist();
-	add_symop(new, v(0,1,0,0), v(-1,0,0,0), v(0,0,0,-1), 4); /* -4 // l */
-	add_symop(new, v(0,1,0,0), v(1,0,0,0), v(0,0,0,-1), 2); /* 2 // h+k */
-	new->name = strdup("-4m2");
-	return expand_ops(new);
-}
-
-
-static SymOpList *make_4mmm()
-{
-	SymOpList *new = new_symoplist();
-	add_symop(new, v(0,-1,0,0), v(1,0,0,0), v(0,0,0,1), 4); /* 4 // l */
-	add_symop(new, v(-1,0,0,0), v(0,1,0,0), v(0,0,0,1), 2); /* m -| k */
-	add_symop(new, v(1,0,0,0), v(0,1,0,0), v(0,0,0,-1), 2); /* m -| l */
-	new->name = strdup("4/mmm");
-	return expand_ops(new);
-}
-
-
-/************************** Trigonal (Rhombohedral) ***************************/
-
-static SymOpList *make_3_R()
-{
-	SymOpList *new = new_symoplist();
-	add_symop(new, v(0,0,0,1), v(1,0,0,0), v(0,1,0,0), 3); /* 3 // h+k+l */
-	new->name = strdup("3_R");
-	return expand_ops(new);
-}
-
-
-static SymOpList *make_3bar_R()
-{
-	SymOpList *new = new_symoplist();
-	add_symop(new, v(0,0,0,1), v(1,0,0,0), v(0,1,0,0), 3); /* -3 // h+k+l */
-	add_symop(new, v(-1,0,0,0), v(0,-1,0,0), v(0,0,0,-1), 2); /* -I */
-	new->name = strdup("-3_R");
-	return expand_ops(new);
-}
-
-
-static SymOpList *make_32_R()
-{
-	SymOpList *new = new_symoplist();
-	add_symop(new, v(0,0,0,1), v(1,0,0,0), v(0,1,0,0), 3); /* 3 // h+k+l */
-	add_symop(new, v(0,-1,0,0), v(-1,0,0,0), v(0,0,0,-1), 2); /* 2 -| 3 */
-	new->name = strdup("32_R");
-	return expand_ops(new);
-}
-
-
-static SymOpList *make_3m_R()
-{
-	SymOpList *new = new_symoplist();
-	add_symop(new, v(0,0,0,1), v(1,0,0,0), v(0,1,0,0), 3); /* 3 // h+k+l */
-	add_symop(new, v(0,1,0,0), v(1,0,0,0), v(0,0,0,1), 2); /* m */
-	new->name = strdup("3m_R");
-	return expand_ops(new);
-}
-
-
-static SymOpList *make_3barm_R()
-{
-	SymOpList *new = new_symoplist();
-	add_symop(new, v(0,0,0,1), v(1,0,0,0), v(0,1,0,0), 3); /* -3 // h+k+l */
-	add_symop(new, v(-1,0,0,0), v(0,-1,0,0), v(0,0,0,-1), 2); /* -I */
-	add_symop(new, v(0,1,0,0), v(1,0,0,0), v(0,0,0,1), 2);    /* m */
-	new->name = strdup("-3m_R");
-	return expand_ops(new);
-}
-
-
-/*************************** Trigonal (Hexagonal) *****************************/
-
-static SymOpList *make_3_H()
-{
-	SymOpList *new = new_symoplist();
-	add_symop(new, v(0,0,1,0), v(1,0,0,0), v(0,0,0,1), 3); /* 3 // l */
-	new->name = strdup("3_H");
-	return expand_ops(new);
-}
-
-
-static SymOpList *make_3bar_H()
-{
-	SymOpList *new = new_symoplist();
-	add_symop(new, v(0,0,1,0), v(1,0,0,0), v(0,0,0,1), 3);    /* 3 // l */
-	add_symop(new, v(-1,0,0,0), v(0,-1,0,0), v(0,0,0,-1), 2); /* -I */
-	new->name = strdup("-3_H");
-	return expand_ops(new);
-}
-
-
-static SymOpList *make_321_H()
-{
-	SymOpList *new = new_symoplist();
-	add_symop(new, v(0,0,1,0), v(1,0,0,0), v(0,0,0,1), 3);  /* 3 // l */
-	add_symop(new, v(0,1,0,0), v(1,0,0,0), v(0,0,0,-1), 2); /* 2 // h */
-	new->name = strdup("321_H");
-	return expand_ops(new);
-}
-
-
-static SymOpList *make_312_H()
-{
-	SymOpList *new = new_symoplist();
-	add_symop(new, v(0,0,1,0), v(1,0,0,0), v(0,0,0,1), 3);    /* 3 // l */
-	add_symop(new, v(0,-1,0,0), v(-1,0,0,0), v(0,0,0,-1), 2); /* 2 // h+k */
-	new->name = strdup("312_H");
-	return expand_ops(new);
-}
-
-
-static SymOpList *make_3m1_H()
-{
-	SymOpList *new = new_symoplist();
-	add_symop(new, v(0,0,1,0), v(1,0,0,0), v(0,0,0,1), 3); /* 3 // l */
-	add_symop(new, v(0,-1,0,0), v(-1,0,0,0), v(0,0,0,1), 2); /* m -| i */
-	new->name = strdup("3m1_H");
-	return expand_ops(new);
-}
-
-
-static SymOpList *make_31m_H()
-{
-	SymOpList *new = new_symoplist();
-	add_symop(new, v(0,0,1,0), v(1,0,0,0), v(0,0,0,1), 3); /* 3 // l */
-	add_symop(new, v(0,1,0,0), v(1,0,0,0), v(0,0,0,1), 2); /* m -| (k+i) */
-	new->name = strdup("31m_H");
-	return expand_ops(new);
-}
-
-
-static SymOpList *make_3barm1_H()
-{
-	SymOpList *new = new_symoplist();
-	add_symop(new, v(0,0,1,0), v(1,0,0,0), v(0,0,0,1), 3);    /* 3 // l */
-	add_symop(new, v(-1,0,0,0), v(0,-1,0,0), v(0,0,0,-1), 2); /* -I */
-	add_symop(new, v(0,1,0,0), v(1,0,0,0), v(0,0,0,-1), 2);   /* 2 // h */
-	new->name = strdup("-3m1_H");
-	return expand_ops(new);
-}
-
-
-static SymOpList *make_3bar1m_H()
-{
-	SymOpList *new = new_symoplist();
-	add_symop(new, v(0,0,1,0), v(1,0,0,0), v(0,0,0,1), 3);    /* 3 // l */
-	add_symop(new, v(-1,0,0,0), v(0,-1,0,0), v(0,0,0,-1), 2); /* -I */
-	add_symop(new, v(0,-1,0,0), v(-1,0,0,0), v(0,0,0,-1), 2); /* 2 // h+k */
-	new->name = strdup("-31m_H");
-	return expand_ops(new);
-}
-
-
-/********************************** Hexgonal **********************************/
-
-static SymOpList *make_6()
-{
-	SymOpList *new = new_symoplist();
-	add_symop(new, v(0,0,-1,0), v(-1,0,0,0), v(0,0,0,1), 6); /* 6 // l */
-	new->name = strdup("6");
-	return expand_ops(new);
-}
-
-
-static SymOpList *make_6bar()
-{
-	SymOpList *new = new_symoplist();
-	add_symop(new, v(0,0,1,0), v(1,0,0,0), v(0,0,0,-1), 6); /* -6 // l */
-	new->name = strdup("-6");
-	return expand_ops(new);
-}
-
-
-static SymOpList *make_6m()
-{
-	SymOpList *new = new_symoplist();
-	add_symop(new, v(0,0,-1,0), v(-1,0,0,0), v(0,0,0,1), 6); /* 6 // l */
-	add_symop(new, v(1,0,0,0), v(0,1,0,0), v(0,0,0,-1), 2);  /* m -| l */
-	new->name = strdup("6/m");
-	return expand_ops(new);
-}
-
-
-static SymOpList *make_622()
-{
-	SymOpList *new = new_symoplist();
-	add_symop(new, v(0,0,-1,0), v(-1,0,0,0), v(0,0,0,1), 6); /* 6 // l */
-	add_symop(new, v(0,1,0,0), v(1,0,0,0), v(0,0,0,-1), 2);   /* 2 // h */
-	new->name = strdup("622");
-	return expand_ops(new);
-}
-
-
-static SymOpList *make_6mm()
-{
-	SymOpList *new = new_symoplist();
-	add_symop(new, v(0,0,-1,0), v(-1,0,0,0), v(0,0,0,1), 6); /* 6 // l */
-	add_symop(new, v(0,-1,0,0), v(-1,0,0,0), v(0,0,0,1), 2); /* m -| i */
-	new->name = strdup("6mm");
-	return expand_ops(new);
-}
-
-
-static SymOpList *make_6barm2()
-{
-	SymOpList *new = new_symoplist();
-	add_symop(new, v(0,0,1,0), v(1,0,0,0), v(0,0,0,-1), 6); /* -6 // l */
-	add_symop(new, v(0,-1,0,0), v(-1,0,0,0), v(0,0,0,1), 2); /* m -| i */
-	new->name = strdup("-6m2");
-	return expand_ops(new);
-}
-
-
-static SymOpList *make_6bar2m()
-{
-	SymOpList *new = new_symoplist();
-	add_symop(new, v(0,0,1,0), v(1,0,0,0), v(0,0,0,-1), 6); /* -6 // l */
-	add_symop(new, v(0,1,0,0), v(1,0,0,0), v(0,0,0,1), 2);  /* m -| (k+i) */
-	new->name = strdup("-62m");
-	return expand_ops(new);
-}
-
-
-static SymOpList *make_6mmm()
-{
-	SymOpList *new = new_symoplist();
-	add_symop(new, v(0,0,1,0), v(1,0,0,0), v(0,0,0,-1), 6); /* -6 // l */
-	add_symop(new, v(0,-1,0,0), v(-1,0,0,0), v(0,0,0,1), 2); /* m -| i */
-	add_symop(new, v(-1,0,0,0), v(0,-1,0,0), v(0,0,0,-1), 2); /* -I */
-	new->name = strdup("6/mmm");
-	return expand_ops(new);
-}
-
-
-/************************************ Cubic ***********************************/
-
-static SymOpList *make_23()
-{
-	SymOpList *new = new_symoplist();
-	add_symop(new, v(-1,0,0,0), v(0,-1,0,0), v(0,0,0,1), 2); /* 2// l */
-	add_symop(new, v(-1,0,0,0), v(0,1,0,0), v(0,0,0,-1), 2); /* 2// k */
-	add_symop(new, v(0,1,0,0), v(0,0,0,1), v(1,0,0,0), 3); /* 3// h+k+l */
-	new->name = strdup("23");
-	return expand_ops(new);
-}
-
-
-static SymOpList *make_m3bar()
-{
-	SymOpList *new = new_symoplist();
-	add_symop(new, v(-1,0,0,0), v(0,-1,0,0), v(0,0,0,1), 2); /* 2// l */
-	add_symop(new, v(-1,0,0,0), v(0,1,0,0), v(0,0,0,-1), 2); /* 2// k */
-	add_symop(new, v(0,1,0,0), v(0,0,0,1), v(1,0,0,0), 3); /* 3// h+k+l */
-	add_symop(new, v(-1,0,0,0), v(0,-1,0,0), v(0,0,0,-1), 2); /* -I */
-	new->name = strdup("m-3");
-	return expand_ops(new);
-}
-
-
-static SymOpList *make_432()
-{
-	SymOpList *new = new_symoplist();
-	add_symop(new, v(0,-1,0,0), v(1,0,0,0), v(0,0,0,1), 4); /* 4 // l */
-	add_symop(new, v(-1,0,0,0), v(0,1,0,0), v(0,0,0,-1), 2);/* 2 // k */
-	add_symop(new, v(0,1,0,0), v(0,0,0,1), v(1,0,0,0), 3);  /* 3 // h+k+l */
-	new->name = strdup("432");
-	return expand_ops(new);
-}
-
-
-static SymOpList *make_4bar3m()
-{
-	SymOpList *new = new_symoplist();
-	add_symop(new, v(0,1,0,0), v(-1,0,0,0), v(0,0,0,-1), 4); /* -4 // l */
-	add_symop(new, v(-1,0,0,0), v(0,1,0,0), v(0,0,0,-1), 2);/* 2 // k */
-	add_symop(new, v(0,1,0,0), v(0,0,0,1), v(1,0,0,0), 3);  /* 3 // h+k+l */
-	new->name = strdup("-43m");
-	return expand_ops(new);
-}
-
-
-static SymOpList *make_m3barm()
-{
-	SymOpList *new = new_symoplist();
-	add_symop(new, v(0,-1,0,0), v(1,0,0,0), v(0,0,0,1), 4); /* 4 // l */
-	add_symop(new, v(-1,0,0,0), v(0,1,0,0), v(0,0,0,-1), 2);/* 2 // k */
-	add_symop(new, v(0,1,0,0), v(0,0,0,1), v(1,0,0,0), 3);  /* 3 // h+k+l */
-	add_symop(new, v(-1,0,0,0), v(0,-1,0,0), v(0,0,0,-1), 2); /* -I */
-	new->name = strdup("m-3m");
-	return expand_ops(new);
-}
-
-
-/**
- * get_pointgroup:
- * @sym: A string representation of a point group
- *
- * This function parses @sym and returns the corresponding %SymOpList.
- * In the string representation of the point group, use a preceding minus sign
- * for any character which would have a "bar".  Trigonal groups must be suffixed
- * with either "_H" or "_R" for a hexagonal or rhombohedral lattice
- * respectively.
- *
- * Examples: -1 1 2/m 2 m mmm 222 mm2 4/m 4 -4 4/mmm 422 -42m -4m2 4mm
- * 3_R -3_R 32_R 3m_R -3m_R 3_H -3_H 321_H 312_H 3m1_H 31m_H -3m1_H -31m_H
- * 6/m 6 -6 6/mmm 622 -62m -6m2 6mm 23 m-3 432 -43m m-3m.
- **/
-SymOpList *get_pointgroup(const char *sym)
-{
-	/* Triclinic */
-	if ( strcmp(sym, "-1") == 0 ) return make_1bar();
-	if ( strcmp(sym, "1") == 0 ) return make_1();
-
-	/* Monoclinic */
-	if ( strcmp(sym, "2/m") == 0 ) return make_2m();
-	if ( strcmp(sym, "2") == 0 ) return make_2();
-	if ( strcmp(sym, "2_uab") == 0 ) return make_2_uab();
-	if ( strcmp(sym, "m") == 0 ) return make_m();
-
-	/* Orthorhombic */
-	if ( strcmp(sym, "mmm") == 0 ) return make_mmm();
-	if ( strcmp(sym, "222") == 0 ) return make_222();
-	if ( strcmp(sym, "mm2") == 0 ) return make_mm2();
-
-	/* Tetragonal */
-	if ( strcmp(sym, "4/m") == 0 ) return make_4m();
-	if ( strcmp(sym, "4") == 0 ) return make_4();
-	if ( strcmp(sym, "-4") == 0 ) return make_4bar();
-	if ( strcmp(sym, "4/mmm") == 0 ) return make_4mmm();
-	if ( strcmp(sym, "422") == 0 ) return make_422();
-	if ( strcmp(sym, "-42m") == 0 ) return make_4bar2m();
-	if ( strcmp(sym, "-4m2") == 0 ) return make_4barm2();
-	if ( strcmp(sym, "4mm") == 0 ) return make_4mm();
-
-	/* Trigonal (rhombohedral) */
-	if ( strcmp(sym, "3_R") == 0 ) return make_3_R();
-	if ( strcmp(sym, "-3_R") == 0 ) return make_3bar_R();
-	if ( strcmp(sym, "32_R") == 0 ) return make_32_R();
-	if ( strcmp(sym, "3m_R") == 0 ) return make_3m_R();
-	if ( strcmp(sym, "-3m_R") == 0 ) return make_3barm_R();
-
-	/* Trigonal (hexagonal) */
-	if ( strcmp(sym, "3_H") == 0 ) return make_3_H();
-	if ( strcmp(sym, "-3_H") == 0 ) return make_3bar_H();
-	if ( strcmp(sym, "321_H") == 0 ) return make_321_H();
-	if ( strcmp(sym, "312_H") == 0 ) return make_312_H();
-	if ( strcmp(sym, "3m1_H") == 0 ) return make_3m1_H();
-	if ( strcmp(sym, "31m_H") == 0 ) return make_31m_H();
-	if ( strcmp(sym, "-3m1_H") == 0 ) return make_3barm1_H();
-	if ( strcmp(sym, "-31m_H") == 0 ) return make_3bar1m_H();
-
-	/* Hexagonal */
-	if ( strcmp(sym, "6/m") == 0 ) return make_6m();
-	if ( strcmp(sym, "6") == 0 ) return make_6();
-	if ( strcmp(sym, "-6") == 0 ) return make_6bar();
-	if ( strcmp(sym, "6/mmm") == 0 ) return make_6mmm();
-	if ( strcmp(sym, "622") == 0 ) return make_622();
-	if ( strcmp(sym, "-62m") == 0 ) return make_6bar2m();
-	if ( strcmp(sym, "-6m2") == 0 ) return make_6barm2();
-	if ( strcmp(sym, "6mm") == 0 ) return make_6mm();
-
-	/* Cubic */
-	if ( strcmp(sym, "23") == 0 ) return make_23();
-	if ( strcmp(sym, "m-3") == 0 ) return make_m3bar();
-	if ( strcmp(sym, "432") == 0 ) return make_432();
-	if ( strcmp(sym, "-43m") == 0 ) return make_4bar3m();
-	if ( strcmp(sym, "m-3m") == 0 ) return make_m3barm();
-
-	ERROR("Unknown point group '%s'\n", sym);
-	return NULL;
-}
-
-
-static void do_op(const struct sym_op *op,
-                  signed int h, signed int k, signed int l,
-                  signed int *he, signed int *ke, signed int *le)
-{
-	*he = h*op->h[0] + k*op->h[1] + l*op->h[2];
-	*ke = h*op->k[0] + k*op->k[1] + l*op->k[2];
-	*le = h*op->l[0] + k*op->l[1] + l*op->l[2];
-}
-
-
-/**
- * get_equiv:
- * @ops: A %SymOpList
- * @m: A %SymOpMask, which has been shown to special_position()
- * @idx: Index of the operation to use
- * @h: index of reflection
- * @k: index of reflection
- * @l: index of reflection
- * @he: location to store h index of equivalent reflection
- * @ke: location to store k index of equivalent reflection
- * @le: location to store l index of equivalent reflection
- *
- * This function applies the @idx-th symmetry operation from @ops to the
- * reflection @h, @k, @l, and stores the result at @he, @ke and @le.
- *
- * If you don't mind that the same equivalent might appear twice, simply call
- * this function the number of times returned by num_ops(), using the actual
- * point group.  If repeating the same equivalent twice (for example, if the
- * given reflection is a special high-symmetry one), call special_position()
- * first to get a "specialised" SymOpList and use that instead.
- **/
-void get_equiv(const SymOpList *ops, const SymOpMask *m, int idx,
-               signed int h, signed int k, signed int l,
-               signed int *he, signed int *ke, signed int *le)
-{
-	const int n = num_ops(ops);
-
-	if ( m != NULL ) {
-
-		int i, c;
-
-		c = 0;
-		for ( i=0; i<n; i++ ) {
-
-			if ( (c == idx) && m->mask[i] ) {
-				do_op(&ops->ops[i], h, k, l, he, ke, le);
-				return;
-			}
-
-			if ( m->mask[i] ) {
-				c++;
-			}
-
-		}
-
-		ERROR("Index %i out of range for point group '%s' with"
-		      " reflection %i %i %i\n",
-		      idx, symmetry_name(ops), h, k, l);
-
-		*he = 0;  *ke = 0;  *le = 0;
-
-		return;
-
-	}
-
-
-
-	if ( idx >= n ) {
-
-		ERROR("Index %i out of range for point group '%s'\n", idx,
-		      symmetry_name(ops));
-
-		*he = 0;  *ke = 0;  *le = 0;
-		return;
-
-	}
-
-	do_op(&ops->ops[idx], h, k, l, he, ke, le);
-}
-
-
-/**
- * special_position:
- * @ops: A %SymOpList, usually corresponding to a point group
- * @m: A %SymOpMask created with new_symopmask()
- * @h: index of a reflection
- * @k: index of a reflection
- * @l: index of a reflection
- *
- * This function determines which operations in @ops map the reflection @h, @k,
- * @l onto itself, and uses @m to flag the operations in @ops which cause this.
- *
- **/
-void special_position(const SymOpList *ops, SymOpMask *m,
-                      signed int h, signed int k, signed int l)
-{
-	int i, n;
-	signed int *htest;
-	signed int *ktest;
-	signed int *ltest;
-
-	assert(m->list = ops);
-
-	n = num_equivs(ops, NULL);
-	htest = malloc(n*sizeof(signed int));
-	ktest = malloc(n*sizeof(signed int));
-	ltest = malloc(n*sizeof(signed int));
-
-	for ( i=0; i<n; i++ ) {
-
-		signed int he, ke, le;
-		int j;
-
-		get_equiv(ops, NULL, i, h, k, l, &he, &ke, &le);
-
-		m->mask[i] = 1;
-		for ( j=0; j<i; j++ ) {
-			if ( (he==htest[j]) && (ke==ktest[j])
-			  && (le==ltest[j]) )
-			{
-				m->mask[i] = 0;
-				break;  /* Only need to find one */
-			}
-		}
-
-		htest[i] = he;
-		ktest[i] = ke;
-		ltest[i] = le;
-
-	}
-
-	free(htest);
-	free(ktest);
-	free(ltest);
-}
-
-
-static int any_negative(signed int h, signed int k, signed int l)
-{
-	if ( h < 0 ) return 1;
-	if ( k < 0 ) return 1;
-	if ( l < 0 ) return 1;
-	return 0;
-}
-
-
-/**
- * get_asymm:
- * @ops: A %SymOpList, usually corresponding to a point group
- * @h: index of a reflection
- * @k: index of a reflection
- * @l: index of a reflection
- * @hp: location for asymmetric index of reflection
- * @kp: location for asymmetric index of reflection
- * @lp: location for asymmetric index of reflection
- *
- * This function determines the asymmetric version of the reflection @h, @k, @l
- * in symmetry group @ops, and puts the result in @hp, @kp, @lp.
- *
- * This is a relatively expensive operation because of its generality.
- * Therefore, if you know you'll need to make repeated use of the asymmetric
- * indices, consider creating a new %RefList indexed according to the asymmetric
- * indices themselves with asymmetric_indices().  If you do that, you'll still
- * be able to get the original versions of the indices with
- * get_symmetric_indices().
- *
- **/
-void get_asymm(const SymOpList *ops,
-               signed int h, signed int k, signed int l,
-               signed int *hp, signed int *kp, signed int *lp)
-{
-	int nequiv;
-	int p;
-	signed int best_h, best_k, best_l;
-	int have_negs;
-
-	nequiv = num_equivs(ops, NULL);
-
-	best_h = h;  best_k = k;  best_l = l;
-	have_negs = any_negative(best_h, best_k, best_l);
-	for ( p=0; p<nequiv; p++ ) {
-
-		int will_have_negs;
-
-		get_equiv(ops, NULL, p, h, k, l, hp, kp, lp);
-
-		will_have_negs = any_negative(*hp, *kp, *lp);
-
-		/* Don't lose "no negs" status */
-		if ( !have_negs && will_have_negs ) continue;
-
-		if ( have_negs && !will_have_negs ) {
-			best_h = *hp;  best_k = *kp;  best_l = *lp;
-			have_negs = 0;
-			continue;
-		}
-
-		if ( *hp > best_h ) {
-			best_h = *hp;  best_k = *kp;  best_l = *lp;
-			have_negs = any_negative(best_h, best_k, best_l);
-			continue;
-		}
-		if ( *hp < best_h ) continue;
-
-		if ( *kp > best_k ) {
-			best_h = *hp;  best_k = *kp;  best_l = *lp;
-			have_negs = any_negative(best_h, best_k, best_l);
-			continue;
-		}
-		if ( *kp < best_k ) continue;
-
-		if ( *lp > best_l ) {
-			best_h = *hp;  best_k = *kp;  best_l = *lp;
-			have_negs = any_negative(best_h, best_k, best_l);
-			continue;
-		}
-
-	}
-
-	*hp = best_h;  *kp = best_k;  *lp = best_l;
-}
-
-
-static int is_inversion(const struct sym_op *op)
-{
-	if ( (op->h[0]!=-1) || (op->h[1]!=0) || (op->h[2]!=0) ) return 0;
-	if ( (op->k[0]!=0) || (op->k[1]!=-1) || (op->k[2]!=0) ) return 0;
-	if ( (op->l[0]!=0) || (op->l[1]!=0) || (op->l[2]!=-1) ) return 0;
-	return 1;
-}
-
-
-static int is_identity(const struct sym_op *op)
-{
-	if ( (op->h[0]!=1) || (op->h[1]!=0) || (op->h[2]!=0) ) return 0;
-	if ( (op->k[0]!=0) || (op->k[1]!=1) || (op->k[2]!=0) ) return 0;
-	if ( (op->l[0]!=0) || (op->l[1]!=0) || (op->l[2]!=1) ) return 0;
-	return 1;
-}
-
-
-static signed int determinant(const struct sym_op *op)
-{
-	signed int det = 0;
-
-	det += op->h[0] * (op->k[1]*op->l[2] - op->k[2]*op->l[1]);
-	det -= op->h[1] * (op->k[0]*op->l[2] - op->k[2]*op->l[0]);
-	det += op->h[2] * (op->k[0]*op->l[1] - op->k[1]*op->l[0]);
-
-	return det;
-}
-
-
-/**
- * is_centrosymmetric:
- * @s: A %SymOpList
- *
- * Returns: non-zero if @s contains an inversion operation
- */
-int is_centrosymmetric(const SymOpList *s)
-{
-	int i, n;
-
-	n = num_ops(s);
-	for ( i=0; i<n; i++ ) {
-		if ( is_inversion(&s->ops[i]) ) return 1;
-	}
-
-	return 0;
-}
-
-
-static int ops_equal(const struct sym_op *op,
-                     signed int *h, signed int *k, signed int *l)
-{
-	if ( (op->h[0]!=h[0]) || (op->h[1]!=h[1]) || (op->h[2]!=h[2]) )
-		return 0;
-	if ( (op->k[0]!=k[0]) || (op->k[1]!=k[1]) || (op->k[2]!=k[2]) )
-		return 0;
-	if ( (op->l[0]!=l[0]) || (op->l[1]!=l[1]) || (op->l[2]!=l[2]) )
-		return 0;
-	return 1;
-}
-
-
-static int struct_ops_equal(const struct sym_op *op1, const struct sym_op *op2)
-{
-	return ops_equal(op1, op2->h, op2->k, op2->l);
-}
-
-
-/* Return true if a*b = ans */
-static int check_mult(const struct sym_op *ans,
-                      const struct sym_op *a, const struct sym_op *b)
-{
-	signed int *ans_h, *ans_k, *ans_l;
-	int val;
-
-	ans_h = malloc(3*sizeof(signed int));
-	ans_k = malloc(3*sizeof(signed int));
-	ans_l = malloc(3*sizeof(signed int));
-
-	combine_ops(a->h, a->k, a->l, b->h, b->k, b->l, ans_h, ans_k, ans_l);
-	val = ops_equal(ans, ans_h, ans_k, ans_l);
-
-	free(ans_h);
-	free(ans_k);
-	free(ans_l);
-
-	return val;
-}
-
-
-/**
- * is_subgroup:
- * @source: A %SymOpList
- * @target: Another %SymOpList, which might be a subgroup of @source.
- *
- * Returns: non-zero if every operation in @target is also in @source.
- **/
-int is_subgroup(const SymOpList *source, const SymOpList *target)
-{
-	int n_src, n_tgt;
-	int i;
-
-	n_src = num_ops(source);
-	n_tgt = num_ops(target);
-
-	for ( i=0; i<n_tgt; i++ ) {
-
-		int j;
-		int found = 0;
-
-		for ( j=0; j<n_src; j++ ) {
-
-			if ( struct_ops_equal(&target->ops[i],
-			                      &source->ops[j] ) )
-			{
-				found = 1;
-				break;
-			}
-
-		}
-
-		if ( !found ) return 0;
-
-	}
-
-	return 1;
-}
-
-
-/**
- * get_ambiguities:
- * @source: The "source" symmetry, a %SymOpList
- * @target: The "target" symmetry, a %SymOpList
-
- * Calculates twinning laws.  Returns a %SymOpList containing the twinning
- * operators, which are the symmetry operations which can be added to @target
- * to generate @source.  Only rotations are allowable - no mirrors nor
- * inversions.
- * To count the number of possibilities, use num_ops() on the result.
- *
- * Returns: A %SymOpList containing the twinning operators, or NULL if the
- * source symmetry cannot be generated from that target symmetry without using
- * mirror or inversion operations.
- */
-SymOpList *get_ambiguities(const SymOpList *source, const SymOpList *target)
-{
-	int n_src, n_tgt;
-	int i;
-	SymOpList *twins;
-	SymOpList *src_reordered;
-	SymOpMask *used;
-	char *name;
-	int index;
-
-	n_src = num_ops(source);
-	n_tgt = num_ops(target);
-
-	if ( !is_subgroup(source, target) ) {
-		ERROR("'%s' is not a subgroup of '%s'\n",
-		      symmetry_name(target), symmetry_name(source));
-		return NULL;
-	}
-
-	if ( n_src % n_tgt != 0 ) {
-		ERROR("Subgroup index would be fractional.\n");
-		return NULL;
-	}
-	index = n_src / n_tgt;
-
-	src_reordered = new_symoplist();
-	used = new_symopmask(source);
-
-	/* Find identity */
-	for ( i=0; i<n_src; i++ ) {
-		if ( used->mask[i] == 0 ) continue;
-		if ( is_identity(&source->ops[i]) ) {
-			add_copied_op(src_reordered, &source->ops[i]);
-			used->mask[i] = 0;
-		}
-	}
-
-	/* Find binary options (order=2) of first kind (determinant positive) */
-	for ( i=0; i<n_src; i++ ) {
-		if ( used->mask[i] == 0 ) continue;
-		if ( (source->ops[i].order == 2)
-		  && (determinant(&source->ops[i]) > 0) ) {
-			add_copied_op(src_reordered, &source->ops[i]);
-			used->mask[i] = 0;
-		}
-	}
-
-	/* Find other operations of first kind (determinant positive) */
-	for ( i=0; i<n_src; i++ ) {
-		if ( used->mask[i] == 0 ) continue;
-		if ( determinant(&source->ops[i]) > 0 ) {
-			add_copied_op(src_reordered, &source->ops[i]);
-			used->mask[i] = 0;
-		}
-	}
-
-	/* Find inversion */
-	for ( i=0; i<n_src; i++ ) {
-		if ( used->mask[i] == 0 ) continue;
-		if ( is_inversion(&source->ops[i]) ) {
-			add_copied_op(src_reordered, &source->ops[i]);
-			used->mask[i] = 0;
-		}
-	}
-
-	/* Find binary options of second kind (determinant negative) */
-	for ( i=0; i<n_src; i++ ) {
-		if ( used->mask[i] == 0 ) continue;
-		if ( (source->ops[i].order == 2)
-		  && (determinant(&source->ops[i]) < 0) ) {
-			add_copied_op(src_reordered, &source->ops[i]);
-			used->mask[i] = 0;
-		}
-	}
-
-	/* Find other operations of second kind (determinant negative) */
-	for ( i=0; i<n_src; i++ ) {
-		if ( used->mask[i] == 0 ) continue;
-		if ( determinant(&source->ops[i]) < 0 ) {
-			add_copied_op(src_reordered, &source->ops[i]);
-			used->mask[i] = 0;
-		}
-	}
-
-	int n_left_over = 0;
-	for ( i=0; i<n_src; i++ ) {
-		if ( used->mask[i] == 0 ) continue;
-		n_left_over++;
-	}
-	if ( n_left_over != 0 ) {
-		ERROR("%i operations left over after rearranging for"
-		      " left coset decomposition.\n", n_left_over);
-	}
-
-	if ( num_ops(src_reordered) != num_ops(source) ) {
-		ERROR("%i ops went to %i after rearranging.\n",
-		      num_ops(src_reordered), num_ops(source));
-	}
-
-	free_symopmask(used);
-	used = new_symopmask(src_reordered);
-
-	/* This is the first method from Flack (1987) */
-	for ( i=0; i<n_src; i++ ) {
-
-		int j;
-		if ( used->mask[i] == 0 ) continue;
-
-		for ( j=1; j<n_tgt; j++ ) {
-
-			int k;
-			for ( k=i+1; k<n_src; k++ ) {
-				if ( check_mult(&src_reordered->ops[k],
-				                &src_reordered->ops[i],
-				                &target->ops[j]) )
-				{
-					used->mask[k] = 0;
-				}
-			}
-
-		}
-
-	}
-
-	twins = new_symoplist();
-	for ( i=0; i<n_src; i++ ) {
-		if ( used->mask[i] == 0 ) continue;
-		if ( determinant(&src_reordered->ops[i]) < 0 ) {
-			/* A mirror or inversion turned up in the list.
-			 * That means that no pure rotational ambiguity can
-			 * account for this subgroup relationship. */
-			free_symoplist(twins);
-			free_symopmask(used);
-			free_symoplist(src_reordered);
-			return NULL;
-		}
-		add_copied_op(twins, &src_reordered->ops[i]);
-	}
-
-	free_symopmask(used);
-	free_symoplist(src_reordered);
-
-	name = malloc(64);
-	snprintf(name, 63, "%s -> %s", symmetry_name(source),
-	                               symmetry_name(target));
-	twins->name = name;
-
-	return twins;
-}
-
-
-static void add_chars(char *t, const char *s, int max_len)
-{
-	char *tmp;
-
-	tmp = strdup(t);
-
-	snprintf(t, max_len, "%s%s", tmp, s);
-	free(tmp);
-}
-
-
-static char *get_matrix_name(signed int *v)
-{
-	char *text;
-	const int max_len = 9;
-	int i;
-	int printed = 0;
-
-	text = malloc(max_len+1);
-	text[0] = '\0';
-
-	for ( i=0; i<3; i++ ) {
-
-		if ( v[i] == 0 ) continue;
-
-		if ( (i==0) && (v[0]==v[1]) ) {
-			if ( v[i]>0 ) add_chars(text, "-", max_len);
-			add_chars(text, "i", max_len);
-			v[1] -= v[0];
-			continue;
-		}
-
-		if ( printed ) add_chars(text, "+", max_len);
-
-		if ( v[i]<0 ) add_chars(text, "-", max_len);
-
-		if ( abs(v[i])>1 ) {
-			char num[3];
-			snprintf(num, 2, "%i", abs(v[i]));
-			add_chars(text, num, max_len);
-		}
-
-		switch ( i )
-		{
-			case 0  : add_chars(text, "h", max_len); break;
-			case 1  : add_chars(text, "k", max_len); break;
-			case 2  : add_chars(text, "l", max_len); break;
-			default : add_chars(text, "X", max_len); break;
-		}
-
-		printed = 1;
-
-	}
-
-	return text;
-}
-
-
-static char *name_equiv(const struct sym_op *op)
-{
-	char *h, *k, *l;
-	char *name;
-
-	h = get_matrix_name(op->h);
-	k = get_matrix_name(op->k);
-	l = get_matrix_name(op->l);
-	name = malloc(32);
-
-	snprintf(name, 31, "%s%s%s", h, k, l);
-	free(h);
-	free(k);
-	free(l);
-
-	return name;
-}
-
-
-/**
- * describe_symmetry:
- * @s: A %SymOpList
- *
- * Writes the name and a list of operations to stderr.
- */
-void describe_symmetry(const SymOpList *s)
-{
-	int i, n;
-
-	n = num_equivs(s, NULL);
-
-	STATUS("%15s :", symmetry_name(s));
-
-	for ( i=0; i<n; i++ ) {
-		char *name = name_equiv(&s->ops[i]);
-		STATUS(" %6s", name);
-		free(name);
-		if ( (i!=0) && (i%8==0) ) STATUS("\n%15s  ", "");
-	}
-	STATUS("\n");
-}
-
-
-/**
- * symmetry_name:
- * @ops: A %SymOpList
- *
- * Returns: a text description of @ops.
- */
-const char *symmetry_name(const SymOpList *ops)
-{
-	return ops->name;
-}
diff --git a/src/symmetry.h b/src/symmetry.h
deleted file mode 100644
index 071ebbde..00000000
--- a/src/symmetry.h
+++ /dev/null
@@ -1,63 +0,0 @@
-/*
- * symmetry.h
- *
- * Symmetry
- *
- * (c) 2006-2010 Thomas White <taw@physics.org>
- *
- * Part of CrystFEL - crystallography with a FEL
- *
- */
-
-
-#ifndef SYMMETRY_H
-#define SYMMETRY_H
-
-#ifdef HAVE_CONFIG_H
-#include <config.h>
-#endif
-
-/**
- * SymOpList
- *
- * The SymOpList is an opaque data structure containing a list of point symmetry
- * operations.  It could represent an point group or a list of indexing
- * ambiguities (twin laws), or similar.
- **/
-typedef struct _symoplist SymOpList;
-
-/**
- * SymOpMask
- *
- * The SymOpMask is an opaque data structure containing a list of flags
- * associated with point symmetry operations in a specific %SymOpList.  It is
- * used to filter the operations in the %SymOpList to avoid duplicating
- * equivalent reflections when the reflection is somehow special (e.g. 'hk0').
- **/
-typedef struct _symopmask SymOpMask;
-
-extern void free_symoplist(SymOpList *ops);
-extern SymOpList *get_pointgroup(const char *sym);
-
-extern SymOpMask *new_symopmask(const SymOpList *list);
-extern void free_symopmask(SymOpMask *m);
-
-extern void special_position(const SymOpList *ops, SymOpMask *m,
-                             signed int h, signed int k, signed int l);
-extern void get_asymm(const SymOpList *ops,
-                      signed int h, signed int k, signed int l,
-                      signed int *hp, signed int *kp, signed int *lp);
-extern int num_equivs(const SymOpList *ops, const SymOpMask *m);
-extern void get_equiv(const SymOpList *ops, const SymOpMask *m, int idx,
-                      signed int h, signed int k, signed int l,
-                      signed int *he, signed int *ke, signed int *le);
-
-extern SymOpList *get_ambiguities(const SymOpList *source,
-                                  const SymOpList *target);
-extern int is_subgroup(const SymOpList *source, const SymOpList *target);
-
-extern int is_centrosymmetric(const SymOpList *s);
-extern const char *symmetry_name(const SymOpList *ops);
-extern void describe_symmetry(const SymOpList *s);
-
-#endif	/* SYMMETRY_H */