Move a load more stuff into libcrystfel

20 files changed, 5832 insertions, 0 deletions

diff --git a/libcrystfel/src/cl-utils.c b/libcrystfel/src/cl-utils.c
new file mode 100644
index 00000000..65a09363
--- /dev/null
+++ b/libcrystfel/src/cl-utils.c
@@ -0,0 +1,200 @@
+/*
+ * 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/libcrystfel/src/cl-utils.h b/libcrystfel/src/cl-utils.h
new file mode 100644
index 00000000..21a7ecd2
--- /dev/null
+++ b/libcrystfel/src/cl-utils.h
@@ -0,0 +1,27 @@
+/*
+ * 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/libcrystfel/src/diffraction-gpu.c b/libcrystfel/src/diffraction-gpu.c
new file mode 100644
index 00000000..605b1514
--- /dev/null
+++ b/libcrystfel/src/diffraction-gpu.c
@@ -0,0 +1,529 @@
+/*
+ * 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/libcrystfel/src/diffraction-gpu.h b/libcrystfel/src/diffraction-gpu.h
new file mode 100644
index 00000000..a3bde4e1
--- /dev/null
+++ b/libcrystfel/src/diffraction-gpu.h
@@ -0,0 +1,57 @@
+/*
+ * 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/libcrystfel/src/diffraction.c b/libcrystfel/src/diffraction.c
new file mode 100644
index 00000000..9532a6ce
--- /dev/null
+++ b/libcrystfel/src/diffraction.c
@@ -0,0 +1,463 @@
+/*
+ * 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/libcrystfel/src/diffraction.h b/libcrystfel/src/diffraction.h
new file mode 100644
index 00000000..f71d3cce
--- /dev/null
+++ b/libcrystfel/src/diffraction.h
@@ -0,0 +1,34 @@
+/*
+ * 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/libcrystfel/src/filters.c b/libcrystfel/src/filters.c
new file mode 100644
index 00000000..c4e409df
--- /dev/null
+++ b/libcrystfel/src/filters.c
@@ -0,0 +1,130 @@
+/*
+ * 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/libcrystfel/src/filters.h b/libcrystfel/src/filters.h
new file mode 100644
index 00000000..6b35d7e8
--- /dev/null
+++ b/libcrystfel/src/filters.h
@@ -0,0 +1,25 @@
+/*
+ * 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/libcrystfel/src/geometry.c b/libcrystfel/src/geometry.c
new file mode 100644
index 00000000..485abba3
--- /dev/null
+++ b/libcrystfel/src/geometry.c
@@ -0,0 +1,341 @@
+/*
+ * 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/libcrystfel/src/geometry.h b/libcrystfel/src/geometry.h
new file mode 100644
index 00000000..ddf04b80
--- /dev/null
+++ b/libcrystfel/src/geometry.h
@@ -0,0 +1,26 @@
+/*
+ * 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/libcrystfel/src/peaks.c b/libcrystfel/src/peaks.c
new file mode 100644
index 00000000..ad524c61
--- /dev/null
+++ b/libcrystfel/src/peaks.c
@@ -0,0 +1,593 @@
+/*
+ * 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 "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/libcrystfel/src/peaks.h b/libcrystfel/src/peaks.h
new file mode 100644
index 00000000..9d475ea9
--- /dev/null
+++ b/libcrystfel/src/peaks.h
@@ -0,0 +1,38 @@
+/*
+ * 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/libcrystfel/src/reflist-utils.c b/libcrystfel/src/reflist-utils.c
new file mode 100644
index 00000000..b64e9979
--- /dev/null
+++ b/libcrystfel/src/reflist-utils.c
@@ -0,0 +1,502 @@
+/*
+ * 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/libcrystfel/src/reflist-utils.h b/libcrystfel/src/reflist-utils.h
new file mode 100644
index 00000000..d14e5f8e
--- /dev/null
+++ b/libcrystfel/src/reflist-utils.h
@@ -0,0 +1,52 @@
+/*
+ * 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/libcrystfel/src/statistics.c b/libcrystfel/src/statistics.c
new file mode 100644
index 00000000..7990672d
--- /dev/null
+++ b/libcrystfel/src/statistics.c
@@ -0,0 +1,668 @@
+/*
+ * 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/libcrystfel/src/statistics.h b/libcrystfel/src/statistics.h
new file mode 100644
index 00000000..8fa78ea6
--- /dev/null
+++ b/libcrystfel/src/statistics.h
@@ -0,0 +1,45 @@
+/*
+ * 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/libcrystfel/src/stream.c b/libcrystfel/src/stream.c
new file mode 100644
index 00000000..a7cdc2d9
--- /dev/null
+++ b/libcrystfel/src/stream.c
@@ -0,0 +1,487 @@
+/*
+ * 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/libcrystfel/src/stream.h b/libcrystfel/src/stream.h
new file mode 100644
index 00000000..ba218fb9
--- /dev/null
+++ b/libcrystfel/src/stream.h
@@ -0,0 +1,49 @@
+/*
+ * 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/libcrystfel/src/symmetry.c b/libcrystfel/src/symmetry.c
new file mode 100644
index 00000000..f0b24146
--- /dev/null
+++ b/libcrystfel/src/symmetry.c
@@ -0,0 +1,1503 @@
+/*
+ * 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/libcrystfel/src/symmetry.h b/libcrystfel/src/symmetry.h
new file mode 100644
index 00000000..071ebbde
--- /dev/null
+++ b/libcrystfel/src/symmetry.h
@@ -0,0 +1,63 @@
+/*
+ * 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 */