Introduce "libcrystfel"

17 files changed, 6352 insertions, 0 deletions

diff --git a/libcrystfel/src/beam-parameters.c b/libcrystfel/src/beam-parameters.c
new file mode 100644
index 00000000..ac30ad04
--- /dev/null
+++ b/libcrystfel/src/beam-parameters.c
@@ -0,0 +1,133 @@
+/*
+ * beam-parameters.c
+ *
+ * Beam parameters
+ *
+ * (c) 2006-2010 Thomas White <taw@physics.org>
+ *
+ * Part of CrystFEL - crystallography with a FEL
+ *
+ */
+
+
+#ifdef HAVE_CONFIG_H
+#include <config.h>
+#endif
+
+#include <stdio.h>
+#include <stdlib.h>
+
+#include "beam-parameters.h"
+#include "utils.h"
+
+
+struct beam_params *get_beam_parameters(const char *filename)
+{
+	FILE *fh;
+	struct beam_params *b;
+	char *rval;
+	int reject;
+
+	fh = fopen(filename, "r");
+	if ( fh == NULL ) return NULL;
+
+	b = calloc(1, sizeof(struct beam_params));
+	if ( b == NULL ) return NULL;
+
+	b->fluence = -1.0;
+	b->beam_radius = -1.0;
+	b->photon_energy = -1.0;
+	b->bandwidth = -1.0;
+	b->divergence = -1.0;
+	b->dqe = -1.0;
+	b->adu_per_photon = -1.0;
+
+	do {
+
+		int n1;
+		char line[1024];
+		int i;
+		char **bits;
+
+		rval = fgets(line, 1023, fh);
+		if ( rval == NULL ) break;
+		chomp(line);
+
+		n1 = assplode(line, " \t", &bits, ASSPLODE_NONE);
+		if ( n1 < 3 ) {
+			for ( i=0; i<n1; i++ ) free(bits[i]);
+			free(bits);
+			continue;
+		}
+
+		if ( bits[1][0] != '=' ) {
+			for ( i=0; i<n1; i++ ) free(bits[i]);
+			free(bits);
+			continue;
+		}
+
+		if ( strcmp(bits[0], "beam/fluence") == 0 ) {
+			b->fluence = atof(bits[2]);
+		} else if ( strcmp(bits[0], "beam/radius") == 0 ) {
+			b->beam_radius = atof(bits[2]);
+		} else if ( strcmp(bits[0], "beam/photon_energy") == 0 ) {
+			b->photon_energy = atof(bits[2]);
+		} else if ( strcmp(bits[0], "beam/bandwidth") == 0 ) {
+			b->bandwidth = atof(bits[2]);
+		} else if ( strcmp(bits[0], "beam/divergence") == 0 ) {
+			b->divergence = atof(bits[2]);
+		} else if ( strcmp(bits[0], "detector/dqe") == 0 ) {
+			b->dqe = atof(bits[2]);
+		} else if ( strcmp(bits[0], "detector/adu_per_photon") == 0 ) {
+			b->adu_per_photon = atof(bits[2]);
+		} else {
+			ERROR("Unrecognised field '%s'\n", bits[0]);
+		}
+
+		for ( i=0; i<n1; i++ ) free(bits[i]);
+		free(bits);
+
+	} while ( rval != NULL );
+	fclose(fh);
+
+	reject = 0;
+
+	if ( b->fluence < 0.0 ) {
+		ERROR("Invalid or unspecified value for 'beam/fluence'.\n");
+		reject = 1;
+	}
+	if ( b->beam_radius < 0.0 ) {
+		ERROR("Invalid or unspecified value for 'beam/radius'.\n");
+		reject = 1;
+	}
+	if ( b->photon_energy < 0.0 ) {
+		ERROR("Invalid or unspecified value for"
+		      " 'beam/photon_energy'.\n");
+		reject = 1;
+	}
+	if ( b->bandwidth < 0.0 ) {
+		ERROR("Invalid or unspecified value for 'beam/bandwidth'.\n");
+		reject = 1;
+	}
+	if ( b->divergence < 0.0 ) {
+		ERROR("Invalid or unspecified value for 'beam/divergence'.\n");
+		reject = 1;
+	}
+	if ( b->dqe < 0.0 ) {
+		ERROR("Invalid or unspecified value for 'detector/dqe'.\n");
+		reject = 1;
+	}
+	if ( b->adu_per_photon < 0.0 ) {
+		ERROR("Invalid or unspecified value for"
+		      " 'detector/adu_per_photon'.\n");
+		reject = 1;
+	}
+
+	if ( reject ) {
+		ERROR("Please fix the above problems with the beam"
+		      " parameters file and try again.\n");
+		return NULL;
+	}
+
+	return b;
+}
diff --git a/libcrystfel/src/beam-parameters.h b/libcrystfel/src/beam-parameters.h
new file mode 100644
index 00000000..411ab449
--- /dev/null
+++ b/libcrystfel/src/beam-parameters.h
@@ -0,0 +1,40 @@
+/*
+ * beam-parameters.h
+ *
+ * Beam parameters
+ *
+ * (c) 2006-2010 Thomas White <taw@physics.org>
+ *
+ * Part of CrystFEL - crystallography with a FEL
+ *
+ */
+
+
+#ifndef BEAM_PARAMETERS_H
+#define BEAM_PARAMETERS_H
+
+#ifdef HAVE_CONFIG_H
+#include <config.h>
+#endif
+
+
+struct beam_params
+{
+	double fluence;        /* photons per pulse */
+	double beam_radius;    /* metres */
+	double photon_energy;  /* eV per photon */
+	double bandwidth;      /* FWHM(wavelength) over wavelength.
+	                        *  Note: current simulation code just uses
+	                        *        a rectangular distribution with this as
+	                        *        its (full) width. */
+	double divergence;     /* divergence (radians) */
+
+	double dqe;            /* Detector DQE (fraction) */
+	double adu_per_photon; /* Detector "gain" */
+};
+
+
+extern struct beam_params *get_beam_parameters(const char *filename);
+
+
+#endif	/* BEAM_PARAMETERS_H */
diff --git a/libcrystfel/src/cell.c b/libcrystfel/src/cell.c
new file mode 100644
index 00000000..c31697bc
--- /dev/null
+++ b/libcrystfel/src/cell.c
@@ -0,0 +1,1164 @@
+/*
+ * cell.c
+ *
+ * Unit Cell Calculations
+ *
+ * (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 <stdio.h>
+#include <string.h>
+#include <gsl/gsl_matrix.h>
+#include <gsl/gsl_blas.h>
+#include <gsl/gsl_linalg.h>
+
+#include "cell.h"
+#include "utils.h"
+#include "image.h"
+
+
+/**
+ * SECTION:unitcell
+ * @short_description: Unit cell
+ * @title: UnitCell
+ * @section_id:
+ * @see_also:
+ * @include: "cell.h"
+ * @Image:
+ *
+ * This structure represents a unit cell.
+ */
+
+
+
+/* Weighting factor of lengths relative to angles */
+#define LWEIGHT (10.0e-9)
+
+typedef enum {
+	CELL_REP_CRYST,
+	CELL_REP_CART,
+	CELL_REP_RECIP
+} CellRepresentation;
+
+struct _unitcell {
+
+	CellRepresentation rep;
+
+	/* Crystallographic representation */
+	double a;	/* m */
+	double b;	/* m */
+	double c;	/* m */
+	double alpha;	/* Radians */
+	double beta;	/* Radians */
+	double gamma;	/* Radians */
+
+	/* Cartesian representation */
+	double ax;	double bx;	double cx;
+	double ay;	double by;	double cy;
+	double az;	double bz;	double cz;
+
+	/* Cartesian representation of reciprocal axes */
+	double axs;	double bxs;	double cxs;
+	double ays;	double bys;	double cys;
+	double azs;	double bzs;	double czs;
+
+	char *pointgroup;
+	char *spacegroup;
+};
+
+
+/************************** Setters and Constructors **************************/
+
+
+/**
+ * cell_new:
+ *
+ * Create a new %UnitCell.
+ *
+ * Returns: the new unit cell, or NULL on failure.
+ *
+ */
+UnitCell *cell_new()
+{
+	UnitCell *cell;
+
+	cell = malloc(sizeof(UnitCell));
+	if ( cell == NULL ) return NULL;
+
+	cell->a = 1.0;
+	cell->b = 1.0;
+	cell->c = 1.0;
+	cell->alpha = M_PI_2;
+	cell->beta = M_PI_2;
+	cell->gamma = M_PI_2;
+
+	cell->rep = CELL_REP_CRYST;
+
+	cell->pointgroup = strdup("1");
+	cell->spacegroup = strdup("P 1");
+
+	return cell;
+}
+
+
+/**
+ * cell_free:
+ * @cell: A %UnitCell to free.
+ *
+ * Frees a %UnitCell, and all internal resources concerning that cell.
+ *
+ */
+void cell_free(UnitCell *cell)
+{
+	if ( cell == NULL ) return;
+	free(cell->pointgroup);
+	free(cell->spacegroup);
+	free(cell);
+}
+
+
+void cell_set_parameters(UnitCell *cell, double a, double b, double c,
+                         double alpha, double beta, double gamma)
+{
+	if ( cell == NULL ) return;
+
+	cell->a = a;
+	cell->b = b;
+	cell->c = c;
+	cell->alpha = alpha;
+	cell->beta = beta;
+	cell->gamma = gamma;
+
+	cell->rep = CELL_REP_CRYST;
+}
+
+
+void cell_set_cartesian(UnitCell *cell,
+			double ax, double ay, double az,
+			double bx, double by, double bz,
+			double cx, double cy, double cz)
+{
+	if ( cell == NULL ) return;
+
+	cell->ax = ax;  cell->ay = ay;  cell->az = az;
+	cell->bx = bx;  cell->by = by;  cell->bz = bz;
+	cell->cx = cx;  cell->cy = cy;  cell->cz = cz;
+
+	cell->rep = CELL_REP_CART;
+}
+
+
+void cell_set_cartesian_a(UnitCell *cell, double ax, double ay, double az)
+{
+	if ( cell == NULL ) return;
+	cell->ax = ax;  cell->ay = ay;  cell->az = az;
+	cell->rep = CELL_REP_CART;
+}
+
+
+void cell_set_cartesian_b(UnitCell *cell, double bx, double by, double bz)
+{
+	if ( cell == NULL ) return;
+	cell->bx = bx;  cell->by = by;  cell->bz = bz;
+	cell->rep = CELL_REP_CART;
+}
+
+
+void cell_set_cartesian_c(UnitCell *cell, double cx, double cy, double cz)
+{
+	if ( cell == NULL ) return;
+	cell->cx = cx;  cell->cy = cy;  cell->cz = cz;
+	cell->rep = CELL_REP_CART;
+}
+
+
+UnitCell *cell_new_from_parameters(double a, double b, double c,
+                                   double alpha, double beta, double gamma)
+{
+	UnitCell *cell;
+
+	cell = cell_new();
+	if ( cell == NULL ) return NULL;
+
+	cell_set_parameters(cell, a, b, c, alpha, beta, gamma);
+
+	return cell;
+}
+
+
+UnitCell *cell_new_from_reciprocal_axes(struct rvec as, struct rvec bs,
+                                        struct rvec cs)
+{
+	UnitCell *cell;
+
+	cell = cell_new();
+	if ( cell == NULL ) return NULL;
+
+	cell->axs = as.u;  cell->ays = as.v;  cell->azs = as.w;
+	cell->bxs = bs.u;  cell->bys = bs.v;  cell->bzs = bs.w;
+	cell->cxs = cs.u;  cell->cys = cs.v;  cell->czs = cs.w;
+
+	cell->rep = CELL_REP_RECIP;
+
+	return cell;
+}
+
+
+UnitCell *cell_new_from_direct_axes(struct rvec a, struct rvec b, struct rvec c)
+{
+	UnitCell *cell;
+
+	cell = cell_new();
+	if ( cell == NULL ) return NULL;
+
+	cell->ax = a.u;  cell->ay = a.v;  cell->az = a.w;
+	cell->bx = b.u;  cell->by = b.v;  cell->bz = b.w;
+	cell->cx = c.u;  cell->cy = c.v;  cell->cz = c.w;
+
+	cell->rep = CELL_REP_CART;
+
+	return cell;
+}
+
+
+UnitCell *cell_new_from_cell(UnitCell *orig)
+{
+	UnitCell *new;
+	double ax, ay, az, bx, by, bz, cx, cy, cz;
+
+	new = cell_new();
+
+	cell_get_cartesian(orig, &ax, &ay, &az, &bx, &by, &bz, &cx, &cy, &cz);
+	cell_set_cartesian(new, ax, ay, az, bx, by, bz, cx, cy, cz);
+	cell_set_pointgroup(new, orig->pointgroup);
+	cell_set_spacegroup(new, orig->spacegroup);
+
+	return new;
+}
+
+
+void cell_set_reciprocal(UnitCell *cell,
+                        double asx, double asy, double asz,
+                        double bsx, double bsy, double bsz,
+                        double csx, double csy, double csz)
+{
+	if ( cell == NULL ) return;
+
+	cell->axs = asx;  cell->ays = asy;  cell->azs = asz;
+	cell->bxs = bsx;  cell->bys = bsy;  cell->bzs = bsz;
+	cell->cxs = csx;  cell->cys = csy;  cell->czs = csz;
+
+	cell->rep = CELL_REP_RECIP;
+}
+
+
+void cell_set_spacegroup(UnitCell *cell, const char *sym)
+{
+	free(cell->spacegroup);
+	cell->spacegroup = strdup(sym);
+}
+
+
+void cell_set_pointgroup(UnitCell *cell, const char *sym)
+{
+	free(cell->pointgroup);
+	cell->pointgroup = strdup(sym);
+}
+
+
+/************************* Getter helper functions ****************************/
+
+static int cell_crystallographic_to_cartesian(UnitCell *cell,
+                                             double *ax, double *ay, double *az,
+                                             double *bx, double *by, double *bz,
+                                             double *cx, double *cy, double *cz)
+{
+	double tmp, V, cosalphastar, cstar;
+
+	/* Firstly: Get a in terms of x, y and z
+	 * +a (cryst) is defined to lie along +x (cart) */
+	*ax = cell->a;
+	*ay = 0.0;
+	*az = 0.0;
+
+	/* b in terms of x, y and z
+	 * b (cryst) is defined to lie in the xy (cart) plane */
+	*bx = cell->b*cos(cell->gamma);
+	*by = cell->b*sin(cell->gamma);
+	*bz = 0.0;
+
+	tmp = cos(cell->alpha)*cos(cell->alpha)
+		+ cos(cell->beta)*cos(cell->beta)
+		+ cos(cell->gamma)*cos(cell->gamma)
+		- 2.0*cos(cell->alpha)*cos(cell->beta)*cos(cell->gamma);
+	V = cell->a * cell->b * cell->c * sqrt(1.0 - tmp);
+
+	cosalphastar = cos(cell->beta)*cos(cell->gamma) - cos(cell->alpha);
+	cosalphastar /= sin(cell->beta)*sin(cell->gamma);
+
+	cstar = (cell->a * cell->b * sin(cell->gamma))/V;
+
+	/* c in terms of x, y and z */
+	*cx = cell->c*cos(cell->beta);
+	*cy = -cell->c*sin(cell->beta)*cosalphastar;
+	*cz = 1.0/cstar;
+
+	return 0;
+}
+
+
+/* Why yes, I do enjoy long argument lists...! */
+static int cell_invert(double ax, double ay, double az,
+                       double bx, double by, double bz,
+                       double cx, double cy, double cz,
+                       double *asx, double *asy, double *asz,
+                       double *bsx, double *bsy, double *bsz,
+                       double *csx, double *csy, double *csz)
+{
+	int s;
+	gsl_matrix *m;
+	gsl_matrix *inv;
+	gsl_permutation *perm;
+
+	m = gsl_matrix_alloc(3, 3);
+	if ( m == NULL ) {
+		ERROR("Couldn't allocate memory for matrix\n");
+		return 1;
+	}
+	gsl_matrix_set(m, 0, 0, ax);
+	gsl_matrix_set(m, 0, 1, bx);
+	gsl_matrix_set(m, 0, 2, cx);
+	gsl_matrix_set(m, 1, 0, ay);
+	gsl_matrix_set(m, 1, 1, by);
+	gsl_matrix_set(m, 1, 2, cy);
+	gsl_matrix_set(m, 2, 0, az);
+	gsl_matrix_set(m, 2, 1, bz);
+	gsl_matrix_set(m, 2, 2, cz);
+
+	/* Invert */
+	perm = gsl_permutation_alloc(m->size1);
+	if ( perm == NULL ) {
+		ERROR("Couldn't allocate permutation\n");
+		gsl_matrix_free(m);
+		return 1;
+	}
+	inv = gsl_matrix_alloc(m->size1, m->size2);
+	if ( inv == NULL ) {
+		ERROR("Couldn't allocate inverse\n");
+		gsl_matrix_free(m);
+		gsl_permutation_free(perm);
+		return 1;
+	}
+	if ( gsl_linalg_LU_decomp(m, perm, &s) ) {
+		ERROR("Couldn't decompose matrix\n");
+		gsl_matrix_free(m);
+		gsl_permutation_free(perm);
+		return 1;
+	}
+	if ( gsl_linalg_LU_invert(m, perm, inv)  ) {
+		ERROR("Couldn't invert matrix\n");
+		gsl_matrix_free(m);
+		gsl_permutation_free(perm);
+		return 1;
+	}
+	gsl_permutation_free(perm);
+	gsl_matrix_free(m);
+
+	/* Transpose */
+	gsl_matrix_transpose(inv);
+
+	*asx = gsl_matrix_get(inv, 0, 0);
+	*bsx = gsl_matrix_get(inv, 0, 1);
+	*csx = gsl_matrix_get(inv, 0, 2);
+	*asy = gsl_matrix_get(inv, 1, 0);
+	*bsy = gsl_matrix_get(inv, 1, 1);
+	*csy = gsl_matrix_get(inv, 1, 2);
+	*asz = gsl_matrix_get(inv, 2, 0);
+	*bsz = gsl_matrix_get(inv, 2, 1);
+	*csz = gsl_matrix_get(inv, 2, 2);
+
+	gsl_matrix_free(inv);
+
+	return 0;
+}
+
+
+/********************************** Getters ***********************************/
+
+int cell_get_parameters(UnitCell *cell, double *a, double *b, double *c,
+                        double *alpha, double *beta, double *gamma)
+{
+	double ax, ay, az, bx, by, bz, cx, cy, cz;
+
+	if ( cell == NULL ) return 1;
+
+	switch ( cell->rep ) {
+
+	case CELL_REP_CRYST:
+		/* Direct response */
+		*a = cell->a;
+		*b = cell->b;
+		*c = cell->c;
+		*alpha = cell->alpha;
+		*beta = cell->beta;
+		*gamma = cell->gamma;
+		return 0;
+
+	case CELL_REP_CART:
+		/* Convert cartesian -> crystallographic */
+		*a = modulus(cell->ax, cell->ay, cell->az);
+		*b = modulus(cell->bx, cell->by, cell->bz);
+		*c = modulus(cell->cx, cell->cy, cell->cz);
+
+		*alpha = angle_between(cell->bx, cell->by, cell->bz,
+		                       cell->cx, cell->cy, cell->cz);
+		*beta = angle_between(cell->ax, cell->ay, cell->az,
+		                      cell->cx, cell->cy, cell->cz);
+		*gamma = angle_between(cell->ax, cell->ay, cell->az,
+		                       cell->bx, cell->by, cell->bz);
+		return 0;
+
+	case CELL_REP_RECIP:
+		/* Convert reciprocal -> crystallographic.
+                 * Start by converting reciprocal -> cartesian */
+		cell_invert(cell->axs, cell->ays, cell->azs,
+		            cell->bxs, cell->bys, cell->bzs,
+		            cell->cxs, cell->cys, cell->czs,
+		            &ax, &ay, &az, &bx, &by, &bz, &cx, &cy, &cz);
+
+		/* Now convert cartesian -> crystallographic */
+		*a = modulus(ax, ay, az);
+		*b = modulus(bx, by, bz);
+		*c = modulus(cx, cy, cz);
+
+		*alpha = angle_between(bx, by, bz, cx, cy, cz);
+		*beta = angle_between(ax, ay, az, cx, cy, cz);
+		*gamma = angle_between(ax, ay, az, bx, by, bz);
+		return 0;
+	}
+
+	return 1;
+}
+
+
+int cell_get_cartesian(UnitCell *cell,
+                       double *ax, double *ay, double *az,
+                       double *bx, double *by, double *bz,
+                       double *cx, double *cy, double *cz)
+{
+	if ( cell == NULL ) return 1;
+
+	switch ( cell->rep ) {
+
+	case CELL_REP_CRYST:
+		/* Convert crystallographic -> cartesian. */
+		return cell_crystallographic_to_cartesian(cell,
+		                                          ax, ay, az,
+		                                          bx, by, bz,
+		                                          cx, cy, cz);
+
+	case CELL_REP_CART:
+		/* Direct response */
+		*ax = cell->ax;  *ay = cell->ay;  *az = cell->az;
+		*bx = cell->bx;  *by = cell->by;  *bz = cell->bz;
+		*cx = cell->cx;  *cy = cell->cy;  *cz = cell->cz;
+		return 0;
+
+	case CELL_REP_RECIP:
+		/* Convert reciprocal -> cartesian */
+		return cell_invert(cell->axs, cell->ays, cell->azs,
+		                   cell->bxs, cell->bys, cell->bzs,
+		                   cell->cxs, cell->cys, cell->czs,
+		                   ax, ay, az, bx, by, bz, cx, cy, cz);
+
+	}
+
+	return 1;
+}
+
+
+int cell_get_reciprocal(UnitCell *cell,
+                        double *asx, double *asy, double *asz,
+                        double *bsx, double *bsy, double *bsz,
+                        double *csx, double *csy, double *csz)
+{
+	int r;
+	double ax, ay, az, bx, by, bz, cx, cy, cz;
+	if ( cell == NULL ) return 1;
+
+	switch ( cell->rep ) {
+
+	case CELL_REP_CRYST:
+		/* Convert crystallographic -> reciprocal */
+		r = cell_crystallographic_to_cartesian(cell,
+		                                       &ax, &ay, &az,
+		                                       &bx, &by, &bz,
+		                                       &cx, &cy, &cz);
+		if ( r ) return r;
+		return cell_invert(ax, ay, az,bx, by, bz, cx, cy, cz,
+		                   asx, asy, asz, bsx, bsy, bsz, csx, csy, csz);
+
+	case CELL_REP_CART:
+		/* Convert cartesian -> reciprocal */
+		cell_invert(cell->ax, cell->ay, cell->az,
+		            cell->bx, cell->by, cell->bz,
+		            cell->cx, cell->cy, cell->cz,
+		            asx, asy, asz, bsx, bsy, bsz, csx, csy, csz);
+		return 0;
+
+	case CELL_REP_RECIP:
+		/* Direct response */
+		*asx = cell->axs;  *asy = cell->ays;  *asz = cell->azs;
+		*bsx = cell->bxs;  *bsy = cell->bys;  *bsz = cell->bzs;
+		*csx = cell->cxs;  *csy = cell->cys;  *csz = cell->czs;
+		return 0;
+
+	}
+
+	return 1;
+}
+
+
+const char *cell_get_pointgroup(UnitCell *cell)
+{
+	return cell->pointgroup;
+}
+
+
+const char *cell_get_spacegroup(UnitCell *cell)
+{
+	return cell->spacegroup;
+}
+
+
+
+
+
+/********************************* Utilities **********************************/
+
+static const char *cell_rep(UnitCell *cell)
+{
+	switch ( cell->rep ) {
+	case CELL_REP_CRYST:
+		return "crystallographic, direct space";
+	case CELL_REP_CART:
+		return "cartesian, direct space";
+	case CELL_REP_RECIP:
+		return "cartesian, reciprocal space";
+	}
+
+	return "unknown";
+}
+
+
+/**
+ * cell_rotate:
+ * @in: A %UnitCell to rotate
+ * @quat: A %quaternion
+ *
+ * Rotate a %UnitCell using a %quaternion.
+ *
+ * Returns: a newly allocated rotated copy of @in.
+ *
+ */
+UnitCell *cell_rotate(UnitCell *in, struct quaternion quat)
+{
+	struct rvec a, b, c;
+	struct rvec an, bn, cn;
+	UnitCell *out = cell_new_from_cell(in);
+
+	cell_get_cartesian(in, &a.u, &a.v, &a.w,
+	                       &b.u, &b.v, &b.w,
+	                       &c.u, &c.v, &c.w);
+
+	an = quat_rot(a, quat);
+	bn = quat_rot(b, quat);
+	cn = quat_rot(c, quat);
+
+	cell_set_cartesian(out, an.u, an.v, an.w,
+	                        bn.u, bn.v, bn.w,
+	                        cn.u, cn.v, cn.w);
+
+	return out;
+}
+
+
+void cell_print(UnitCell *cell)
+{
+	double asx, asy, asz;
+	double bsx, bsy, bsz;
+	double csx, csy, csz;
+	double a, b, c, alpha, beta, gamma;
+	double ax, ay, az, bx, by, bz, cx, cy, cz;
+
+	cell_get_parameters(cell, &a, &b, &c, &alpha, &beta, &gamma);
+
+	STATUS("  a     b     c         alpha   beta  gamma\n");
+	STATUS("%5.2f %5.2f %5.2f nm    %6.2f %6.2f %6.2f deg\n",
+	       a*1e9, b*1e9, c*1e9,
+	       rad2deg(alpha), rad2deg(beta), rad2deg(gamma));
+
+	cell_get_cartesian(cell, &ax, &ay, &az, &bx, &by, &bz, &cx, &cy, &cz);
+
+	STATUS("a = %10.3e %10.3e %10.3e m\n", ax, ay, az);
+	STATUS("b = %10.3e %10.3e %10.3e m\n", bx, by, bz);
+	STATUS("c = %10.3e %10.3e %10.3e m\n", cx, cy, cz);
+
+	cell_get_reciprocal(cell, &asx, &asy, &asz,
+	                          &bsx, &bsy, &bsz,
+	                          &csx, &csy, &csz);
+
+	STATUS("astar = %10.3e %10.3e %10.3e m^-1 (modulus = %10.3e m^-1)\n",
+	                             asx, asy, asz, modulus(asx, asy, asz));
+	STATUS("bstar = %10.3e %10.3e %10.3e m^-1 (modulus = %10.3e m^-1)\n",
+	                             bsx, bsy, bsz, modulus(bsx, bsy, bsz));
+	STATUS("cstar = %10.3e %10.3e %10.3e m^-1 (modulus = %10.3e m^-1)\n",
+	                             csx, csy, csz, modulus(csx, csy, csz));
+
+	STATUS("Point group: %s\n", cell_get_pointgroup(cell));
+	STATUS("Cell representation is %s.\n", cell_rep(cell));
+}
+
+
+#define MAX_CAND (1024)
+
+static int right_handed(struct rvec a, struct rvec b, struct rvec c)
+{
+	struct rvec aCb;
+	double aCb_dot_c;
+
+	/* "a" cross "b" */
+	aCb.u = a.v*b.w - a.w*b.v;
+	aCb.v = - (a.u*b.w - a.w*b.u);
+	aCb.w = a.u*b.v - a.v*b.u;
+
+	/* "a cross b" dot "c" */
+	aCb_dot_c = aCb.u*c.u + aCb.v*c.v + aCb.w*c.w;
+
+	if ( aCb_dot_c > 0.0 ) return 1;
+	return 0;
+}
+
+
+struct cvec {
+	struct rvec vec;
+	float na;
+	float nb;
+	float nc;
+	float fom;
+};
+
+
+static int same_vector(struct cvec a, struct cvec b)
+{
+	if ( a.na != b.na ) return 0;
+	if ( a.nb != b.nb ) return 0;
+	if ( a.nc != b.nc ) return 0;
+	return 1;
+}
+
+
+/* Attempt to make 'cell' fit into 'template' somehow */
+UnitCell *match_cell(UnitCell *cell, UnitCell *template, int verbose,
+                     int reduce)
+{
+	signed int n1l, n2l, n3l;
+	double asx, asy, asz;
+	double bsx, bsy, bsz;
+	double csx, csy, csz;
+	int i, j;
+	double lengths[3];
+	double angles[3];
+	struct cvec *cand[3];
+	UnitCell *new_cell = NULL;
+	float best_fom = +999999999.9; /* Large number.. */
+	int ncand[3] = {0,0,0};
+	signed int ilow, ihigh;
+	float ltl = 5.0;     /* percent */
+	float angtol = deg2rad(1.5);
+
+	if ( cell_get_reciprocal(template, &asx, &asy, &asz,
+	                         &bsx, &bsy, &bsz,
+	                         &csx, &csy, &csz) ) {
+		ERROR("Couldn't get reciprocal cell for template.\n");
+		return NULL;
+	}
+
+	lengths[0] = modulus(asx, asy, asz);
+	lengths[1] = modulus(bsx, bsy, bsz);
+	lengths[2] = modulus(csx, csy, csz);
+
+	angles[0] = angle_between(bsx, bsy, bsz, csx, csy, csz);
+	angles[1] = angle_between(asx, asy, asz, csx, csy, csz);
+	angles[2] = angle_between(asx, asy, asz, bsx, bsy, bsz);
+
+	cand[0] = malloc(MAX_CAND*sizeof(struct cvec));
+	cand[1] = malloc(MAX_CAND*sizeof(struct cvec));
+	cand[2] = malloc(MAX_CAND*sizeof(struct cvec));
+
+	if ( cell_get_reciprocal(cell, &asx, &asy, &asz,
+	                         &bsx, &bsy, &bsz,
+	                         &csx, &csy, &csz) ) {
+		ERROR("Couldn't get reciprocal cell.\n");
+		return NULL;
+	}
+
+	if ( reduce ) {
+		ilow = -2;  ihigh = 4;
+	} else {
+		ilow = 0;  ihigh = 1;
+	}
+
+	/* Negative values mean 1/n, positive means n, zero means zero */
+	for ( n1l=ilow; n1l<=ihigh; n1l++ ) {
+	for ( n2l=ilow; n2l<=ihigh; n2l++ ) {
+	for ( n3l=ilow; n3l<=ihigh; n3l++ ) {
+
+		float n1, n2, n3;
+		signed int b1, b2, b3;
+
+		n1 = (n1l>=0) ? (n1l) : (1.0/n1l);
+		n2 = (n2l>=0) ? (n2l) : (1.0/n2l);
+		n3 = (n3l>=0) ? (n3l) : (1.0/n3l);
+
+		if ( !reduce ) {
+			if ( n1l + n2l + n3l > 1 ) continue;
+		}
+
+		/* 'bit' values can be +1 or -1 */
+		for ( b1=-1; b1<=1; b1+=2 ) {
+		for ( b2=-1; b2<=1; b2+=2 ) {
+		for ( b3=-1; b3<=1; b3+=2 ) {
+
+			double tx, ty, tz;
+			double tlen;
+			int i;
+
+			n1 *= b1;  n2 *= b2;  n3 *= b3;
+
+			tx = n1*asx + n2*bsx + n3*csx;
+			ty = n1*asy + n2*bsy + n3*csy;
+			tz = n1*asz + n2*bsz + n3*csz;
+			tlen = modulus(tx, ty, tz);
+
+			/* Test modulus for agreement with moduli of template */
+			for ( i=0; i<3; i++ ) {
+
+				if ( !within_tolerance(lengths[i], tlen, ltl) )
+					continue;
+
+				cand[i][ncand[i]].vec.u = tx;
+				cand[i][ncand[i]].vec.v = ty;
+				cand[i][ncand[i]].vec.w = tz;
+				cand[i][ncand[i]].na = n1;
+				cand[i][ncand[i]].nb = n2;
+				cand[i][ncand[i]].nc = n3;
+				cand[i][ncand[i]].fom = fabs(lengths[i] - tlen);
+				if ( ncand[i] == MAX_CAND ) {
+					ERROR("Too many candidates\n");
+				} else {
+					ncand[i]++;
+				}
+			}
+
+		}
+		}
+		}
+
+	}
+	}
+	}
+
+	if ( verbose ) {
+		STATUS("Candidates: %i %i %i\n", ncand[0], ncand[1], ncand[2]);
+	}
+
+	for ( i=0; i<ncand[0]; i++ ) {
+	for ( j=0; j<ncand[1]; j++ ) {
+
+		double ang;
+		int k;
+		float fom1;
+
+		if ( same_vector(cand[0][i], cand[1][j]) ) continue;
+
+		/* Measure the angle between the ith candidate for axis 0
+		 * and the jth candidate for axis 1 */
+		ang = angle_between(cand[0][i].vec.u, cand[0][i].vec.v,
+		                    cand[0][i].vec.w, cand[1][j].vec.u,
+		                    cand[1][j].vec.v, cand[1][j].vec.w);
+
+		/* Angle between axes 0 and 1 should be angle 2 */
+		if ( fabs(ang - angles[2]) > angtol ) continue;
+
+		fom1 = fabs(ang - angles[2]);
+
+		for ( k=0; k<ncand[2]; k++ ) {
+
+			float fom2, fom3;
+
+			if ( same_vector(cand[1][j], cand[2][k]) ) continue;
+
+			/* Measure the angle between the current candidate for
+			 * axis 0 and the kth candidate for axis 2 */
+			ang = angle_between(cand[0][i].vec.u, cand[0][i].vec.v,
+			                    cand[0][i].vec.w, cand[2][k].vec.u,
+			                    cand[2][k].vec.v, cand[2][k].vec.w);
+
+			/* ... it should be angle 1 ... */
+			if ( fabs(ang - angles[1]) > angtol ) continue;
+
+			fom2 = fom1 + fabs(ang - angles[1]);
+
+			/* Finally, the angle between the current candidate for
+			 * axis 1 and the kth candidate for axis 2 */
+			ang = angle_between(cand[1][j].vec.u, cand[1][j].vec.v,
+			                    cand[1][j].vec.w, cand[2][k].vec.u,
+			                    cand[2][k].vec.v, cand[2][k].vec.w);
+
+			/* ... it should be angle 0 ... */
+			if ( fabs(ang - angles[0]) > angtol ) continue;
+
+			/* Unit cell must be right-handed */
+			if ( !right_handed(cand[0][i].vec, cand[1][j].vec,
+			                   cand[2][k].vec) ) continue;
+
+			fom3 = fom2 + fabs(ang - angles[0]);
+			fom3 += LWEIGHT * (cand[0][i].fom + cand[1][j].fom
+			                   + cand[2][k].fom);
+
+			if ( fom3 < best_fom ) {
+				if ( new_cell != NULL ) free(new_cell);
+				new_cell = cell_new_from_reciprocal_axes(
+				                 cand[0][i].vec, cand[1][j].vec,
+			                         cand[2][k].vec);
+				best_fom = fom3;
+			}
+
+		}
+
+	}
+	}
+
+	free(cand[0]);
+	free(cand[1]);
+	free(cand[2]);
+
+	return new_cell;
+}
+
+
+
+UnitCell *match_cell_ab(UnitCell *cell, UnitCell *template)
+{
+	double ax, ay, az;
+	double bx, by, bz;
+	double cx, cy, cz;
+	int i;
+	double lengths[3];
+	int used[3];
+	struct rvec real_a, real_b, real_c;
+	struct rvec params[3];
+	double alen, blen, clen;
+	float ltl = 5.0;     /* percent */
+	int have_real_a;
+	int have_real_b;
+	int have_real_c;
+
+	/* Get the lengths to match */
+	if ( cell_get_cartesian(template, &ax, &ay, &az,
+	                                  &bx, &by, &bz,
+	                                  &cx, &cy, &cz) )
+	{
+		ERROR("Couldn't get cell for template.\n");
+		return NULL;
+	}
+	alen = modulus(ax, ay, az);
+	blen = modulus(bx, by, bz);
+	clen = modulus(cx, cy, cz);
+
+	/* Get the lengths from the cell and turn them into anonymous vectors */
+	if ( cell_get_cartesian(cell, &ax, &ay, &az,
+	                              &bx, &by, &bz,
+	                              &cx, &cy, &cz) )
+	{
+		ERROR("Couldn't get cell.\n");
+		return NULL;
+	}
+	lengths[0] = modulus(ax, ay, az);
+	lengths[1] = modulus(bx, by, bz);
+	lengths[2] = modulus(cx, cy, cz);
+	used[0] = 0;  used[1] = 0;  used[2] = 0;
+	params[0].u = ax;  params[0].v = ay;  params[0].w = az;
+	params[1].u = bx;  params[1].v = by;  params[1].w = bz;
+	params[2].u = cx;  params[2].v = cy;  params[2].w = cz;
+
+	real_a.u = 0.0;  real_a.v = 0.0;  real_a.w = 0.0;
+	real_b.u = 0.0;  real_b.v = 0.0;  real_b.w = 0.0;
+	real_c.u = 0.0;  real_c.v = 0.0;  real_c.w = 0.0;
+
+	/* Check each vector against a and b */
+	have_real_a = 0;
+	have_real_b = 0;
+	for ( i=0; i<3; i++ ) {
+		if ( within_tolerance(lengths[i], alen, ltl)
+		     && !used[i] && !have_real_a )
+		{
+			used[i] = 1;
+			memcpy(&real_a, &params[i], sizeof(struct rvec));
+			have_real_a = 1;
+		}
+		if ( within_tolerance(lengths[i], blen, ltl)
+		     && !used[i] && !have_real_b )
+		{
+			used[i] = 1;
+			memcpy(&real_b, &params[i], sizeof(struct rvec));
+			have_real_b = 1;
+		}
+	}
+
+	/* Have we matched both a and b? */
+	if ( !(have_real_a && have_real_b) ) return NULL;
+
+	/* "c" is "the other one" */
+	have_real_c = 0;
+	for ( i=0; i<3; i++ ) {
+		if ( !used[i] ) {
+			memcpy(&real_c, &params[i], sizeof(struct rvec));
+			have_real_c = 1;
+		}
+	}
+
+	if ( !have_real_c ) {
+		ERROR("Huh?  Couldn't find the third vector.\n");
+		ERROR("Matches: %i %i %i\n", used[0], used[1], used[2]);
+		return NULL;
+	}
+
+	/* Flip c if not right-handed */
+	if ( !right_handed(real_a, real_b, real_c) ) {
+		real_c.u = -real_c.u;
+		real_c.v = -real_c.v;
+		real_c.w = -real_c.w;
+	}
+
+	return cell_new_from_direct_axes(real_a, real_b, real_c);
+}
+
+
+/* Return sin(theta)/lambda = 1/2d.  Multiply by two if you want 1/d */
+double resolution(UnitCell *cell, signed int h, signed int k, signed int l)
+{
+	double a, b, c, alpha, beta, gamma;
+
+	cell_get_parameters(cell, &a, &b, &c, &alpha, &beta, &gamma);
+
+	const double Vsq = a*a*b*b*c*c*(1 - cos(alpha)*cos(alpha)
+	                                  - cos(beta)*cos(beta)
+	                                  - cos(gamma)*cos(gamma)
+				          + 2*cos(alpha)*cos(beta)*cos(gamma) );
+
+	const double S11 = b*b*c*c*sin(alpha)*sin(alpha);
+	const double S22 = a*a*c*c*sin(beta)*sin(beta);
+	const double S33 = a*a*b*b*sin(gamma)*sin(gamma);
+	const double S12 = a*b*c*c*(cos(alpha)*cos(beta) - cos(gamma));
+	const double S23 = a*a*b*c*(cos(beta)*cos(gamma) - cos(alpha));
+	const double S13 = a*b*b*c*(cos(gamma)*cos(alpha) - cos(beta));
+
+	const double brackets = S11*h*h + S22*k*k + S33*l*l
+	                         + 2*S12*h*k + 2*S23*k*l + 2*S13*h*l;
+	const double oneoverdsq = brackets / Vsq;
+	const double oneoverd = sqrt(oneoverdsq);
+
+	return oneoverd / 2;
+}
+
+
+static void cell_set_pointgroup_from_pdb(UnitCell *cell, const char *sym)
+{
+	char *new = NULL;
+
+	if ( strcmp(sym, "P 1") == 0 ) new = "1";
+	if ( strcmp(sym, "P 63") == 0 ) new = "6";
+	if ( strcmp(sym, "P 21 21 21") == 0 ) new = "222";
+	if ( strcmp(sym, "P 2 2 2") == 0 ) new = "222";
+	if ( strcmp(sym, "P 43 21 2") == 0 ) new = "422";
+
+	if ( new != NULL ) {
+		cell_set_pointgroup(cell, new);
+	} else {
+		ERROR("Can't determine point group for '%s'\n", sym);
+	}
+}
+
+
+UnitCell *load_cell_from_pdb(const char *filename)
+{
+	FILE *fh;
+	char *rval;
+	UnitCell *cell = NULL;
+
+	fh = fopen(filename, "r");
+	if ( fh == NULL ) {
+		ERROR("Couldn't open '%s'\n", filename);
+		return NULL;
+	}
+
+	do {
+
+		char line[1024];
+
+		rval = fgets(line, 1023, fh);
+
+		if ( strncmp(line, "CRYST1", 6) == 0 ) {
+
+			float a, b, c, al, be, ga;
+			char as[10], bs[10], cs[10];
+			char als[8], bes[8], gas[8];
+			char *sym;
+			int r;
+
+			memcpy(as, line+6, 9);    as[9] = '\0';
+			memcpy(bs, line+15, 9);   bs[9] = '\0';
+			memcpy(cs, line+24, 9);   cs[9] = '\0';
+			memcpy(als, line+33, 7);  als[7] = '\0';
+			memcpy(bes, line+40, 7);  bes[7] = '\0';
+			memcpy(gas, line+47, 7);  gas[7] = '\0';
+
+			r = sscanf(as, "%f", &a);
+			r += sscanf(bs, "%f", &b);
+			r += sscanf(cs, "%f", &c);
+			r += sscanf(als, "%f", &al);
+			r += sscanf(bes, "%f", &be);
+			r += sscanf(gas, "%f", &ga);
+
+			if ( r != 6 ) {
+				STATUS("Couldn't understand CRYST1 line.\n");
+				continue;
+			}
+
+			cell = cell_new_from_parameters(a*1e-10,
+			                                b*1e-10, c*1e-10,
+	                                                deg2rad(al),
+	                                                deg2rad(be),
+	                                                deg2rad(ga));
+
+			if ( strlen(line) > 65 ) {
+				sym = strndup(line+55, 10);
+				notrail(sym);
+				cell_set_pointgroup_from_pdb(cell, sym);
+				cell_set_spacegroup(cell, sym);
+				free(sym);
+			} else {
+				cell_set_pointgroup_from_pdb(cell, "P 1");
+				cell_set_spacegroup(cell, "P 1");
+				ERROR("CRYST1 line without space group.\n");
+			}
+
+			break;  /* Done */
+		}
+
+	} while ( rval != NULL );
+
+	fclose(fh);
+
+	return cell;
+}
+
+
+#ifdef GSL_FUDGE
+/* Force the linker to bring in CBLAS to make GSL happy */
+void cell_fudge_gslcblas()
+{
+        STATUS("%p\n", cblas_sgemm);
+}
+#endif
+
+
+UnitCell *rotate_cell(UnitCell *in, double omega, double phi, double rot)
+{
+	UnitCell *out;
+	double asx, asy, asz;
+	double bsx, bsy, bsz;
+	double csx, csy, csz;
+	double xnew, ynew, znew;
+
+	cell_get_reciprocal(in, &asx, &asy, &asz, &bsx, &bsy,
+	                        &bsz, &csx, &csy, &csz);
+
+	/* Rotate by "omega" about +z (parallel to c* and c unless triclinic) */
+	xnew = asx*cos(omega) + asy*sin(omega);
+	ynew = -asx*sin(omega) + asy*cos(omega);
+	znew = asz;
+	asx = xnew;  asy = ynew;  asz = znew;
+	xnew = bsx*cos(omega) + bsy*sin(omega);
+	ynew = -bsx*sin(omega) + bsy*cos(omega);
+	znew = bsz;
+	bsx = xnew;  bsy = ynew;  bsz = znew;
+	xnew = csx*cos(omega) + csy*sin(omega);
+	ynew = -csx*sin(omega) + csy*cos(omega);
+	znew = csz;
+	csx = xnew;  csy = ynew;  csz = znew;
+
+	/* Rotate by "phi" about +x (not parallel to anything specific) */
+	xnew = asx;
+	ynew = asy*cos(phi) + asz*sin(phi);
+	znew = -asy*sin(phi) + asz*cos(phi);
+	asx = xnew;  asy = ynew;  asz = znew;
+	xnew = bsx;
+	ynew = bsy*cos(phi) + bsz*sin(phi);
+	znew = -bsy*sin(phi) + bsz*cos(phi);
+	bsx = xnew;  bsy = ynew;  bsz = znew;
+	xnew = csx;
+	ynew = csy*cos(phi) + csz*sin(phi);
+	znew = -csy*sin(phi) + csz*cos(phi);
+	csx = xnew;  csy = ynew;  csz = znew;
+
+	/* Rotate by "rot" about the new +z (in-plane rotation) */
+	xnew = asx*cos(rot) + asy*sin(rot);
+	ynew = -asx*sin(rot) + asy*cos(rot);
+	znew = asz;
+	asx = xnew;  asy = ynew;  asz = znew;
+	xnew = bsx*cos(rot) + bsy*sin(rot);
+	ynew = -bsx*sin(rot) + bsy*cos(rot);
+	znew = bsz;
+	bsx = xnew;  bsy = ynew;  bsz = znew;
+	xnew = csx*cos(rot) + csy*sin(rot);
+	ynew = -csx*sin(rot) + csy*cos(rot);
+	znew = csz;
+	csx = xnew;  csy = ynew;  csz = znew;
+
+	out = cell_new_from_cell(in);
+	cell_set_reciprocal(out, asx, asy, asz, bsx, bsy, bsz, csx, csy, csz);
+
+	return out;
+}
+
+
+int cell_is_sensible(UnitCell *cell)
+{
+	double a, b, c, al, be, ga;
+
+	cell_get_parameters(cell, &a, &b, &c, &al, &be, &ga);
+	if (   al + be + ga >= 2.0*M_PI ) return 0;
+	if (   al + be - ga >= 2.0*M_PI ) return 0;
+	if (   al - be + ga >= 2.0*M_PI ) return 0;
+	if ( - al + be + ga >= 2.0*M_PI ) return 0;
+	if (   al + be + ga <= 0.0 ) return 0;
+	if (   al + be - ga <= 0.0 ) return 0;
+	if (   al - be + ga <= 0.0 ) return 0;
+	if ( - al + be + ga <= 0.0 ) return 0;
+	if ( isnan(al) ) return 0;
+	if ( isnan(be) ) return 0;
+	if ( isnan(ga) ) return 0;
+	return 1;
+}
diff --git a/libcrystfel/src/cell.h b/libcrystfel/src/cell.h
new file mode 100644
index 00000000..b5d31fc6
--- /dev/null
+++ b/libcrystfel/src/cell.h
@@ -0,0 +1,107 @@
+/*
+ * cell.h
+ *
+ * Unit Cell Calculations
+ *
+ * (c) 2006-2010 Thomas White <taw@physics.org>
+ *
+ * Part of CrystFEL - crystallography with a FEL
+ *
+ */
+
+#ifndef CELL_H
+#define CELL_H
+
+#ifdef HAVE_CONFIG_H
+#include <config.h>
+#endif
+
+#include "utils.h"
+
+/* A 3D vector in reciprocal space */
+struct rvec
+{
+	double   u;
+	double   v;
+	double   w;
+};
+
+
+/**
+ * UnitCell:
+ *
+ * This data structure is opaque.  You must use the available accessor functions
+ * to read and write its contents.
+ **/
+typedef struct _unitcell UnitCell;
+
+extern UnitCell *cell_new(void);
+extern UnitCell *cell_new_from_cell(UnitCell *orig);
+extern void cell_free(UnitCell *cell);
+
+/* Lengths in m, angles in radians */
+extern UnitCell *cell_new_from_parameters(double a, double b, double c,
+				double alpha, double beta, double gamma);
+
+extern UnitCell *cell_new_from_reciprocal_axes(struct rvec as, struct rvec bs,
+                                               struct rvec cs);
+
+extern UnitCell *cell_new_from_direct_axes(struct rvec as, struct rvec bs,
+                                           struct rvec cs);
+
+extern void cell_set_cartesian(UnitCell *cell,
+                               double ax, double ay, double az,
+                               double bx, double by, double bz,
+                               double cx, double cy, double cz);
+
+extern void cell_set_parameters(UnitCell *cell, double a, double b, double c,
+				double alpha, double beta, double gamma);
+
+extern void cell_set_cartesian_a(UnitCell *cell, double ax, double ay, double az);
+extern void cell_set_cartesian_b(UnitCell *cell, double bx, double by, double bz);
+extern void cell_set_cartesian_c(UnitCell *cell, double cx, double cy, double cz);
+extern void cell_set_spacegroup(UnitCell *cell, const char *sym);
+extern void cell_set_pointgroup(UnitCell *cell, const char *sym);
+
+
+extern int cell_get_parameters(UnitCell *cell, double *a, double *b, double *c,
+                               double *alpha, double *beta, double *gamma);
+
+extern int cell_get_cartesian(UnitCell *cell,
+                              double *ax, double *ay, double *az,
+                              double *bx, double *by, double *bz,
+                              double *cx, double *cy, double *cz);
+
+extern int cell_get_reciprocal(UnitCell *cell,
+                               double *asx, double *asy, double *asz,
+                               double *bsx, double *bsy, double *bsz,
+                               double *csx, double *csy, double *csz);
+
+extern void cell_set_reciprocal(UnitCell *cell,
+                                double asx, double asy, double asz,
+                                double bsx, double bsy, double bsz,
+                                double csx, double csy, double csz);
+
+extern const char *cell_get_pointgroup(UnitCell *cell);
+
+extern const char *cell_get_spacegroup(UnitCell *cell);
+
+extern double resolution(UnitCell *cell,
+                         signed int h, signed int k, signed int l);
+
+extern UnitCell *cell_rotate(UnitCell *in, struct quaternion quat);
+extern UnitCell *rotate_cell(UnitCell *in, double omega, double phi,
+                             double rot);
+
+extern void cell_print(UnitCell *cell);
+
+extern UnitCell *match_cell(UnitCell *cell, UnitCell *template, int verbose,
+                            int reduce);
+
+extern UnitCell *match_cell_ab(UnitCell *cell, UnitCell *template);
+
+extern UnitCell *load_cell_from_pdb(const char *filename);
+
+extern int cell_is_sensible(UnitCell *cell);
+
+#endif	/* CELL_H */
diff --git a/libcrystfel/src/detector.c b/libcrystfel/src/detector.c
new file mode 100644
index 00000000..54bae1d7
--- /dev/null
+++ b/libcrystfel/src/detector.c
@@ -0,0 +1,1233 @@
+/*
+ * detector.c
+ *
+ * Detector properties
+ *
+ * (c) 2006-2011 Thomas White <taw@physics.org>
+ * (c) 2011 Rick Kirian <rkirian@asu.edu>
+ *
+ * Part of CrystFEL - crystallography with a FEL
+ *
+ */
+
+
+#include <stdlib.h>
+#include <math.h>
+#include <stdio.h>
+#include <string.h>
+#include <assert.h>
+#include <ctype.h>
+
+#include "image.h"
+#include "utils.h"
+#include "detector.h"
+#include "beam-parameters.h"
+#include "hdf5-file.h"
+
+
+/**
+ * SECTION:detector
+ * @short_description: Detector geometry
+ * @title: Detector
+ * @section_id:
+ * @see_also:
+ * @include: "detector.h"
+ * @Image:
+ *
+ * This structure represents the detector geometry
+ */
+
+
+static int atob(const char *a)
+{
+	if ( strcasecmp(a, "true") == 0 ) return 1;
+	if ( strcasecmp(a, "false") == 0 ) return 0;
+	return atoi(a);
+}
+
+
+static int assplode_algebraic(const char *a_orig, char ***pbits)
+{
+	int len, i;
+	int nexp;
+	char **bits;
+	char *a;
+	int idx, istr;
+
+	len = strlen(a_orig);
+
+	/* Add plus at start if no sign there already */
+	if ( (a_orig[0] != '+') && (a_orig[0] != '-') ) {
+		len += 1;
+		a = malloc(len+1);
+		snprintf(a, len+1, "+%s", a_orig);
+		a[len] = '\0';
+
+	} else {
+		a = strdup(a_orig);
+	}
+
+	/* Count the expressions */
+	nexp = 0;
+	for ( i=0; i<len; i++ ) {
+		if ( (a[i] == '+') || (a[i] == '-') ) nexp++;
+	}
+
+	bits = calloc(nexp, sizeof(char *));
+
+	/* Break the string up */
+	idx = -1;
+	istr = 0;
+	assert((a[0] == '+') || (a[0] == '-'));
+	for ( i=0; i<len; i++ ) {
+
+		char ch;
+
+		ch = a[i];
+
+		if ( (ch == '+') || (ch == '-') ) {
+			if ( idx >= 0 ) bits[idx][istr] = '\0';
+			idx++;
+			bits[idx] = malloc(len+1);
+			istr = 0;
+		}
+
+		if ( !isdigit(ch) && (ch != '.') && (ch != 'x') && (ch != 'y')
+		  && (ch != '+') && (ch != '-') )
+		{
+			ERROR("Invalid character '%C' found.\n", ch);
+			return 0;
+		}
+
+		assert(idx >= 0);
+		bits[idx][istr++] = ch;
+
+	}
+	if ( idx >= 0 ) bits[idx][istr] = '\0';
+
+	*pbits = bits;
+	free(a);
+
+	return nexp;
+}
+
+
+/* Parses the scan directions (accounting for possible rotation)
+ * Assumes all white spaces have been already removed */
+static int dir_conv(const char *a, double *sx, double *sy)
+{
+	int n;
+	char **bits;
+	int i;
+
+	*sx = 0.0;  *sy = 0.0;
+
+	n = assplode_algebraic(a, &bits);
+
+	if ( n == 0 ) {
+		ERROR("Invalid direction '%s'\n", a);
+		return 1;
+	}
+
+	for ( i=0; i<n; i++ ) {
+
+		int len;
+		double val;
+		char axis;
+		int j;
+
+		len = strlen(bits[i]);
+		assert(len != 0);
+		axis = bits[i][len-1];
+		if ( (axis != 'x') && (axis != 'y') ) {
+			ERROR("Invalid symbol '%C' - must be x or y.\n", axis);
+			return 1;
+		}
+
+		/* Chop off the symbol now it's dealt with */
+		bits[i][len-1] = '\0';
+
+		/* Check for anything that isn't part of a number */
+		for ( j=0; j<strlen(bits[i]); j++ ) {
+			if ( isdigit(bits[i][j]) ) continue;
+			if ( bits[i][j] == '+' ) continue;
+			if ( bits[i][j] == '-' ) continue;
+			if ( bits[i][j] == '.' ) continue;
+			ERROR("Invalid coefficient '%s'\n", bits[i]);
+		}
+
+		if ( strlen(bits[i]) == 0 ) {
+			val = 1.0;
+		} else {
+			val = atof(bits[i]);
+		}
+		if ( strlen(bits[i]) == 1 ) {
+			if ( bits[i][0] == '+' ) val = 1.0;
+			if ( bits[i][0] == '-' ) val = -1.0;
+		}
+		if ( axis == 'x' ) {
+			*sx += val;
+		} else if ( axis == 'y' ) {
+			*sy += val;
+		}
+
+		free(bits[i]);
+
+	}
+	free(bits);
+
+	//STATUS("'%s' -> %5.2fx + %5.2fy\n", a, *sx, *sy);
+
+	return 0;
+}
+
+
+struct rvec get_q_for_panel(struct panel *p, double fs, double ss,
+                            double *ttp, double k)
+{
+	struct rvec q;
+	double twotheta, r, az;
+	double rx, ry;
+	double xs, ys;
+
+	/* Convert xs and ys, which are in fast scan/slow scan coordinates,
+	 * to x and y */
+	xs = fs*p->fsx + ss*p->ssx;
+	ys = fs*p->fsy + ss*p->ssy;
+
+	rx = (xs + p->cnx) / p->res;
+	ry = (ys + p->cny) / p->res;
+
+	/* Calculate q-vector for this sub-pixel */
+	r = sqrt(pow(rx, 2.0) + pow(ry, 2.0));
+
+	twotheta = atan2(r, p->clen);
+	az = atan2(ry, rx);
+	if ( ttp != NULL ) *ttp = twotheta;
+
+	q.u = k * sin(twotheta)*cos(az);
+	q.v = k * sin(twotheta)*sin(az);
+	q.w = k * (cos(twotheta) - 1.0);
+
+	return q;
+}
+
+
+struct rvec get_q(struct image *image, double fs, double ss,
+                  double *ttp, double k)
+{
+	struct panel *p;
+	const unsigned int fsi = fs;
+	const unsigned int ssi = ss;  /* Explicit rounding */
+
+	/* Determine which panel to use */
+	p = find_panel(image->det, fsi, ssi);
+	assert(p != NULL);
+
+	return get_q_for_panel(p, fs-(double)p->min_fs, ss-(double)p->min_ss,
+	                       ttp, k);
+}
+
+
+int in_bad_region(struct detector *det, double fs, double ss)
+{
+	double rx, ry;
+	struct panel *p;
+	double xs, ys;
+	int i;
+
+	/* Determine which panel to use */
+	p = find_panel(det, fs, ss);
+
+	/* No panel found -> definitely bad! */
+	if ( p == NULL ) return 1;
+
+	/* Convert xs and ys, which are in fast scan/slow scan coordinates,
+	 * to x and y */
+	xs = (fs-(double)p->min_fs)*p->fsx + (ss-(double)p->min_ss)*p->ssx;
+	ys = (fs-(double)p->min_fs)*p->fsy + (ss-(double)p->min_ss)*p->ssy;
+
+	rx = xs + p->cnx;
+	ry = ys + p->cny;
+
+	for ( i=0; i<det->n_bad; i++ ) {
+		struct badregion *b = &det->bad[i];
+		if ( rx < b->min_x ) continue;
+		if ( rx > b->max_x ) continue;
+		if ( ry < b->min_y ) continue;
+		if ( ry > b->max_y ) continue;
+		return 1;
+	}
+
+	return 0;
+}
+
+
+double get_tt(struct image *image, double fs, double ss)
+{
+	double r, rx, ry;
+	struct panel *p;
+	double xs, ys;
+
+	p = find_panel(image->det, fs, ss);
+
+	/* Convert xs and ys, which are in fast scan/slow scan coordinates,
+	 * to x and y */
+	xs = (fs-p->min_fs)*p->fsx + (ss-p->min_ss)*p->ssx;
+	ys = (fs-p->min_fs)*p->fsy + (ss-p->min_ss)*p->ssy;
+
+	rx = (xs + p->cnx) / p->res;
+	ry = (ys + p->cny) / p->res;
+
+	r = sqrt(pow(rx, 2.0) + pow(ry, 2.0));
+
+	return atan2(r, p->clen);
+}
+
+
+void record_image(struct image *image, int do_poisson)
+{
+	int x, y;
+	double total_energy, energy_density;
+	double ph_per_e;
+	double area;
+	double max_tt = 0.0;
+	int n_inf1 = 0;
+	int n_neg1 = 0;
+	int n_nan1 = 0;
+	int n_inf2 = 0;
+	int n_neg2 = 0;
+	int n_nan2 = 0;
+
+	/* How many photons are scattered per electron? */
+	area = M_PI*pow(image->beam->beam_radius, 2.0);
+	total_energy = image->beam->fluence * ph_lambda_to_en(image->lambda);
+	energy_density = total_energy / area;
+	ph_per_e = (image->beam->fluence /area) * pow(THOMSON_LENGTH, 2.0);
+	STATUS("Fluence = %8.2e photons, "
+	       "Energy density = %5.3f kJ/cm^2, "
+	       "Total energy = %5.3f microJ\n",
+	       image->beam->fluence, energy_density/1e7, total_energy*1e6);
+
+	for ( x=0; x<image->width; x++ ) {
+	for ( y=0; y<image->height; y++ ) {
+
+		double counts;
+		double cf;
+		double intensity, sa;
+		double pix_area, Lsq;
+		double xs, ys, rx, ry;
+		double dsq, proj_area;
+		struct panel *p;
+
+		intensity = (double)image->data[x + image->width*y];
+		if ( isinf(intensity) ) n_inf1++;
+		if ( intensity < 0.0 ) n_neg1++;
+		if ( isnan(intensity) ) n_nan1++;
+
+		p = find_panel(image->det, x, y);
+
+		/* Area of one pixel */
+		pix_area = pow(1.0/p->res, 2.0);
+		Lsq = pow(p->clen, 2.0);
+
+		/* Area of pixel as seen from crystal (approximate) */
+		proj_area = pix_area * cos(image->twotheta[x + image->width*y]);
+
+		/* Calculate distance from crystal to pixel */
+		xs = (x-p->min_fs)*p->fsx + (y-p->min_ss)*p->ssx;
+		ys = (x-p->min_fs)*p->fsy + (y-p->min_ss)*p->ssy;
+		rx = (xs + p->cnx) / p->res;
+		ry = (ys + p->cny) / p->res;
+		dsq = sqrt(pow(rx, 2.0) + pow(ry, 2.0));
+
+		/* Projected area of pixel divided by distance squared */
+		sa = proj_area / (dsq + Lsq);
+
+		if ( do_poisson ) {
+			counts = poisson_noise(intensity * ph_per_e
+			                              * sa * image->beam->dqe );
+		} else {
+			cf = intensity * ph_per_e * sa * image->beam->dqe;
+			counts = cf;
+		}
+
+		image->data[x + image->width*y] = counts
+		                                  * image->beam->adu_per_photon;
+
+		/* Sanity checks */
+		if ( isinf(image->data[x+image->width*y]) ) n_inf2++;
+		if ( isnan(image->data[x+image->width*y]) ) n_nan2++;
+		if ( image->data[x+image->width*y] < 0.0 ) n_neg2++;
+
+		if ( image->twotheta[x + image->width*y] > max_tt ) {
+			max_tt = image->twotheta[x + image->width*y];
+		}
+
+	}
+	progress_bar(x, image->width-1, "Post-processing");
+	}
+
+	STATUS("Max 2theta = %.2f deg, min d = %.2f nm\n",
+	        rad2deg(max_tt), (image->lambda/(2.0*sin(max_tt/2.0)))/1e-9);
+
+	double tt_side = image->twotheta[(image->width/2)+image->width*0];
+	STATUS("At middle of bottom edge: %.2f deg, min d = %.2f nm\n",
+	        rad2deg(tt_side), (image->lambda/(2.0*sin(tt_side/2.0)))/1e-9);
+
+	tt_side = image->twotheta[0+image->width*(image->height/2)];
+	STATUS("At middle of left edge: %.2f deg, min d = %.2f nm\n",
+	        rad2deg(tt_side), (image->lambda/(2.0*sin(tt_side/2.0)))/1e-9);
+
+	STATUS("Halve the d values to get the voxel size for a synthesis.\n");
+
+	if ( n_neg1 + n_inf1 + n_nan1 + n_neg2 + n_inf2 + n_nan2 ) {
+		ERROR("WARNING: The raw calculation produced %i negative"
+		      " values, %i infinities and %i NaNs.\n",
+		      n_neg1, n_inf1, n_nan1);
+		ERROR("WARNING: After processing, there were %i negative"
+		      " values, %i infinities and %i NaNs.\n",
+		      n_neg2, n_inf2, n_nan2);
+	}
+}
+
+
+struct panel *find_panel(struct detector *det, double fs, double ss)
+{
+	int p;
+
+	for ( p=0; p<det->n_panels; p++ ) {
+		if ( (fs >= det->panels[p].min_fs)
+		  && (fs <= det->panels[p].max_fs)
+		  && (ss >= det->panels[p].min_ss)
+		  && (ss <= det->panels[p].max_ss) ) {
+			return &det->panels[p];
+		}
+	}
+
+	return NULL;
+}
+
+
+int find_panel_number(struct detector *det, int fs, int ss)
+{
+	int p;
+
+	for ( p=0; p<det->n_panels; p++ ) {
+		if ( (fs >= det->panels[p].min_fs)
+		  && (fs <= det->panels[p].max_fs)
+		  && (ss >= det->panels[p].min_ss)
+		  && (ss <= det->panels[p].max_ss) ) return p;
+	}
+
+	return -1;
+}
+
+
+void fill_in_values(struct detector *det, struct hdfile *f)
+{
+	int i;
+
+	for ( i=0; i<det->n_panels; i++ ) {
+
+		struct panel *p = &det->panels[i];
+
+		if ( p->clen_from != NULL ) {
+			p->clen = get_value(f, p->clen_from) * 1.0e-3;
+			free(p->clen_from);
+			p->clen_from = NULL;
+		}
+
+		p->clen += p->coffset;
+
+	}
+}
+
+
+static struct panel *new_panel(struct detector *det, const char *name)
+{
+	struct panel *new;
+
+	det->n_panels++;
+	det->panels = realloc(det->panels, det->n_panels*sizeof(struct panel));
+
+	new = &det->panels[det->n_panels-1];
+	memcpy(new, &det->defaults, sizeof(struct panel));
+
+	/* Create a new copy of the camera length location if needed */
+	if ( new->clen_from != NULL ) {
+		new->clen_from = strdup(new->clen_from);
+	}
+	strcpy(new->name, name);
+
+	return new;
+}
+
+
+static struct badregion *new_bad_region(struct detector *det, const char *name)
+{
+	struct badregion *new;
+
+	det->n_bad++;
+	det->bad = realloc(det->bad, det->n_bad*sizeof(struct badregion));
+
+	new = &det->bad[det->n_bad-1];
+	new->min_x = NAN;
+	new->max_x = NAN;
+	new->min_y = NAN;
+	new->max_y = NAN;
+	strcpy(new->name, name);
+
+	return new;
+}
+
+
+struct panel *find_panel_by_name(struct detector *det, const char *name)
+{
+	int i;
+
+	for ( i=0; i<det->n_panels; i++ ) {
+		if ( strcmp(det->panels[i].name, name) == 0 ) {
+			return &det->panels[i];
+		}
+	}
+
+	return NULL;
+}
+
+
+static struct badregion *find_bad_region_by_name(struct detector *det,
+                                                 const char *name)
+{
+	int i;
+
+	for ( i=0; i<det->n_bad; i++ ) {
+		if ( strcmp(det->bad[i].name, name) == 0 ) {
+			return &det->bad[i];
+		}
+	}
+
+	return NULL;
+}
+
+
+static char *find_or_add_rg(struct detector *det, const char *name)
+{
+	int i;
+	char **new;
+	char *tmp;
+
+	for ( i=0; i<det->num_rigid_groups; i++ ) {
+
+		if ( strcmp(det->rigid_groups[i], name) == 0 ) {
+			return det->rigid_groups[i];
+		}
+
+	}
+
+	new = realloc(det->rigid_groups,
+	              (1+det->num_rigid_groups)*sizeof(char *));
+	if ( new == NULL ) return NULL;
+
+	det->rigid_groups = new;
+
+	tmp = strdup(name);
+	det->rigid_groups[det->num_rigid_groups] = tmp;
+
+	det->num_rigid_groups++;
+
+	return tmp;
+}
+
+
+static int parse_field_for_panel(struct panel *panel, const char *key,
+                                 const char *val, struct detector *det)
+{
+	int reject = 0;
+
+	if ( strcmp(key, "min_fs") == 0 ) {
+		panel->min_fs = atof(val);
+	} else if ( strcmp(key, "max_fs") == 0 ) {
+		panel->max_fs = atof(val);
+	} else if ( strcmp(key, "min_ss") == 0 ) {
+		panel->min_ss = atof(val);
+	} else if ( strcmp(key, "max_ss") == 0 ) {
+		panel->max_ss = atof(val);
+	} else if ( strcmp(key, "corner_x") == 0 ) {
+		panel->cnx = atof(val);
+	} else if ( strcmp(key, "corner_y") == 0 ) {
+		panel->cny = atof(val);
+	} else if ( strcmp(key, "rigid_group") == 0 ) {
+		panel->rigid_group = find_or_add_rg(det, val);
+	} else if ( strcmp(key, "clen") == 0 ) {
+
+		char *end;
+		double v = strtod(val, &end);
+		if ( end == val ) {
+			/* This means "fill in later" */
+			panel->clen = -1.0;
+			panel->clen_from = strdup(val);
+		} else {
+			panel->clen = v;
+			panel->clen_from = NULL;
+		}
+
+	} else if ( strcmp(key, "coffset") == 0) {
+		panel->coffset = atof(val);
+	} else if ( strcmp(key, "res") == 0 ) {
+		panel->res = atof(val);
+	} else if ( strcmp(key, "peak_sep") == 0 ) {
+		panel->peak_sep = atof(val);
+	} else if ( strcmp(key, "integr_radius") == 0 ) {
+		panel->integr_radius = atof(val);
+	} else if ( strcmp(key, "badrow_direction") == 0 ) {
+		panel->badrow = val[0]; /* First character only */
+		if ( (panel->badrow != 'x') && (panel->badrow != 'y')
+		  && (panel->badrow != 'f') && (panel->badrow != 's')
+		  && (panel->badrow != '-') ) {
+			ERROR("badrow_direction must be x, y, f, s or '-'\n");
+			ERROR("Assuming '-'\n.");
+			panel->badrow = '-';
+		}
+		if ( panel->badrow == 'x' ) panel->badrow = 'f';
+		if ( panel->badrow == 'y' ) panel->badrow = 's';
+	} else if ( strcmp(key, "no_index") == 0 ) {
+		panel->no_index = atob(val);
+	} else if ( strcmp(key, "fs") == 0 ) {
+		if ( dir_conv(val, &panel->fsx, &panel->fsy) != 0 ) {
+			ERROR("Invalid fast scan direction '%s'\n", val);
+			reject = 1;
+		}
+	} else if ( strcmp(key, "ss") == 0 ) {
+		if ( dir_conv(val, &panel->ssx, &panel->ssy) != 0 ) {
+			ERROR("Invalid slow scan direction '%s'\n", val);
+			reject = 1;
+		}
+	} else {
+		ERROR("Unrecognised field '%s'\n", key);
+	}
+
+	return reject;
+}
+
+
+static int parse_field_bad(struct badregion *panel, const char *key,
+                                  const char *val)
+{
+	int reject = 0;
+
+	if ( strcmp(key, "min_x") == 0 ) {
+		panel->min_x = atof(val);
+	} else if ( strcmp(key, "max_x") == 0 ) {
+		panel->max_x = atof(val);
+	} else if ( strcmp(key, "min_y") == 0 ) {
+		panel->min_y = atof(val);
+	} else if ( strcmp(key, "max_y") == 0 ) {
+		panel->max_y = atof(val);
+	} else {
+		ERROR("Unrecognised field '%s'\n", key);
+	}
+
+	return reject;
+}
+
+
+static void parse_toplevel(struct detector *det, const char *key,
+                           const char *val)
+{
+	if ( strcmp(key, "mask") == 0 ) {
+
+		det->mask = strdup(val);
+
+	} else if ( strcmp(key, "mask_bad") == 0 ) {
+
+		char *end;
+		double v = strtod(val, &end);
+
+		if ( end != val ) {
+			det->mask_bad = v;
+		}
+
+	} else if ( strcmp(key, "mask_good") == 0 ) {
+
+		char *end;
+		double v = strtod(val, &end);
+
+		if ( end != val ) {
+			det->mask_good = v;
+		}
+
+	} else if ( strcmp(key, "peak_sep") == 0 ) {
+		det->defaults.peak_sep = atof(val);
+	} else if ( strcmp(key, "integr_radius") == 0 ) {
+		det->defaults.integr_radius = atof(val);
+	} else if ( parse_field_for_panel(&det->defaults, key, val, det) ) {
+		ERROR("Unrecognised top level field '%s'\n", key);
+	}
+}
+
+
+struct detector *get_detector_geometry(const char *filename)
+{
+	FILE *fh;
+	struct detector *det;
+	char *rval;
+	char **bits;
+	int i;
+	int reject = 0;
+	int x, y, max_fs, max_ss;
+
+	fh = fopen(filename, "r");
+	if ( fh == NULL ) return NULL;
+
+	det = calloc(1, sizeof(struct detector));
+	if ( det == NULL ) {
+		fclose(fh);
+		return NULL;
+	}
+
+	det->n_panels = 0;
+	det->panels = NULL;
+	det->n_bad = 0;
+	det->bad = NULL;
+	det->mask_good = 0;
+	det->mask_bad = 0;
+	det->mask = NULL;
+	det->num_rigid_groups = 0;
+	det->rigid_groups = NULL;
+
+	/* The default defaults... */
+	det->defaults.min_fs = -1;
+	det->defaults.min_ss = -1;
+	det->defaults.max_fs = -1;
+	det->defaults.max_ss = -1;
+	det->defaults.cnx = NAN;
+	det->defaults.cny = NAN;
+	det->defaults.clen = -1.0;
+	det->defaults.coffset = 0.0;
+	det->defaults.res = -1.0;
+	det->defaults.badrow = '-';
+	det->defaults.no_index = 0;
+	det->defaults.peak_sep = 50.0;
+	det->defaults.integr_radius = 3.0;
+	det->defaults.fsx = 1.0;
+	det->defaults.fsy = 0.0;
+	det->defaults.ssx = 0.0;
+	det->defaults.ssy = 1.0;
+	det->defaults.rigid_group = NULL;
+	strncpy(det->defaults.name, "", 1023);
+
+	do {
+
+		int n1, n2;
+		char **path;
+		char line[1024];
+		struct badregion *badregion = NULL;
+		struct panel *panel = NULL;
+		char *key;
+		char wholeval[1024];
+
+		rval = fgets(line, 1023, fh);
+		if ( rval == NULL ) break;
+		chomp(line);
+
+		if ( line[0] == ';' ) continue;
+
+		n1 = assplode(line, " \t", &bits, ASSPLODE_NONE);
+		if ( n1 < 3 ) {
+			for ( i=0; i<n1; i++ ) free(bits[i]);
+			free(bits);
+			continue;
+		}
+
+		/* Stitch the pieces of the "value" back together */
+		wholeval[0] = '\0';  /* Empty string */
+		for ( i=2; i<n1; i++ ) {
+			if ( bits[i][0] == ';' ) break;  /* Stop on comment */
+			strncat(wholeval, bits[i], 1023);
+		}
+
+		if ( bits[1][0] != '=' ) {
+			for ( i=0; i<n1; i++ ) free(bits[i]);
+			free(bits);
+			continue;
+		}
+
+		n2 = assplode(bits[0], "/\\.", &path, ASSPLODE_NONE);
+		if ( n2 < 2 ) {
+			/* This was a top-level option, not handled above. */
+			parse_toplevel(det, bits[0], bits[2]);
+			for ( i=0; i<n1; i++ ) free(bits[i]);
+			free(bits);
+			for ( i=0; i<n2; i++ ) free(path[i]);
+			free(path);
+			continue;
+		}
+
+		if ( strncmp(path[0], "bad", 3) == 0 ) {
+			badregion = find_bad_region_by_name(det, path[0]);
+			if ( badregion == NULL ) {
+				badregion = new_bad_region(det, path[0]);
+			}
+		} else {
+			panel = find_panel_by_name(det, path[0]);
+			if ( panel == NULL ) {
+				panel = new_panel(det, path[0]);
+			}
+		}
+
+		key = path[1];
+
+		if ( panel != NULL ) {
+			if ( parse_field_for_panel(panel, path[1],
+			                           wholeval, det) )
+			{
+				reject = 1;
+			}
+		} else {
+			if ( parse_field_bad(badregion, path[1], wholeval) ) {
+				reject = 1;
+			}
+		}
+
+		for ( i=0; i<n1; i++ ) free(bits[i]);
+		for ( i=0; i<n2; i++ ) free(path[i]);
+		free(bits);
+		free(path);
+
+	} while ( rval != NULL );
+
+	if ( det->n_panels == -1 ) {
+		ERROR("No panel descriptions in geometry file.\n");
+		fclose(fh);
+		free(det);
+		return NULL;
+	}
+
+	max_fs = 0;
+	max_ss = 0;
+	for ( i=0; i<det->n_panels; i++ ) {
+
+		if ( det->panels[i].min_fs < 0 ) {
+			ERROR("Please specify the minimum FS coordinate for"
+			      " panel %s\n", det->panels[i].name);
+			reject = 1;
+		}
+		if ( det->panels[i].max_fs < 0 ) {
+			ERROR("Please specify the maximum FS coordinate for"
+			      " panel %s\n", det->panels[i].name);
+			reject = 1;
+		}
+		if ( det->panels[i].min_ss < 0 ) {
+			ERROR("Please specify the minimum SS coordinate for"
+			      " panel %s\n", det->panels[i].name);
+			reject = 1;
+		}
+		if ( det->panels[i].max_ss < 0 ) {
+			ERROR("Please specify the maximum SS coordinate for"
+			      " panel %s\n", det->panels[i].name);
+			reject = 1;
+		}
+		if ( isnan(det->panels[i].cnx)  ) {
+			ERROR("Please specify the corner X coordinate for"
+			      " panel %s\n", det->panels[i].name);
+			reject = 1;
+		}
+		if ( isnan(det->panels[i].cny) ) {
+			ERROR("Please specify the corner Y coordinate for"
+			      " panel %s\n", det->panels[i].name);
+			reject = 1;
+		}
+		if ( (det->panels[i].clen < 0.0)
+		  && (det->panels[i].clen_from == NULL) ) {
+			ERROR("Please specify the camera length for"
+			      " panel %s\n", det->panels[i].name);
+			reject = 1;
+		}
+		if ( det->panels[i].res < 0 ) {
+			ERROR("Please specify the resolution for"
+			      " panel %s\n", det->panels[i].name);
+			reject = 1;
+		}
+		/* It's OK if the badrow direction is '0' */
+		/* It's not a problem if "no_index" is still zero */
+		/* The default peak_sep is OK (maybe) */
+		/* The default transformation matrix is at least valid */
+
+		if ( det->panels[i].max_fs > max_fs ) {
+			max_fs = det->panels[i].max_fs;
+		}
+		if ( det->panels[i].max_ss > max_ss ) {
+			max_ss = det->panels[i].max_ss;
+		}
+
+	}
+
+	for ( i=0; i<det->n_bad; i++ ) {
+
+		if ( isnan(det->bad[i].min_x) ) {
+			ERROR("Please specify the minimum x coordinate for"
+			      " bad region %s\n", det->bad[i].name);
+			reject = 1;
+		}
+		if ( isnan(det->bad[i].min_y) ) {
+			ERROR("Please specify the minimum y coordinate for"
+			      " bad region %s\n", det->bad[i].name);
+			reject = 1;
+		}
+		if ( isnan(det->bad[i].max_x) ) {
+			ERROR("Please specify the maximum x coordinate for"
+			      " bad region %s\n", det->bad[i].name);
+			reject = 1;
+		}
+		if ( isnan(det->bad[i].max_y) ) {
+			ERROR("Please specify the maximum y coordinate for"
+			      " bad region %s\n", det->bad[i].name);
+			reject = 1;
+		}
+	}
+
+	for ( x=0; x<=max_fs; x++ ) {
+	for ( y=0; y<=max_ss; y++ ) {
+		if ( find_panel(det, x, y) == NULL ) {
+			ERROR("Detector geometry invalid: contains gaps.\n");
+			reject = 1;
+			goto out;
+		}
+	}
+	}
+
+out:
+	det->max_fs = max_fs;
+	det->max_ss = max_ss;
+
+	/* Calculate matrix inverse */
+	for ( i=0; i<det->n_panels; i++ ) {
+
+		struct panel *p;
+		double d;
+
+		p = &det->panels[i];
+
+		if ( p->fsx*p->ssy == p->ssx*p->fsy ) {
+			ERROR("Panel %i transformation singular.\n", i);
+			reject = 1;
+		}
+
+		d = (double)p->fsx*p->ssy - p->ssx*p->fsy;
+		p->xfs = p->ssy / d;
+		p->yfs = -p->ssx / d;
+		p->xss = -p->fsy / d;
+		p->yss = p->fsx / d;
+
+	}
+
+	if ( reject ) return NULL;
+
+	fclose(fh);
+
+	return det;
+}
+
+
+void free_detector_geometry(struct detector *det)
+{
+	int i;
+
+	for ( i=0; i<det->num_rigid_groups; i++ ) {
+		free(det->rigid_groups[i]);
+	}
+	free(det->rigid_groups);
+
+	free(det->panels);
+	free(det->bad);
+	free(det->mask);
+	free(det);
+}
+
+
+struct detector *copy_geom(const struct detector *in)
+{
+	struct detector *out;
+	int i;
+
+	out = malloc(sizeof(struct detector));
+	memcpy(out, in, sizeof(struct detector));
+
+	if ( in->mask != NULL ) {
+		out->mask = strdup(in->mask);
+	} else {
+		out->mask = NULL;  /* = in->mask */
+	}
+
+	out->panels = malloc(out->n_panels * sizeof(struct panel));
+	memcpy(out->panels, in->panels, out->n_panels * sizeof(struct panel));
+
+	out->bad = malloc(out->n_bad * sizeof(struct badregion));
+	memcpy(out->bad, in->bad, out->n_bad * sizeof(struct badregion));
+
+	if ( in->rigid_groups != NULL ) {
+
+		out->rigid_groups = malloc(out->num_rigid_groups*sizeof(char *));
+		memcpy(out->rigid_groups, in->rigid_groups,
+		       out->num_rigid_groups*sizeof(char *));
+
+		for ( i=0; i<in->num_rigid_groups; i++ ) {
+			out->rigid_groups[i] = strdup(in->rigid_groups[i]);
+		}
+
+	}
+
+	for ( i=0; i<out->n_panels; i++ ) {
+
+		struct panel *p;
+
+		p = &out->panels[i];
+
+		if ( p->clen_from != NULL ) {
+			/* Make a copy of the clen_from fields unique to this
+			 * copy of the structure. */
+			p->clen_from = strdup(p->clen_from);
+		}
+
+	}
+
+	for ( i=0; i<in->num_rigid_groups; i++ ) {
+
+		int j;
+		char *rg = in->rigid_groups[i];
+		char *rgn = out->rigid_groups[i];
+
+		for ( j=0; j<in->n_panels; j++ ) {
+
+			if ( in->panels[j].rigid_group == rg ) {
+				out->panels[j].rigid_group = rgn;
+			}
+
+		}
+
+	}
+
+	return out;
+}
+
+
+struct detector *simple_geometry(const struct image *image)
+{
+	struct detector *geom;
+
+	geom = calloc(1, sizeof(struct detector));
+
+	geom->n_panels = 1;
+	geom->panels = calloc(1, sizeof(struct panel));
+
+	geom->panels[0].min_fs = 0;
+	geom->panels[0].max_fs = image->width-1;
+	geom->panels[0].min_ss = 0;
+	geom->panels[0].max_ss = image->height-1;
+	geom->panels[0].cnx = -image->width / 2.0;
+	geom->panels[0].cny = -image->height / 2.0;
+	geom->panels[0].rigid_group = NULL;
+
+	geom->panels[0].fsx = 1;
+	geom->panels[0].fsy = 0;
+	geom->panels[0].ssx = 0;
+	geom->panels[0].ssy = 1;
+
+	geom->panels[0].xfs = 1;
+	geom->panels[0].xss = 0;
+	geom->panels[0].yfs = 0;
+	geom->panels[0].yss = 1;
+
+	return geom;
+}
+
+
+int reverse_2d_mapping(double x, double y, double *pfs, double *pss,
+                       struct detector *det)
+{
+	int i;
+
+	for ( i=0; i<det->n_panels; i++ ) {
+
+		struct panel *p = &det->panels[i];
+		double cx, cy, fs, ss;
+
+		/* Get position relative to corner */
+		cx = x - p->cnx;
+		cy = y - p->cny;
+
+		/* Reverse the transformation matrix */
+		fs = cx*p->xfs + cy*p->yfs;
+		ss = cx*p->xss + cy*p->yss;
+
+		/* In range? */
+		if ( fs < 0 ) continue;
+		if ( ss < 0 ) continue;
+		if ( fs > (p->max_fs-p->min_fs+1) ) continue;
+		if ( ss > (p->max_ss-p->min_ss+1) ) continue;
+
+		*pfs = fs + p->min_fs;
+		*pss = ss + p->min_ss;
+		return 0;
+
+	}
+
+	return 1;
+}
+
+
+static void check_extents(struct panel p, double *min_x, double *min_y,
+                          double *max_x, double *max_y, double fs, double ss)
+{
+	double xs, ys, rx, ry;
+
+	xs = fs*p.fsx + ss*p.ssx;
+	ys = fs*p.fsy + ss*p.ssy;
+
+	rx = xs + p.cnx;
+	ry = ys + p.cny;
+
+	if ( rx > *max_x ) *max_x = rx;
+	if ( ry > *max_y ) *max_y = ry;
+	if ( rx < *min_x ) *min_x = rx;
+	if ( ry < *min_y ) *min_y = ry;
+}
+
+
+double largest_q(struct image *image)
+{
+	int fs, ss;
+	double ttm = 0.0;
+	double qmax = 0.0;
+
+	for ( fs=0; fs<image->width; fs++ ) {
+	for ( ss=0; ss<image->height; ss++ ) {
+
+		struct rvec q;
+		double tt;
+
+		q = get_q(image, fs, ss, &tt, 1.0/image->lambda);
+
+		if ( tt > ttm ) {
+			qmax = modulus(q.u, q.v, q.w);
+			ttm = tt;
+		}
+
+	}
+	}
+
+	return qmax;
+}
+
+
+double smallest_q(struct image *image)
+{
+	int fs, ss;
+	double ttm  = +INFINITY;
+	double qmin = +INFINITY;
+	for ( fs=0; fs<image->width; fs++ ) {
+	for ( ss=0; ss<image->height; ss++ ) {
+
+		struct rvec q;
+		double tt;
+
+		q = get_q(image, fs, ss, &tt, 1.0/image->lambda);
+
+		if ( tt < ttm ) {
+			qmin = modulus(q.u, q.v, q.w);
+			ttm = tt;
+		}
+
+	}
+	}
+
+	return qmin;
+}
+
+
+void get_pixel_extents(struct detector *det,
+                       double *min_x, double *min_y,
+                       double *max_x, double *max_y)
+{
+	int i;
+
+	*min_x = 0.0;
+	*max_x = 0.0;
+	*min_y = 0.0;
+	*max_y = 0.0;
+
+	/* To determine the maximum extents of the detector, put all four
+	 * corners of each panel through the transformations and watch for the
+	 * biggest */
+
+	for ( i=0; i<det->n_panels; i++ ) {
+
+		check_extents(det->panels[i], min_x, min_y, max_x, max_y,
+		              0.0,
+		              0.0);
+
+		check_extents(det->panels[i], min_x, min_y, max_x, max_y,
+		              0.0,
+		              det->panels[i].max_ss-det->panels[i].min_ss+1);
+
+		check_extents(det->panels[i], min_x, min_y, max_x, max_y,
+		              det->panels[i].max_fs-det->panels[i].min_fs+1,
+		              0.0);
+
+		check_extents(det->panels[i], min_x, min_y, max_x, max_y,
+		              det->panels[i].max_fs-det->panels[i].min_fs+1,
+		              det->panels[i].max_ss-det->panels[i].min_ss+1);
+
+
+	}
+}
+
+
+int write_detector_geometry(const char *filename, struct detector *det)
+{
+	struct panel *p;
+	int pi;
+	FILE *fh;
+
+	if ( filename == NULL ) return 2;
+	if ( det->n_panels < 1 ) return 3;
+
+	fh = fopen(filename, "w");
+	if ( fh == NULL ) return 1;
+
+	for ( pi=0; pi<det->n_panels; pi++) {
+
+		p = &(det->panels[pi]);
+
+		if ( p == NULL ) return 4;
+
+		if ( pi > 0 ) fprintf(fh, "\n");
+
+		fprintf(fh, "%s/min_fs = %d\n", p->name, p->min_fs);
+		fprintf(fh, "%s/min_ss = %d\n", p->name, p->min_ss);
+		fprintf(fh, "%s/max_fs = %d\n", p->name, p->max_fs);
+		fprintf(fh, "%s/max_ss = %d\n", p->name, p->max_ss);
+		fprintf(fh, "%s/badrow_direction = %C\n", p->name, p->badrow);
+		fprintf(fh, "%s/res = %g\n", p->name, p->res);
+		fprintf(fh, "%s/peak_sep = %g\n", p->name, p->peak_sep);
+		fprintf(fh, "%s/clen = %s\n", p->name, p->clen_from);
+		fprintf(fh, "%s/fs = %+fx %+fy\n", p->name, p->fsx, p->fsy);
+		fprintf(fh, "%s/ss = %+fx %+fy\n", p->name, p->ssx, p->ssy);
+		fprintf(fh, "%s/corner_x = %g\n", p->name, p->cnx);
+		fprintf(fh, "%s/corner_y = %g\n", p->name, p->cny);
+
+		if ( p->no_index ) {
+			fprintf(fh, "%s/no_index = 1\n", p->name);
+		} /* else don't clutter up the file */
+
+		if ( p->rigid_group != NULL ) {
+			fprintf(fh, "%s/rigid_group = %s\n",
+			        p->name, p->rigid_group);
+		}
+
+	}
+	fclose(fh);
+
+	return 0;
+}
diff --git a/libcrystfel/src/detector.h b/libcrystfel/src/detector.h
new file mode 100644
index 00000000..dd5dccec
--- /dev/null
+++ b/libcrystfel/src/detector.h
@@ -0,0 +1,131 @@
+/*
+ * detector.h
+ *
+ * Detector properties
+ *
+ * (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
+
+#ifndef DETECTOR_H
+#define DETECTOR_H
+
+struct image;
+struct hdfile;
+
+#include "hdf5-file.h"
+#include "image.h"
+
+
+struct panel
+{
+	char     name[1024];  /* Name for this panel */
+
+	int      min_fs;  /* Smallest FS value considered to be in the panel */
+	int      max_fs;  /* Largest FS value considered to be in this panel */
+	int      min_ss;  /* ... and so on */
+	int      max_ss;
+	double   cnx;       /* Location of corner (min_fs,min_ss) in pixels */
+	double   cny;
+	double   coffset;
+	double   clen;     /* Camera length in metres */
+	char    *clen_from;
+	double   res;      /* Resolution in pixels per metre */
+	char     badrow;   /* 'x' or 'y' */
+	int      no_index; /* Don't index peaks in this panel if non-zero */
+	double   peak_sep; /* Characteristic peak separation */
+	double   integr_radius;  /* Peak integration radius */
+	char    *rigid_group;  /* Rigid group, or -1 for none */
+
+	double fsx;
+	double fsy;
+	double ssx;
+	double ssy;
+
+	double xfs;
+	double yfs;
+	double xss;
+	double yss;
+};
+
+
+struct badregion
+{
+	char name[1024];
+	double min_x;
+	double max_x;
+	double min_y;
+	double max_y;
+};
+
+
+struct detector
+{
+	struct panel     *panels;
+	int               n_panels;
+
+	int               max_fs;
+	int               max_ss;  /* Size of overall array needed, minus 1 */
+
+	struct badregion *bad;
+	int               n_bad;
+
+	char              *mask;
+	unsigned int       mask_bad;
+	unsigned int       mask_good;
+
+	char             **rigid_groups;
+	int                num_rigid_groups;
+
+	struct panel       defaults;
+};
+
+
+extern struct rvec get_q(struct image *image, double fs, double ss,
+                         double *ttp, double k);
+
+extern struct rvec get_q_for_panel(struct panel *p, double fs, double ss,
+                                   double *ttp, double k);
+
+extern double get_tt(struct image *image, double xs, double ys);
+
+extern int in_bad_region(struct detector *det, double fs, double ss);
+
+extern void record_image(struct image *image, int do_poisson);
+
+extern struct panel *find_panel(struct detector *det, double fs, double ss);
+
+extern int find_panel_number(struct detector *det, int fs, int ss);
+
+extern struct detector *get_detector_geometry(const char *filename);
+
+extern void free_detector_geometry(struct detector *det);
+
+extern struct detector *simple_geometry(const struct image *image);
+
+extern void get_pixel_extents(struct detector *det,
+                              double *min_x, double *min_y,
+                              double *max_x, double *max_y);
+
+extern void fill_in_values(struct detector *det, struct hdfile *f);
+
+extern struct detector *copy_geom(const struct detector *in);
+
+extern int reverse_2d_mapping(double x, double y, double *pfs, double *pss,
+                              struct detector *det);
+
+extern double largest_q(struct image *image);
+
+extern double smallest_q(struct image *image);
+
+extern struct panel *find_panel_by_name(struct detector *det, const char *name);
+
+
+#endif	/* DETECTOR_H */
diff --git a/libcrystfel/src/hdf5-file.c b/libcrystfel/src/hdf5-file.c
new file mode 100644
index 00000000..355e97f1
--- /dev/null
+++ b/libcrystfel/src/hdf5-file.c
@@ -0,0 +1,817 @@
+/*
+ * hdf5-file.c
+ *
+ * Read/write HDF5 data files
+ *
+ * (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 <stdint.h>
+#include <hdf5.h>
+#include <assert.h>
+
+#include "image.h"
+#include "hdf5-file.h"
+#include "utils.h"
+
+
+struct hdfile {
+
+	const char      *path;  /* Current data path */
+
+	size_t          nx;  /* Image width */
+	size_t          ny;  /* Image height */
+
+	hid_t           fh;  /* HDF file handle */
+	hid_t           dh;  /* Dataset handle */
+
+	int             data_open;  /* True if dh is initialised */
+};
+
+
+struct hdfile *hdfile_open(const char *filename)
+{
+	struct hdfile *f;
+
+	f = malloc(sizeof(struct hdfile));
+	if ( f == NULL ) return NULL;
+
+	/* Please stop spamming my terminal */
+	H5Eset_auto2(H5E_DEFAULT, NULL, NULL);
+
+	f->fh = H5Fopen(filename, H5F_ACC_RDONLY, H5P_DEFAULT);
+	if ( f->fh < 0 ) {
+		ERROR("Couldn't open file: %s\n", filename);
+		free(f);
+		return NULL;
+	}
+
+	f->data_open = 0;
+
+	return f;
+}
+
+
+int hdfile_set_image(struct hdfile *f, const char *path)
+{
+	hsize_t size[2];
+	hsize_t max_size[2];
+	hid_t sh;
+
+	f->dh = H5Dopen2(f->fh, path, H5P_DEFAULT);
+	if ( f->dh < 0 ) {
+		ERROR("Couldn't open dataset\n");
+		return -1;
+	}
+	f->data_open = 1;
+
+	sh = H5Dget_space(f->dh);
+	if ( H5Sget_simple_extent_ndims(sh) != 2 ) {
+		ERROR("Dataset is not two-dimensional\n");
+		return -1;
+	}
+	H5Sget_simple_extent_dims(sh, size, max_size);
+	H5Sclose(sh);
+
+	f->nx = size[0];
+	f->ny = size[1];
+
+	return 0;
+}
+
+
+int get_peaks(struct image *image, struct hdfile *f, const char *p)
+{
+	hid_t dh, sh;
+	hsize_t size[2];
+	hsize_t max_size[2];
+	int i;
+	float *buf;
+	herr_t r;
+	int tw;
+
+	dh = H5Dopen2(f->fh, p, H5P_DEFAULT);
+	if ( dh < 0 ) {
+		ERROR("Peak list (%s) not found.\n", p);
+		return 1;
+	}
+
+	sh = H5Dget_space(dh);
+	if ( sh < 0 ) {
+		H5Dclose(dh);
+		ERROR("Couldn't get dataspace for peak list.\n");
+		return 1;
+	}
+
+	if ( H5Sget_simple_extent_ndims(sh) != 2 ) {
+		ERROR("Peak list has the wrong dimensionality (%i).\n",
+		      H5Sget_simple_extent_ndims(sh));
+		H5Sclose(sh);
+		H5Dclose(dh);
+		return 1;
+	}
+
+	H5Sget_simple_extent_dims(sh, size, max_size);
+
+	tw = size[1];
+	if ( (tw != 3) && (tw != 4) ) {
+		H5Sclose(sh);
+		H5Dclose(dh);
+		ERROR("Peak list has the wrong dimensions.\n");
+		return 1;
+	}
+
+	buf = malloc(sizeof(float)*size[0]*size[1]);
+	if ( buf == NULL ) {
+		H5Sclose(sh);
+		H5Dclose(dh);
+		ERROR("Couldn't reserve memory for the peak list.\n");
+		return 1;
+	}
+	r = H5Dread(dh, H5T_NATIVE_FLOAT, H5S_ALL, H5S_ALL, H5P_DEFAULT, buf);
+	if ( r < 0 ) {
+		ERROR("Couldn't read peak list.\n");
+		free(buf);
+		return 1;
+	}
+
+	if ( image->features != NULL ) {
+		image_feature_list_free(image->features);
+	}
+	image->features = image_feature_list_new();
+
+	for ( i=0; i<size[0]; i++ ) {
+
+		float fs, ss, val;
+		struct panel *p;
+
+		fs = buf[tw*i+0];
+		ss = buf[tw*i+1];
+		val = buf[tw*i+2];
+
+		p = find_panel(image->det, fs, ss);
+		if ( p == NULL ) continue;
+		if ( p->no_index ) continue;
+
+		image_add_feature(image->features, fs, ss, image, val, NULL);
+
+	}
+
+	free(buf);
+	H5Sclose(sh);
+	H5Dclose(dh);
+
+	return 0;
+}
+
+
+static void cleanup(hid_t fh)
+{
+	int n_ids, i;
+	hid_t ids[256];
+
+	n_ids = H5Fget_obj_ids(fh, H5F_OBJ_ALL, 256, ids);
+	for ( i=0; i<n_ids; i++ ) {
+
+		hid_t id;
+		H5I_type_t type;
+
+		id = ids[i];
+		type = H5Iget_type(id);
+
+		if ( type == H5I_GROUP ) H5Gclose(id);
+		if ( type == H5I_DATASET ) H5Dclose(id);
+		if ( type == H5I_DATATYPE ) H5Tclose(id);
+		if ( type == H5I_DATASPACE ) H5Sclose(id);
+		if ( type == H5I_ATTR ) H5Aclose(id);
+
+	}
+}
+
+
+void hdfile_close(struct hdfile *f)
+{
+	if ( f->data_open ) {
+		H5Dclose(f->dh);
+	}
+	cleanup(f->fh);
+
+	H5Fclose(f->fh);
+	free(f);
+}
+
+
+int hdf5_write(const char *filename, const void *data,
+               int width, int height, int type)
+{
+	hid_t fh, gh, sh, dh;	/* File, group, dataspace and data handles */
+	hid_t ph;  /* Property list */
+	herr_t r;
+	hsize_t size[2];
+	hsize_t max_size[2];
+
+	fh = H5Fcreate(filename, H5F_ACC_TRUNC, H5P_DEFAULT, H5P_DEFAULT);
+	if ( fh < 0 ) {
+		ERROR("Couldn't create file: %s\n", filename);
+		return 1;
+	}
+
+	gh = H5Gcreate2(fh, "data", H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT);
+	if ( gh < 0 ) {
+		ERROR("Couldn't create group\n");
+		H5Fclose(fh);
+		return 1;
+	}
+
+	/* Note the "swap" here, according to section 3.2.5,
+	 * "C versus Fortran Dataspaces", of the HDF5 user's guide. */
+	size[0] = height;
+	size[1] = width;
+	max_size[0] = height;
+	max_size[1] = width;
+	sh = H5Screate_simple(2, size, max_size);
+
+	/* Set compression */
+	ph = H5Pcreate(H5P_DATASET_CREATE);
+	H5Pset_chunk(ph, 2, size);
+	H5Pset_deflate(ph, 3);
+
+	dh = H5Dcreate2(gh, "data", type, sh,
+	                H5P_DEFAULT, ph, H5P_DEFAULT);
+	if ( dh < 0 ) {
+		ERROR("Couldn't create dataset\n");
+		H5Fclose(fh);
+		return 1;
+	}
+
+	/* Muppet check */
+	H5Sget_simple_extent_dims(sh, size, max_size);
+
+	r = H5Dwrite(dh, type, H5S_ALL,
+	             H5S_ALL, H5P_DEFAULT, data);
+	if ( r < 0 ) {
+		ERROR("Couldn't write data\n");
+		H5Dclose(dh);
+		H5Fclose(fh);
+		return 1;
+	}
+
+	H5Pclose(ph);
+	H5Gclose(gh);
+	H5Dclose(dh);
+	H5Fclose(fh);
+
+	return 0;
+}
+
+
+static double get_wavelength(struct hdfile *f)
+{
+	herr_t r;
+	hid_t dh;
+	double lambda;
+	int nm = 1;
+
+	dh = H5Dopen2(f->fh, "/LCLS/photon_wavelength_nm", H5P_DEFAULT);
+	if ( dh < 0 ) {
+		dh = H5Dopen2(f->fh, "/LCLS/photon_wavelength_A", H5P_DEFAULT);
+		if ( dh < 0 ) {
+			ERROR("Couldn't get photon wavelength from HDF5 file.\n");
+			return -1.0;
+		}
+		nm = 0;
+	}
+
+	r = H5Dread(dh, H5T_NATIVE_DOUBLE, H5S_ALL, H5S_ALL,
+	            H5P_DEFAULT, &lambda);
+	H5Dclose(dh);
+
+	if ( r < 0 ) return -1.0;
+	if ( isnan(lambda) ) return -1.0;
+
+	/* Convert nm -> m */
+	if ( nm ) return lambda / 1.0e9;
+	return lambda / 1.0e10;
+}
+
+
+static double get_i0(struct hdfile *f)
+{
+	herr_t r;
+	hid_t dh;
+	double i0;
+
+	dh = H5Dopen2(f->fh, "/LCLS/f_11_ENRC", H5P_DEFAULT);
+	if ( dh < 0 ) return -1.0;
+
+	r = H5Dread(dh, H5T_NATIVE_DOUBLE, H5S_ALL, H5S_ALL,
+		            H5P_DEFAULT, &i0);
+	H5Dclose(dh);
+	if ( r < 0 ) return -1.0;
+
+	return i0;
+}
+
+
+static void debodge_saturation(struct hdfile *f, struct image *image)
+{
+	hid_t dh, sh;
+	hsize_t size[2];
+	hsize_t max_size[2];
+	int i;
+	float *buf;
+	herr_t r;
+
+	dh = H5Dopen2(f->fh, "/processing/hitfinder/peakinfo_saturated",
+	              H5P_DEFAULT);
+
+	if ( dh < 0 ) {
+		/* This isn't an error */
+		return;
+	}
+
+	sh = H5Dget_space(dh);
+	if ( sh < 0 ) {
+		H5Dclose(dh);
+		ERROR("Couldn't get dataspace for saturation table.\n");
+		return;
+	}
+
+	if ( H5Sget_simple_extent_ndims(sh) != 2 ) {
+		H5Sclose(sh);
+		H5Dclose(dh);
+		return;
+	}
+
+	H5Sget_simple_extent_dims(sh, size, max_size);
+
+	if ( size[1] != 3 ) {
+		H5Sclose(sh);
+		H5Dclose(dh);
+		ERROR("Saturation table has the wrong dimensions.\n");
+		return;
+	}
+
+	buf = malloc(sizeof(float)*size[0]*size[1]);
+	if ( buf == NULL ) {
+		H5Sclose(sh);
+		H5Dclose(dh);
+		ERROR("Couldn't reserve memory for saturation table.\n");
+		return;
+	}
+	r = H5Dread(dh, H5T_NATIVE_FLOAT, H5S_ALL, H5S_ALL, H5P_DEFAULT, buf);
+	if ( r < 0 ) {
+		ERROR("Couldn't read saturation table.\n");
+		free(buf);
+		return;
+	}
+
+	for ( i=0; i<size[0]; i++ ) {
+
+		unsigned int x, y;
+		float val;
+
+		x = buf[3*i+0];
+		y = buf[3*i+1];
+		val = buf[3*i+2];
+
+		image->data[x+image->width*y] = val / 5.0;
+		image->data[x+1+image->width*y] = val / 5.0;
+		image->data[x-1+image->width*y] = val / 5.0;
+		image->data[x+image->width*(y+1)] = val / 5.0;
+		image->data[x+image->width*(y-1)] = val / 5.0;
+
+	}
+
+	free(buf);
+	H5Sclose(sh);
+	H5Dclose(dh);
+}
+
+
+int hdf5_read(struct hdfile *f, struct image *image, int satcorr)
+{
+	herr_t r;
+	float *buf;
+	uint16_t *flags;
+	hid_t mask_dh;
+
+	/* Note the "swap" here, according to section 3.2.5,
+	 * "C versus Fortran Dataspaces", of the HDF5 user's guide. */
+	image->width = f->ny;
+	image->height = f->nx;
+
+	buf = malloc(sizeof(float)*f->nx*f->ny);
+
+	r = H5Dread(f->dh, H5T_NATIVE_FLOAT, H5S_ALL, H5S_ALL,
+	            H5P_DEFAULT, buf);
+	if ( r < 0 ) {
+		ERROR("Couldn't read data\n");
+		free(buf);
+		return 1;
+	}
+	image->data = buf;
+
+	if ( (image->det != NULL) && (image->det->mask != NULL) ) {
+
+		mask_dh = H5Dopen2(f->fh, image->det->mask, H5P_DEFAULT);
+		if ( mask_dh <= 0 ) {
+			ERROR("Couldn't open flags\n");
+			image->flags = NULL;
+		} else {
+			flags = malloc(sizeof(uint16_t)*f->nx*f->ny);
+			r = H5Dread(mask_dh, H5T_NATIVE_UINT16, H5S_ALL, H5S_ALL,
+				    H5P_DEFAULT, flags);
+			if ( r < 0 ) {
+				ERROR("Couldn't read flags\n");
+				free(flags);
+				image->flags = NULL;
+			} else {
+				image->flags = flags;
+			}
+			H5Dclose(mask_dh);
+		}
+
+	}
+
+	/* Read wavelength from file */
+	image->lambda = get_wavelength(f);
+
+	image->i0 = get_i0(f);
+	if ( image->i0 < 0.0 ) {
+		ERROR("Couldn't read incident intensity - using 1.0.\n");
+		image->i0 = 1.0;
+		image->i0_available = 0;
+	} else {
+		image->i0_available = 1;
+	}
+
+	if ( satcorr ) debodge_saturation(f, image);
+
+	return 0;
+}
+
+
+static int looks_like_image(hid_t h)
+{
+	hid_t sh;
+	hsize_t size[2];
+	hsize_t max_size[2];
+
+	sh = H5Dget_space(h);
+	if ( sh < 0 ) return 0;
+
+	if ( H5Sget_simple_extent_ndims(sh) != 2 ) {
+		return 0;
+	}
+
+	H5Sget_simple_extent_dims(sh, size, max_size);
+
+	if ( ( size[0] > 64 ) && ( size[1] > 64 ) ) return 1;
+
+	return 0;
+}
+
+
+double get_value(struct hdfile *f, const char *name)
+{
+	hid_t dh;
+	hid_t sh;
+	hsize_t size;
+	hsize_t max_size;
+	hid_t type;
+	hid_t class;
+	herr_t r;
+	double buf;
+
+	dh = H5Dopen2(f->fh, name, H5P_DEFAULT);
+	if ( dh < 0 ) {
+		ERROR("Couldn't open data\n");
+		return 0.0;
+	}
+
+	type = H5Dget_type(dh);
+	class = H5Tget_class(type);
+
+	if ( class != H5T_FLOAT ) {
+		ERROR("Not a floating point value.\n");
+		H5Tclose(type);
+		H5Dclose(dh);
+		return 0.0;
+	}
+
+	sh = H5Dget_space(dh);
+	if ( H5Sget_simple_extent_ndims(sh) != 1 ) {
+		ERROR("Not a scalar value.\n");
+		H5Tclose(type);
+		H5Dclose(dh);
+		return 0.0;
+	}
+
+	H5Sget_simple_extent_dims(sh, &size, &max_size);
+	if ( size != 1 ) {
+		ERROR("Not a scalar value.\n");
+		H5Tclose(type);
+		H5Dclose(dh);
+		return 0.0;
+	}
+
+	r = H5Dread(dh, H5T_NATIVE_DOUBLE, H5S_ALL, H5S_ALL,
+	            H5P_DEFAULT, &buf);
+	if ( r < 0 )  {
+		ERROR("Couldn't read value.\n");
+		H5Tclose(type);
+		H5Dclose(dh);
+		return 0.0;
+	}
+
+	return buf;
+}
+
+
+struct copy_hdf5_field
+{
+	char **fields;
+	int n_fields;
+	int max_fields;
+};
+
+
+struct copy_hdf5_field *new_copy_hdf5_field_list()
+{
+	struct copy_hdf5_field *n;
+
+	n = calloc(1, sizeof(struct copy_hdf5_field));
+	if ( n == NULL ) return NULL;
+
+	n->max_fields = 32;
+	n->fields = malloc(n->max_fields*sizeof(char *));
+	if ( n->fields == NULL ) {
+		free(n);
+		return NULL;
+	}
+
+	return n;
+}
+
+
+void free_copy_hdf5_field_list(struct copy_hdf5_field *n)
+{
+	int i;
+	for ( i=0; i<n->n_fields; i++ ) {
+		free(n->fields[i]);
+	}
+	free(n->fields);
+	free(n);
+}
+
+
+void add_copy_hdf5_field(struct copy_hdf5_field *copyme,
+                         const char *name)
+{
+	assert(copyme->n_fields < copyme->max_fields);  /* FIXME */
+
+	copyme->fields[copyme->n_fields] = strdup(name);
+	if ( copyme->fields[copyme->n_fields] == NULL ) {
+		ERROR("Failed to add field for copying '%s'\n", name);
+		return;
+	}
+
+	copyme->n_fields++;
+}
+
+
+void copy_hdf5_fields(struct hdfile *f, const struct copy_hdf5_field *copyme,
+                      FILE *fh)
+{
+	int i;
+
+	if ( copyme == NULL ) return;
+
+	for ( i=0; i<copyme->n_fields; i++ ) {
+
+		char *val;
+		char *field;
+
+		field = copyme->fields[i];
+		val = hdfile_get_string_value(f, field);
+
+		if ( field[0] == '/' ) {
+			fprintf(fh, "hdf5%s = %s\n", field, val);
+		} else {
+			fprintf(fh, "hdf5/%s = %s\n", field, val);
+		}
+
+		free(val);
+
+	}
+}
+
+
+char *hdfile_get_string_value(struct hdfile *f, const char *name)
+{
+	hid_t dh;
+	hid_t sh;
+	hsize_t size;
+	hsize_t max_size;
+	hid_t type;
+	hid_t class;
+
+	dh = H5Dopen2(f->fh, name, H5P_DEFAULT);
+	if ( dh < 0 ) return NULL;
+
+	type = H5Dget_type(dh);
+	class = H5Tget_class(type);
+
+	if ( class == H5T_STRING ) {
+
+		herr_t r;
+		char *tmp;
+		hid_t sh;
+
+		size = H5Tget_size(type);
+		tmp = malloc(size+1);
+
+		sh = H5Screate(H5S_SCALAR);
+
+		r = H5Dread(dh, type, sh, sh, H5P_DEFAULT, tmp);
+		if ( r < 0 ) goto fail;
+
+		/* Two possibilities:
+		 *   String is already zero-terminated
+		 *   String is not terminated.
+		 * Make sure things are done properly... */
+		tmp[size] = '\0';
+		chomp(tmp);
+
+		return tmp;
+
+	}
+
+	sh = H5Dget_space(dh);
+	if ( H5Sget_simple_extent_ndims(sh) != 1 ) goto fail;
+
+	H5Sget_simple_extent_dims(sh, &size, &max_size);
+	if ( size != 1 ) {
+		H5Dclose(dh);
+		goto fail;
+	}
+
+	switch ( class ) {
+	case H5T_FLOAT : {
+
+		herr_t r;
+		double buf;
+		char *tmp;
+
+		r = H5Dread(dh, H5T_NATIVE_DOUBLE, H5S_ALL, H5S_ALL,
+		            H5P_DEFAULT, &buf);
+		if ( r < 0 ) goto fail;
+
+		tmp = malloc(256);
+		snprintf(tmp, 255, "%f", buf);
+
+		return tmp;
+
+	}
+	case H5T_INTEGER : {
+
+		herr_t r;
+		int buf;
+		char *tmp;
+
+		r = H5Dread(dh, H5T_NATIVE_INT, H5S_ALL, H5S_ALL,
+		            H5P_DEFAULT, &buf);
+		if ( r < 0 ) goto fail;
+
+		tmp = malloc(256);
+		snprintf(tmp, 255, "%d", buf);
+
+		return tmp;
+	}
+	default : {
+		goto fail;
+	}
+	}
+
+fail:
+	H5Tclose(type);
+	H5Dclose(dh);
+	return NULL;
+}
+
+
+char **hdfile_read_group(struct hdfile *f, int *n, const char *parent,
+                        int **p_is_group, int **p_is_image)
+{
+	hid_t gh;
+	hsize_t num;
+	char **res;
+	int i;
+	int *is_group;
+	int *is_image;
+	H5G_info_t ginfo;
+
+	gh = H5Gopen2(f->fh, parent, H5P_DEFAULT);
+	if ( gh < 0 ) {
+		*n = 0;
+		return NULL;
+	}
+
+	if ( H5Gget_info(gh, &ginfo) < 0 ) {
+		/* Whoopsie */
+		*n = 0;
+		return NULL;
+	}
+	num = ginfo.nlinks;
+	*n = num;
+	if ( num == 0 ) return NULL;
+
+	res = malloc(num*sizeof(char *));
+	is_image = malloc(num*sizeof(int));
+	is_group = malloc(num*sizeof(int));
+	*p_is_image = is_image;
+	*p_is_group = is_group;
+
+	for ( i=0; i<num; i++ ) {
+
+		char buf[256];
+		hid_t dh;
+		H5I_type_t type;
+
+		H5Lget_name_by_idx(gh, ".", H5_INDEX_NAME, H5_ITER_NATIVE,
+		                   i, buf, 255, H5P_DEFAULT);
+		res[i] = malloc(256);
+		if ( strlen(parent) > 1 ) {
+			snprintf(res[i], 255, "%s/%s", parent, buf);
+		} else {
+			snprintf(res[i], 255, "%s%s", parent, buf);
+		} /* ick */
+
+		is_image[i] = 0;
+		is_group[i] = 0;
+		dh = H5Oopen(gh, buf, H5P_DEFAULT);
+		if ( dh < 0 ) continue;
+		type = H5Iget_type(dh);
+
+		if ( type == H5I_GROUP ) {
+			is_group[i] = 1;
+		} else if ( type == H5I_DATASET ) {
+			is_image[i] = looks_like_image(dh);
+		}
+		H5Oclose(dh);
+
+	}
+
+	return res;
+}
+
+
+int hdfile_set_first_image(struct hdfile *f, const char *group)
+{
+	char **names;
+	int *is_group;
+	int *is_image;
+	int n, i, j;
+
+	names = hdfile_read_group(f, &n, group, &is_group, &is_image);
+	if ( n == 0 ) return 1;
+
+	for ( i=0; i<n; i++ ) {
+
+		if ( is_image[i] ) {
+			hdfile_set_image(f, names[i]);
+			for ( j=0; j<n; j++ ) free(names[j]);
+			free(is_image);
+			free(is_group);
+			free(names);
+			return 0;
+		} else if ( is_group[i] ) {
+			if ( !hdfile_set_first_image(f, names[i]) ) {
+				for ( j=0; j<n; j++ ) free(names[j]);
+				free(is_image);
+				free(is_group);
+				free(names);
+				return 0;
+			}
+		}
+
+	}
+
+	for ( j=0; j<n; j++ ) free(names[j]);
+	free(is_image);
+	free(is_group);
+	free(names);
+
+	return 1;
+}
diff --git a/libcrystfel/src/hdf5-file.h b/libcrystfel/src/hdf5-file.h
new file mode 100644
index 00000000..385e919c
--- /dev/null
+++ b/libcrystfel/src/hdf5-file.h
@@ -0,0 +1,55 @@
+/*
+ * hdf5-file.h
+ *
+ * Read/write HDF5 data files
+ *
+ * (c) 2006-2011 Thomas White <taw@physics.org>
+ *
+ * Part of CrystFEL - crystallography with a FEL
+ *
+ */
+
+#ifdef HAVE_CONFIG_H
+#include <config.h>
+#endif
+
+#ifndef HDF5_H
+#define HDF5_H
+
+#include <stdint.h>
+#include <hdf5.h>
+
+#include "image.h"
+
+struct hdfile;
+
+struct copy_hdf5_field;
+
+
+extern int hdf5_write(const char *filename, const void *data,
+                      int width, int height, int type);
+
+extern int hdf5_read(struct hdfile *f, struct image *image, int satcorr);
+
+extern struct hdfile *hdfile_open(const char *filename);
+extern int hdfile_set_image(struct hdfile *f, const char *path);
+extern int16_t *hdfile_get_image_binned(struct hdfile *hdfile,
+                                         int binning, int16_t *maxp);
+extern char **hdfile_read_group(struct hdfile *f, int *n, const char *parent,
+                                int **p_is_group, int **p_is_image);
+extern int hdfile_set_first_image(struct hdfile *f, const char *group);
+extern void hdfile_close(struct hdfile *f);
+
+extern char *hdfile_get_string_value(struct hdfile *f, const char *name);
+extern int get_peaks(struct image *image, struct hdfile *f, const char *p);
+extern double get_value(struct hdfile *f, const char *name);
+
+extern struct copy_hdf5_field *new_copy_hdf5_field_list(void);
+extern void free_copy_hdf5_field_list(struct copy_hdf5_field *f);
+extern void copy_hdf5_fields(struct hdfile *f,
+                             const struct copy_hdf5_field *copyme, FILE *fh);
+extern void add_copy_hdf5_field(struct copy_hdf5_field *copyme,
+                                const char *name);
+
+
+#endif	/* HDF5_H */
diff --git a/libcrystfel/src/image.c b/libcrystfel/src/image.c
new file mode 100644
index 00000000..75a8af0a
--- /dev/null
+++ b/libcrystfel/src/image.c
@@ -0,0 +1,144 @@
+/*
+ * image.c
+ *
+ * Handle images and image features
+ *
+ * (c) 2006-2010 Thomas White <taw@physics.org>
+ *
+ * Part of CrystFEL - crystallography with a FEL
+ *
+ */
+
+
+#include <stdlib.h>
+#include <assert.h>
+#include <math.h>
+#include <stdio.h>
+
+#include "image.h"
+#include "utils.h"
+
+/**
+ * SECTION:image
+ * @short_description: Data structure representing an image
+ * @title: Image
+ * @section_id:
+ * @see_also:
+ * @include: "image.h"
+ * @Image:
+ *
+ * The <structname>image</structname> structure represents an image, usually one
+ * frame from a large series of diffraction patterns, which might be from the
+ * same or different crystals.
+ */
+
+
+struct _imagefeaturelist
+{
+	struct imagefeature	*features;
+	int			n_features;
+};
+
+
+void image_add_feature(ImageFeatureList *flist, double fs, double ss,
+                       struct image *parent, double intensity, const char *name)
+{
+	if ( flist->features ) {
+		flist->features = realloc(flist->features,
+		                    (flist->n_features+1)
+		                    *sizeof(struct imagefeature));
+	} else {
+		assert(flist->n_features == 0);
+		flist->features = malloc(sizeof(struct imagefeature));
+	}
+
+	flist->features[flist->n_features].fs = fs;
+	flist->features[flist->n_features].ss = ss;
+	flist->features[flist->n_features].intensity = intensity;
+	flist->features[flist->n_features].parent = parent;
+	flist->features[flist->n_features].name = name;
+	flist->features[flist->n_features].valid = 1;
+
+	flist->n_features++;
+
+}
+
+
+ImageFeatureList *image_feature_list_new()
+{
+	ImageFeatureList *flist;
+
+	flist = malloc(sizeof(ImageFeatureList));
+
+	flist->n_features = 0;
+	flist->features = NULL;
+
+	return flist;
+}
+
+
+void image_feature_list_free(ImageFeatureList *flist)
+{
+	if ( !flist ) return;
+
+	if ( flist->features ) free(flist->features);
+	free(flist);
+}
+
+
+struct imagefeature *image_feature_closest(ImageFeatureList *flist,
+                                           double fs, double ss,
+                                           double *d, int *idx)
+{
+	int i;
+	double dmin = +HUGE_VAL;
+	int closest = 0;
+
+	for ( i=0; i<flist->n_features; i++ ) {
+
+		double ds;
+
+		ds = distance(flist->features[i].fs, flist->features[i].ss,
+		              fs, ss);
+
+		if ( ds < dmin ) {
+			dmin = ds;
+			closest = i;
+		}
+
+	}
+
+	if ( dmin < +HUGE_VAL ) {
+		*d = dmin;
+		*idx = closest;
+		return &flist->features[closest];
+	}
+
+	*d = +INFINITY;
+	return NULL;
+}
+
+
+int image_feature_count(ImageFeatureList *flist)
+{
+	if ( flist == NULL ) return 0;
+	return flist->n_features;
+}
+
+
+struct imagefeature *image_get_feature(ImageFeatureList *flist, int idx)
+{
+	/* Sanity check */
+	if ( flist == NULL ) return NULL;
+	if ( idx > flist->n_features ) return NULL;
+
+	if ( flist->features[idx].valid == 0 ) return NULL;
+
+	return &flist->features[idx];
+}
+
+
+void image_remove_feature(ImageFeatureList *flist, int idx)
+{
+	flist->features[idx].valid = 0;
+}
diff --git a/libcrystfel/src/image.h b/libcrystfel/src/image.h
new file mode 100644
index 00000000..4d4ac2d7
--- /dev/null
+++ b/libcrystfel/src/image.h
@@ -0,0 +1,184 @@
+/*
+ * image.h
+ *
+ * Handle images and image features
+ *
+ * (c) 2006-2011 Thomas White <taw@physics.org>
+ *
+ * Part of CrystFEL - crystallography with a FEL
+ *
+ */
+
+
+#ifdef HAVE_CONFIG_H
+#include <config.h>
+#endif
+
+#ifndef IMAGE_H
+#define IMAGE_H
+
+#include <stdint.h>
+#include <complex.h>
+#include <sys/types.h>
+
+#include "utils.h"
+#include "cell.h"
+#include "detector.h"
+#include "reflist.h"
+
+
+#define MAX_CELL_CANDIDATES (32)
+
+
+/* Structure describing a feature in an image */
+struct imagefeature {
+
+	struct image                    *parent;
+	double                          fs;
+	double                          ss;
+	double                          intensity;
+
+	/* Reciprocal space coordinates (m^-1 of course) of this feature */
+	double                          rx;
+	double                          ry;
+	double                          rz;
+
+	/* Internal use only */
+	int                             valid;
+
+	const char                      *name;
+};
+
+/* An opaque type representing a list of image features */
+typedef struct _imagefeaturelist ImageFeatureList;
+
+
+/**
+ * image:
+ *
+ * <programlisting>
+ * struct image
+ * {
+ *    float                   *data;
+ *    uint16_t                *flags;
+ *    double                  *twotheta;
+ *
+ *    UnitCell                *indexed_cell;
+ *    UnitCell                *candidate_cells[MAX_CELL_CANDIDATES];
+ *    int                     ncells;
+
+ *    struct detector         *det;
+ *    struct beam_params      *beam;
+ *    char                    *filename;
+ *    const struct copy_hdf5_field *copyme;
+ *
+ *    int                     id;
+ *
+ *    double                  m;
+ *
+ *    double                  lambda;
+ *    double                  div;
+ *    double                  bw;
+ *    double                  i0;
+ *    int                     i0_available;
+ *    double                  osf;
+ *    double                  profile_radius;
+ *    int                     pr_dud;
+ *
+ *    int                     width;
+ *    int                     height;
+ *
+ *    RefList                 *reflections;
+ *
+ *    ImageFeatureList        *features;
+ * };
+ * </programlisting>
+ *
+ * The field <structfield>data</structfield> contains the raw image data, if it
+ * is currently available.  The data might be available throughout the
+ * processing of an experimental pattern, but it might not be available when
+ * simulating, scaling or merging patterns.  Similarly,
+ * <structfield>flags</structfield> contains an array of the same dimensions
+ * as <structfield>data</structfield> to contain the bad pixel flags.
+ * <structfield>twotheta</structfield> likewise contains an array of 2*theta
+ * (scattering angle) values in radians, since these values are generated as a
+ * by-product of the scattering vector calculation and can be used later for
+ * calculating intensities from differential scattering cross sections.
+ *
+ * <structfield>candidate_cells</structfield> is an array of unit cells directly
+ * returned by the low-level indexing system. <structfield>ncells</structfield>
+ * is the number of candidate unit cells which were found.  The maximum number
+ * of cells which may be returned is <function>MAX_CELL_CANDIDATES</function>.
+ * <structfield>indexed_cell</structfield> contains the "correct" unit cell
+ * after cell reduction or matching has been performed.  The job of the cell
+ * reduction is to convert the list of candidate cells into a single indexed
+ * cell, or <function>NULL</function> on failure.
+ *
+ * <structfield>copyme</structfield> represents a list of HDF5 fields to copy
+ * to the output stream.
+ **/
+struct image;
+
+struct image {
+
+	float                   *data;
+	uint16_t                *flags;
+	double                  *twotheta;
+
+	UnitCell                *indexed_cell;
+	UnitCell                *candidate_cells[MAX_CELL_CANDIDATES];
+	int                     ncells;
+
+	struct detector         *det;
+	struct beam_params      *beam;  /* The nominal beam parameters */
+	char                    *filename;
+	const struct copy_hdf5_field *copyme;
+
+	int                     id;   /* ID number of the thread
+	                               * handling this image */
+
+	/* Information about the crystal */
+	double                  m;  /* Mosaicity in radians */
+
+	/* Per-shot radiation values */
+	double                  lambda;        /* Wavelength in m */
+	double                  div;           /* Divergence in radians */
+	double                  bw;            /* Bandwidth as a fraction */
+	double                  i0;            /* Incident intensity */
+	int                     i0_available;  /* 0 if f0 wasn't available
+	                                        * from the input. */
+	double                  osf;           /* Overall scaling factor */
+	double                  profile_radius; /* Radius of reflection */
+	int                     pr_dud;        /* Post refinement failed */
+
+	int                     width;
+	int                     height;
+
+	/* Integrated (or about-to-be-integrated) reflections */
+	RefList                 *reflections;
+
+	/* Detected peaks */
+	ImageFeatureList        *features;
+
+};
+
+
+/* Feature lists */
+extern ImageFeatureList *image_feature_list_new(void);
+
+extern void image_feature_list_free(ImageFeatureList *flist);
+
+extern void image_add_feature(ImageFeatureList *flist, double x, double y,
+                              struct image *parent, double intensity,
+                              const char *name);
+
+extern void image_remove_feature(ImageFeatureList *flist, int idx);
+
+extern struct imagefeature *image_feature_closest(ImageFeatureList *flist,
+                                                  double fs, double ss,
+                                                  double *d, int *idx);
+
+extern int image_feature_count(ImageFeatureList *flist);
+extern struct imagefeature *image_get_feature(ImageFeatureList *flist, int idx);
+
+#endif	/* IMAGE_H */
diff --git a/libcrystfel/src/list_tmp.h b/libcrystfel/src/list_tmp.h
new file mode 100644
index 00000000..a524b2f9
--- /dev/null
+++ b/libcrystfel/src/list_tmp.h
@@ -0,0 +1,106 @@
+/*
+ * Template for creating indexed 3D lists of a given type, usually indexed
+ * as signed h,k,l values where -INDMAX<={h,k,l}<=+INDMAX.
+ *
+ * These are used, for example, for:
+ *  - a list of 'double complex' values for storing structure factors,
+ *  - a list of 'double' values for storing reflection intensities,
+ *  - a list of 'unsigned int' values for counts of some sort.
+ *
+ * When LABEL and TYPE are #defined appropriately, including this header
+ * defines functions such as:
+ *  - new_list_<LABEL>(), for creating a new list of the given type,
+ *  - set_<LABEL>(), for setting a value in a list,
+ *  - lookup_<LABEL>(), for retrieving values from a list,
+ * ..and so on.
+ *
+ * See src/utils.h for more information.
+ *
+ */
+
+#include <stdio.h>
+
+#define ERROR_T(...) fprintf(stderr, __VA_ARGS__)
+
+static inline void LABEL(integrate)(TYPE *ref, signed int h,
+                                    signed int k, signed int l,
+                                    TYPE i)
+{
+	int idx;
+
+	if ( (abs(h) > INDMAX) || (abs(k) > INDMAX) || (abs(l) > INDMAX) ) {
+		ERROR_T("\nReflection %i %i %i is out of range!\n", h, k, l);
+		ERROR_T("You need to re-configure INDMAX and re-run.\n");
+		exit(1);
+	}
+
+	if ( h < 0 ) h += IDIM;
+	if ( k < 0 ) k += IDIM;
+	if ( l < 0 ) l += IDIM;
+
+	idx = h + (IDIM*k) + (IDIM*IDIM*l);
+	ref[idx] += i;
+}
+
+
+static inline void LABEL(set)(TYPE *ref, signed int h,
+                              signed int k, signed int l,
+                              TYPE i)
+{
+	int idx;
+
+	if ( (abs(h) > INDMAX) || (abs(k) > INDMAX) || (abs(l) > INDMAX) ) {
+		ERROR_T("\nReflection %i %i %i is out of range!\n", h, k, l);
+		ERROR_T("You need to re-configure INDMAX and re-run.\n");
+		exit(1);
+	}
+
+	if ( h < 0 ) h += IDIM;
+	if ( k < 0 ) k += IDIM;
+	if ( l < 0 ) l += IDIM;
+
+	idx = h + (IDIM*k) + (IDIM*IDIM*l);
+	ref[idx] = i;
+}
+
+
+static inline TYPE LABEL(lookup)(const TYPE *ref, signed int h,
+                                 signed int k, signed int l)
+{
+	int idx;
+
+	if ( (abs(h) > INDMAX) || (abs(k) > INDMAX) || (abs(l) > INDMAX) ) {
+		ERROR_T("\nReflection %i %i %i is out of range!\n", h, k, l);
+		ERROR_T("You need to re-configure INDMAX and re-run.\n");
+		exit(1);
+	}
+
+	if ( h < 0 ) h += IDIM;
+	if ( k < 0 ) k += IDIM;
+	if ( l < 0 ) l += IDIM;
+
+	idx = h + (IDIM*k) + (IDIM*IDIM*l);
+	return ref[idx];
+}
+
+
+static inline TYPE *LABEL(new_list)(void)
+{
+	TYPE *r;
+	size_t r_size;
+	r_size = IDIM*IDIM*IDIM*sizeof(TYPE);
+	r = malloc(r_size);
+	memset(r, 0, r_size);
+	return r;
+}
+
+
+static inline void LABEL(zero_list)(TYPE *ref)
+{
+	memset(ref, 0, IDIM*IDIM*IDIM*sizeof(TYPE));
+}
+
+
+#undef LABEL
+#undef TYPE
+#undef ERROR_T
diff --git a/libcrystfel/src/reflist.c b/libcrystfel/src/reflist.c
new file mode 100644
index 00000000..18856c67
--- /dev/null
+++ b/libcrystfel/src/reflist.c
@@ -0,0 +1,1048 @@
+/*
+ * reflist.c
+ *
+ * Fast reflection/peak list
+ *
+ * (c) 2011 Thomas White <taw@physics.org>
+ *
+ * Part of CrystFEL - crystallography with a FEL
+ *
+ */
+
+
+#include <stdlib.h>
+#include <assert.h>
+#include <stdio.h>
+#include <pthread.h>
+
+#include "reflist.h"
+#include "utils.h"
+
+/**
+ * SECTION:reflist
+ * @short_description: The fast reflection list
+ * @title: RefList
+ * @section_id:
+ * @see_also:
+ * @include: "reflist.h"
+ * @Image:
+ *
+ * The fast reflection list stores reflections in an RB-tree indexed
+ * by the Miller indices h, k and l.  Any reflection can be found in a
+ * length of time which scales logarithmically with the number of reflections in
+ * the list.
+ *
+ * A RefList can contain any number of reflections, and can store more than
+ * one reflection with a given set of indices, for example when two distinct
+ * reflections are to be stored according to their asymmetric indices.
+ *
+ * There are getters and setters which can be used to get and set values for an
+ * individual reflection.  The reflection list does not calculate any values,
+ * only stores what it was given earlier.  As such, you will need to carefully
+ * examine which fields your prior processing steps have filled in.
+ */
+
+
+struct _refldata {
+
+	/* Symmetric indices (i.e. the "real" indices) */
+	signed int hs;
+	signed int ks;
+	signed int ls;
+
+	/* Partiality and related geometrical stuff */
+	double r1;  /* First excitation error */
+	double r2;  /* Second excitation error */
+	double p;   /* Partiality */
+	int clamp1; /* Clamp status for r1 */
+	int clamp2; /* Clamp status for r2 */
+
+	/* Location in image */
+	double fs;
+	double ss;
+
+	/* The distance from the exact Bragg position to the coordinates
+	 * given above. */
+	double excitation_error;
+
+	/* Non-zero if this reflection can be used for scaling */
+	int scalable;
+
+	/* Non-zero if this reflection should be used as a "guide star" for
+	 * post refinement */
+	int refinable;
+
+	/* Intensity */
+	double intensity;
+	double esd_i;
+
+	/* Phase */
+	double phase;
+	int have_phase;
+
+	/* Redundancy */
+	int redundancy;
+
+	/* User-specified temporary values */
+	double temp1;
+	double temp2;
+};
+
+
+enum _nodecol {
+	RED,
+	BLACK
+};
+
+
+struct _reflection {
+
+	/* Listy stuff */
+	unsigned int serial;          /* Unique serial number, key */
+	struct _reflection *child[2]; /* Child nodes */
+	struct _reflection *next;     /* Next and previous in doubly linked */
+	struct _reflection *prev;     /*  list of duplicate reflections */
+	enum _nodecol col;            /* Colour (red or black) */
+
+	/* Payload */
+	pthread_mutex_t lock;         /* Protects the contents of "data" */
+	struct _refldata data;
+};
+
+
+struct _reflist {
+
+	struct _reflection *head;
+	struct _reflection *tail;
+
+};
+
+
+#define SERIAL(h, k, l) ((((h)+256)<<18) + (((k)+256)<<9) + ((l)+256))
+#define GET_H(serial) ((((serial) & 0xfffc0000)>>18)-256)
+#define GET_K(serial) ((((serial) & 0x0003fe00)>>9)-256)
+#define GET_L(serial) (((serial) & 0x000001ff)-256)
+
+/**************************** Creation / deletion *****************************/
+
+static Reflection *new_node(unsigned int serial)
+{
+	Reflection *new;
+
+	new = calloc(1, sizeof(struct _reflection));
+	new->serial = serial;
+	new->next = NULL;
+	new->prev = NULL;
+	new->child[0] = NULL;
+	new->child[1] = NULL;
+	new->col = RED;
+	pthread_mutex_init(&new->lock, NULL);
+
+	return new;
+}
+
+
+/**
+ * reflist_new:
+ *
+ * Creates a new reflection list.
+ *
+ * Returns: the new reflection list, or NULL on error.
+ */
+RefList *reflist_new()
+{
+	RefList *new;
+
+	new = malloc(sizeof(struct _reflist));
+	if ( new == NULL ) return NULL;
+
+	new->head = NULL;
+
+	return new;
+}
+
+
+/**
+ * reflection_new:
+ * @h: The h index of the new reflection
+ * @k: The k index of the new reflection
+ * @l: The l index of the new reflection
+ *
+ * Creates a new individual reflection.  You'll probably want to use
+ * add_refl_to_list() at some later point.
+ */
+Reflection *reflection_new(signed int h, signed int k, signed int l)
+{
+	return new_node(SERIAL(h, k, l));
+}
+
+
+/**
+ * reflection_free:
+ * @refl: The reflection to free.
+ *
+ * Destroys an individual reflection.
+ */
+void reflection_free(Reflection *refl)
+{
+	pthread_mutex_destroy(&refl->lock);
+	free(refl);
+}
+
+
+static void recursive_free(Reflection *refl)
+{
+	if ( refl->child[0] != NULL ) recursive_free(refl->child[0]);
+	if ( refl->child[1] != NULL ) recursive_free(refl->child[1]);
+
+	while ( refl != NULL ) {
+		Reflection *next = refl->next;
+		reflection_free(refl);
+		refl = next;
+	}
+}
+
+
+/**
+ * reflist_free:
+ * @list: The reflection list to free.
+ *
+ * Destroys a reflection list.
+ */
+void reflist_free(RefList *list)
+{
+	if ( list == NULL ) return;
+	if ( list->head != NULL ) {
+		recursive_free(list->head);
+	} /* else empty list */
+	free(list);
+}
+
+
+/********************************** Search ************************************/
+
+/**
+ * find_refl:
+ * @list: The reflection list to search in
+ * @h: The 'h' index to search for
+ * @k: The 'k' index to search for
+ * @l: The 'l' index to search for
+ *
+ * This function finds the first reflection in 'list' with the given indices.
+ *
+ * Since a %RefList can contain multiple reflections with the same indices, you
+ * may need to use next_found_refl() to get the other reflections.
+ *
+ * Returns: The found reflection, or NULL if no reflection with the given
+ * indices could be found.
+ **/
+Reflection *find_refl(const RefList *list,
+                      signed int h, signed int k, signed int l)
+{
+	unsigned int search = SERIAL(h, k, l);
+	Reflection *refl;
+
+	if ( list->head == NULL ) return NULL;
+
+	/* Indices greater than or equal to 256 are filtered out when
+	 * reflections are added, so don't even bother looking.
+	 * (also, looking for such reflections causes trouble because the search
+	 * serial number would be invalid) */
+	if ( abs(h) >= 256 ) return NULL;
+	if ( abs(k) >= 256 ) return NULL;
+	if ( abs(l) >= 256 ) return NULL;
+
+	refl = list->head;
+
+	while ( refl != NULL ) {
+
+		if ( refl->serial == search ) {
+
+			assert(search == refl->serial);
+			assert(h == GET_H(refl->serial));
+			assert(k == GET_K(refl->serial));
+			assert(l == GET_L(refl->serial));
+			return refl;
+
+		} else {
+
+			int dir = search > refl->serial;
+			if ( refl->child[dir] != NULL ) {
+				refl = refl->child[dir];
+			} else {
+				/* Hit the bottom of the tree */
+				return NULL;
+			}
+
+		}
+
+	}
+
+	return NULL;
+}
+
+
+/**
+ * next_found_refl:
+ * @refl: A reflection returned by find_refl() or next_found_refl()
+ *
+ * This function returns the next reflection in @refl's list with the same
+ * indices.
+ *
+ * Returns: The found reflection, or NULL if there are no more reflections with
+ * the same indices.
+ **/
+Reflection *next_found_refl(Reflection *refl)
+{
+	if ( refl->next != NULL ) assert(refl->serial == refl->next->serial);
+
+	return refl->next;  /* Well, that was easy... */
+}
+
+
+/********************************** Getters ***********************************/
+
+/**
+ * get_excitation_error:
+ * @refl: A %Reflection
+ *
+ * Returns: The excitation error for the reflection.
+ **/
+double get_excitation_error(const Reflection *refl)
+{
+	return refl->data.excitation_error;
+}
+
+
+/**
+ * get_detector_pos:
+ * @refl: A %Reflection
+ * @fs: Location at which to store the fast scan offset of the reflection
+ * @ss: Location at which to store the slow scan offset of the reflection
+ *
+ **/
+void get_detector_pos(const Reflection *refl, double *fs, double *ss)
+{
+	*fs = refl->data.fs;
+	*ss = refl->data.ss;
+}
+
+
+/**
+ * get_indices:
+ * @refl: A %Reflection
+ * @h: Location at which to store the 'h' index of the reflection
+ * @k: Location at which to store the 'k' index of the reflection
+ * @l: Location at which to store the 'l' index of the reflection
+ *
+ **/
+void get_indices(const Reflection *refl,
+                 signed int *h, signed int *k, signed int *l)
+{
+	*h = GET_H(refl->serial);
+	*k = GET_K(refl->serial);
+	*l = GET_L(refl->serial);
+}
+
+
+/**
+ * get_symmetric_indices:
+ * @refl: A %Reflection
+ * @hs: Location at which to store the 'h' index of the reflection
+ * @ks: Location at which to store the 'k' index of the reflection
+ * @ls: Location at which to store the 'l' index of the reflection
+ *
+ * This function gives the symmetric indices, that is, the "real" indices before
+ * squashing down to the asymmetric reciprocal unit.  This may be useful if the
+ * list is indexed according to the asymmetric indices, but you still need
+ * access to the symmetric version.  This happens during post-refinement.
+ *
+ **/
+void get_symmetric_indices(const Reflection *refl,
+                                  signed int *hs, signed int *ks,
+                                  signed int *ls)
+{
+	*hs = refl->data.hs;
+	*ks = refl->data.ks;
+	*ls = refl->data.ls;
+}
+
+
+/**
+ * get_partiality:
+ * @refl: A %Reflection
+ *
+ * Returns: The partiality of the reflection.
+ **/
+double get_partiality(const Reflection *refl)
+{
+	return refl->data.p;
+}
+
+
+/**
+ * get_intensity:
+ * @refl: A %Reflection
+ *
+ * Returns: The intensity of the reflection.
+ **/
+double get_intensity(const Reflection *refl)
+{
+	return refl->data.intensity;
+}
+
+
+/**
+ * get_partial:
+ * @refl: A %Reflection
+ * @r1: Location at which to store the first excitation error
+ * @r2: Location at which to store the second excitation error
+ * @p: Location at which to store the partiality
+ * @clamp_low: Location at which to store the first clamp status
+ * @clamp_high: Location at which to store the second clamp status
+ *
+ * This function is used during post refinement (in conjunction with
+ * set_partial()) to get access to the details of the partiality calculation.
+ *
+ **/
+void get_partial(const Reflection *refl, double *r1, double *r2, double *p,
+                 int *clamp_low, int *clamp_high)
+{
+	*r1 = refl->data.r1;
+	*r2 = refl->data.r2;
+	*p = get_partiality(refl);
+	*clamp_low = refl->data.clamp1;
+	*clamp_high = refl->data.clamp2;
+}
+
+
+/**
+ * get_scalable:
+ * @refl: A %Reflection
+ *
+ * Returns: non-zero if this reflection can be scaled.
+ *
+ **/
+int get_scalable(const Reflection *refl)
+{
+	return refl->data.scalable;
+}
+
+
+/**
+ * get_refinable:
+ * @refl: A %Reflection
+ *
+ * Returns: non-zero if this reflection can be used for post refinement.
+ *
+ **/
+int get_refinable(const Reflection *refl)
+{
+	return refl->data.refinable;
+}
+
+
+/**
+ * get_redundancy:
+ * @refl: A %Reflection
+ *
+ * The redundancy of the reflection is the number of measurements that have been
+ * made of it.  Note that a redundancy of zero may have a special meaning, such
+ * as that the reflection was impossible to integrate.  Note further that each
+ * reflection in the list has its own redundancy, even if there are multiple
+ * copies of the reflection in the list.  The total number of reflection
+ * measurements should always be the sum of the redundancies in the entire list.
+ *
+ * Returns: the number of measurements of this reflection.
+ *
+ **/
+int get_redundancy(const Reflection *refl)
+{
+	return refl->data.redundancy;
+}
+
+
+/**
+ * get_esd_intensity:
+ * @refl: A %Reflection
+ *
+ * Returns: the standard error in the intensity measurement (as returned by
+ * get_intensity()) for this reflection.
+ *
+ **/
+double get_esd_intensity(const Reflection *refl)
+{
+	return refl->data.esd_i;
+}
+
+
+/**
+ * get_phase:
+ * @refl: A %Reflection
+ * @have_phase: Place to store a non-zero value if the phase is set, or NULL.
+ *
+ * Returns: the phase for this reflection.
+ *
+ **/
+double get_phase(const Reflection *refl, int *have_phase)
+{
+	if ( have_phase != NULL ) *have_phase = refl->data.have_phase;
+	return refl->data.phase;
+}
+
+
+/**
+ * get_temp1:
+ * @refl: A %Reflection
+ *
+ * The temporary values can be used according to the needs of the calling
+ * program.
+ *
+ * Returns: the first temporary value for this reflection.
+ *
+ **/
+double get_temp1(const Reflection *refl)
+{
+	return refl->data.temp1;
+}
+
+
+/**
+ * get_temp2:
+ * @refl: A %Reflection
+ *
+ * The temporary values can be used according to the needs of the calling
+ * program.
+ *
+ * Returns: the second temporary value for this reflection.
+ *
+ **/
+double get_temp2(const Reflection *refl)
+{
+	return refl->data.temp2;
+}
+
+
+/********************************** Setters ***********************************/
+
+/**
+ * copy_data:
+ * @to: %Reflection to copy data into
+ * @from: %Reflection to copy data from
+ *
+ * This function is used to copy the data (which is everything listed above in
+ * the list of getters and setters, apart from the indices themselves) from one
+ * reflection to another.  This might be used when creating a new list from an
+ * old one, perhaps using the asymmetric indices instead of the raw indices for
+ * the new list.
+ *
+ **/
+void copy_data(Reflection *to, const Reflection *from)
+{
+	memcpy(&to->data, &from->data, sizeof(struct _refldata));
+}
+
+
+/**
+ * set_detector_pos:
+ * @refl: A %Reflection
+ * @exerr: The excitation error for this reflection
+ * @fs: The fast scan offset of the reflection
+ * @ss: The slow scan offset of the reflection
+ *
+ **/
+void set_detector_pos(Reflection *refl, double exerr, double fs, double ss)
+{
+	refl->data.excitation_error = exerr;
+	refl->data.fs = fs;
+	refl->data.ss = ss;
+}
+
+
+/**
+ * set_partial:
+ * @refl: A %Reflection
+ * @r1: The first excitation error
+ * @r2: The second excitation error
+ * @p: The partiality
+ * @clamp_low: The first clamp status
+ * @clamp_high: The second clamp status
+ *
+ * This function is used during post refinement (in conjunction with
+ * get_partial()) to get access to the details of the partiality calculation.
+ *
+ **/
+void set_partial(Reflection *refl, double r1, double r2, double p,
+                 double clamp_low, double clamp_high)
+{
+	refl->data.r1 = r1;
+	refl->data.r2 = r2;
+	refl->data.p = p;
+	refl->data.clamp1 = clamp_low;
+	refl->data.clamp2 = clamp_high;
+}
+
+
+/**
+ * set_int:
+ * @refl: A %Reflection
+ * @intensity: The intensity for the reflection.
+ *
+ * Set the intensity for the reflection.  Note that retrieval is done with
+ * get_intensity().
+ **/
+void set_int(Reflection *refl, double intensity)
+{
+	refl->data.intensity = intensity;
+}
+
+
+/**
+ * set_scalable:
+ * @refl: A %Reflection
+ * @scalable: Non-zero if this reflection should be scaled.
+ *
+ **/
+void set_scalable(Reflection *refl, int scalable)
+{
+	refl->data.scalable = scalable;
+}
+
+
+/**
+ * set_refinable:
+ * @refl: A %Reflection
+ * @refinable: Non-zero if this reflection can be used for post refinement.
+ *
+ **/
+void set_refinable(Reflection *refl, int refinable)
+{
+	refl->data.refinable = refinable;
+}
+
+
+/**
+ * set_redundancy:
+ * @refl: A %Reflection
+ * @red: New redundancy for the reflection
+ *
+ * The redundancy of the reflection is the number of measurements that have been
+ * made of it.  Note that a redundancy of zero may have a special meaning, such
+ * as that the reflection was impossible to integrate.  Note further that each
+ * reflection in the list has its own redundancy, even if there are multiple
+ * copies of the reflection in the list.  The total number of reflection
+ * measurements should always be the sum of the redundancies in the entire list.
+ *
+ **/
+void set_redundancy(Reflection *refl, int red)
+{
+	refl->data.redundancy = red;
+}
+
+
+/**
+ * set_esd_intensity:
+ * @refl: A %Reflection
+ * @esd: New standard error for this reflection's intensity measurement
+ *
+ **/
+void set_esd_intensity(Reflection *refl, double esd)
+{
+	refl->data.esd_i = esd;
+}
+
+
+/**
+ * set_ph:
+ * @refl: A %Reflection
+ * @phase: New phase for the reflection
+ *
+ **/
+void set_ph(Reflection *refl, double phase)
+{
+	refl->data.phase = phase;
+	refl->data.have_phase = 1;
+}
+
+
+/**
+ * set_symmetric_indices:
+ * @refl: A %Reflection
+ * @hs: The 'h' index of the reflection
+ * @ks: The 'k' index of the reflection
+ * @ls: The 'l' index of the reflection
+ *
+ * This function gives the symmetric indices, that is, the "real" indices before
+ * squashing down to the asymmetric reciprocal unit.  This may be useful if the
+ * list is indexed according to the asymmetric indices, but you still need
+ * access to the symmetric version.  This happens during post-refinement.
+ *
+ **/
+void set_symmetric_indices(Reflection *refl,
+                           signed int hs, signed int ks, signed int ls)
+{
+	refl->data.hs = hs;
+	refl->data.ks = ks;
+	refl->data.ls = ls;
+}
+
+
+/**
+ * set_temp1
+ * @refl: A %Reflection
+ * @temp: New temporary value for the reflection
+ *
+ * The temporary values can be used according to the needs of the calling
+ * program.
+ *
+ **/
+void set_temp1(Reflection *refl, double temp)
+{
+	refl->data.temp1 = temp;
+}
+
+
+/**
+ * set_temp2
+ * @refl: A %Reflection
+ * @temp: New temporary value for the reflection
+ *
+ * The temporary values can be used according to the needs of the calling
+ * program.
+ *
+ **/
+void set_temp2(Reflection *refl, double temp)
+{
+	refl->data.temp2 = temp;
+}
+
+
+/********************************* Insertion **********************************/
+
+static Reflection *rotate_once(Reflection *refl, int dir)
+{
+	Reflection *s = refl->child[!dir];
+
+	refl->child[!dir] = s->child[dir];
+	s->child[dir] = refl;
+
+	refl->col = RED;
+	s->col = BLACK;
+
+	return s;
+}
+
+
+static Reflection *rotate_twice(Reflection *refl, int dir)
+{
+	refl->child[!dir] = rotate_once(refl->child[!dir], !dir);
+	return rotate_once(refl, dir);
+}
+
+
+static int is_red(Reflection *refl)
+{
+	return (refl != NULL) && (refl->col == RED);
+}
+
+
+static Reflection *insert_node(Reflection *refl, Reflection *new)
+{
+	if ( refl == NULL ) {
+
+		refl = new;
+
+	} else {
+
+		int dir;
+		Reflection *rcd;
+
+		assert(new->serial != refl->serial);
+		dir = new->serial > refl->serial;
+		refl->child[dir] = insert_node(refl->child[dir], new);
+
+		rcd = refl->child[dir];
+		if ( is_red(rcd) ) {
+
+			if ( is_red(refl->child[!dir]) ) {
+
+				refl->col = RED;
+				refl->child[0]->col = BLACK;
+				refl->child[1]->col = BLACK;
+
+			} else {
+
+				if ( is_red(rcd->child[dir] ) ) {
+					refl = rotate_once(refl, !dir);
+				} else if ( is_red(rcd->child[!dir] ) ) {
+					refl = rotate_twice(refl, !dir);
+				}
+
+			}
+		}
+
+	}
+
+	return refl;
+}
+
+
+/**
+ * add_refl
+ * @list: A %RefList
+ * @h: The 'h' index of the reflection
+ * @k: The 'k' index of the reflection
+ * @l: The 'l' index of the reflection
+ *
+ * Adds a new reflection to @list.  Note that the implementation allows there to
+ * be multiple reflections with the same indices in the list, so this function
+ * should succeed even if the given indices already feature in the list.
+ *
+ * Returns: The newly created reflection, or NULL on failure.
+ *
+ **/
+Reflection *add_refl(RefList *list, signed int h, signed int k, signed int l)
+{
+	Reflection *new;
+	Reflection *f;
+
+	assert(abs(h)<256);
+	assert(abs(k)<256);
+	assert(abs(l)<256);
+
+	new = new_node(SERIAL(h, k, l));
+	if ( new == NULL ) return NULL;
+
+	f = find_refl(list, h, k, l);
+	if ( f == NULL ) {
+
+		list->head = insert_node(list->head, new);
+		list->head->col = BLACK;
+
+	} else {
+
+		/* New reflection is identical to a previous one */
+		while ( f->next != NULL ) {
+			f = f->next;
+		}
+		f->next = new;
+		new->prev = f;
+	}
+
+	return new;
+}
+
+
+/**
+ * add_refl_to_list
+ * @refl: A %Reflection
+ * @list: A %RefList
+ *
+ * Adds a reflection to @list.  The reflection that actually gets added will be
+ * a newly created one, and all the data will be copied across.  The original
+ * reflection will be destroyed and the new reflection returned.
+ *
+ * Returns: The newly created reflection, or NULL on failure.
+ *
+ **/
+Reflection *add_refl_to_list(Reflection *refl, RefList *list)
+{
+	signed int h, k, l;
+	Reflection *r_added;
+
+	get_indices(refl, &h, &k, &l);
+	r_added = add_refl(list, h, k, l);
+	if ( r_added == NULL ) return NULL;
+
+	copy_data(r_added, refl);
+	reflection_free(refl);
+
+	return r_added;
+}
+
+
+/********************************* Iteration **********************************/
+
+struct _reflistiterator {
+
+	int stack_size;
+	int stack_ptr;
+	Reflection **stack;
+
+};
+
+
+/**
+ * first_refl:
+ * @list: A %RefList to iterate over
+ * @piter: Address at which to store a %RefListIterator
+ *
+ * This function sets up the state required for iteration over the entire list,
+ * and then returns the first reflection in the list.  An iterator object will
+ * be created and its address stored at the location given in piter.
+ *
+ * Returns: the first reflection in the list.
+ *
+ **/
+Reflection *first_refl(RefList *list, RefListIterator **piter)
+{
+	RefListIterator *iter;
+
+	iter = malloc(sizeof(struct _reflistiterator));
+	iter->stack_size = 32;
+	iter->stack = malloc(iter->stack_size*sizeof(Reflection *));
+	iter->stack_ptr = 0;
+	*piter = iter;
+
+	Reflection *refl = list->head;
+
+	do {
+
+		if ( refl != NULL ) {
+			iter->stack[iter->stack_ptr++] = refl;
+			if ( iter->stack_ptr == iter->stack_size ) {
+				iter->stack_size += 32;
+				iter->stack = realloc(iter->stack,
+				         iter->stack_size*sizeof(Reflection *));
+			}
+			refl = refl->child[0];
+			continue;
+		}
+
+		if ( iter->stack_ptr == 0 ) {
+			free(iter->stack);
+			free(iter);
+			return NULL;
+		}
+
+		refl = iter->stack[--iter->stack_ptr];
+
+		return refl;
+
+	} while ( 1 );
+}
+
+
+/**
+ * next_refl:
+ * @refl: A reflection
+ * @iter: A %RefListIterator
+ *
+ * This function looks up the next reflection in the list that was given earlier
+ * to first_refl().
+ *
+ * Returns: the next reflection in the list, or NULL if no more.
+ *
+ **/
+Reflection *next_refl(Reflection *refl, RefListIterator *iter)
+{
+	int returned = 1;
+
+	do {
+
+		if ( returned ) refl = refl->child[1];
+		returned = 0;
+
+		if ( refl != NULL ) {
+
+			iter->stack[iter->stack_ptr++] = refl;
+			if ( iter->stack_ptr == iter->stack_size ) {
+				iter->stack_size += 32;
+				iter->stack = realloc(iter->stack,
+				         iter->stack_size*sizeof(Reflection *));
+			}
+			refl = refl->child[0];
+			continue;
+
+		}
+		if ( iter->stack_ptr == 0 ) {
+			free(iter->stack);
+			free(iter);
+			return NULL;
+		}
+
+		return iter->stack[--iter->stack_ptr];
+
+	} while ( 1 );
+}
+
+
+/*********************************** Voodoo ***********************************/
+
+static int recursive_depth(Reflection *refl)
+{
+	int depth_left, depth_right;
+
+	if ( refl == NULL ) return 0;
+
+	depth_left = recursive_depth(refl->child[0]);
+	depth_right = recursive_depth(refl->child[1]);
+
+	return 1 + biggest(depth_left, depth_right);
+}
+
+
+static int recursive_count(Reflection *refl)
+{
+	int count_left, count_right;
+
+	if ( refl == NULL ) return 0;
+
+	count_left = recursive_count(refl->child[0]);
+	count_right = recursive_count(refl->child[1]);
+
+	return 1 + count_left + count_right;
+}
+
+
+/**
+ * num_reflections:
+ * @list: A %RefList
+ *
+ * Returns: the number of reflections in @list.
+ *
+ **/
+int num_reflections(RefList *list)
+{
+	return recursive_count(list->head);
+}
+
+
+/**
+ * tree_depth:
+ * @list: A %RefList
+ *
+ * If the depth of the tree is more than about 20, access to the list will be
+ * slow.  This should never happen.
+ *
+ * Returns: the depth of the RB-tree used internally to represent @list.
+ *
+ **/
+int tree_depth(RefList *list)
+{
+	return recursive_depth(list->head);
+}
+
+
+/**
+ * lock_reflection:
+ * @refl: A %Reflection
+ *
+ * Acquires a lock on the reflection.
+ */
+void lock_reflection(Reflection *refl)
+{
+	pthread_mutex_lock(&refl->lock);
+}
+
+
+/**
+ * unlock_reflection:
+ * @refl: A %Reflection
+ *
+ * Releases a lock on the reflection.
+ */
+void unlock_reflection(Reflection *refl)
+{
+	pthread_mutex_unlock(&refl->lock);
+}
diff --git a/libcrystfel/src/reflist.h b/libcrystfel/src/reflist.h
new file mode 100644
index 00000000..955db69b
--- /dev/null
+++ b/libcrystfel/src/reflist.h
@@ -0,0 +1,112 @@
+/*
+ * reflist.h
+ *
+ * Fast reflection/peak list
+ *
+ * (c) 2011 Thomas White <taw@physics.org>
+ *
+ * Part of CrystFEL - crystallography with a FEL
+ *
+ */
+
+#ifndef REFLIST_H
+#define REFLIST_H
+
+#ifdef HAVE_CONFIG_H
+#include <config.h>
+#endif
+
+/**
+ * RefList:
+ *
+ * A %RefList represents a list of Bragg reflections.
+ *
+ * This data structure is opaque.  You must use the available accessor functions
+ * to read and write its contents.
+ *
+ **/
+typedef struct _reflist RefList;
+
+/**
+ * Reflection:
+ *
+ * A %Reflection represents a single Bragg reflection.
+ *
+ * This data structure is opaque.  You must use the available accessor functions
+ * to read and write its contents.
+ *
+ **/
+typedef struct _reflection Reflection;
+
+/**
+ * RefListIterator:
+ *
+ * A %RefListIterator is an opaque data type used when iterating over a
+ * %RefList.
+ *
+ **/
+typedef struct _reflistiterator RefListIterator;
+
+/* Creation/deletion */
+extern RefList *reflist_new(void);
+extern void reflist_free(RefList *list);
+extern Reflection *reflection_new(signed int h, signed int k, signed int l);
+extern void reflection_free(Reflection *refl);
+
+/* Search */
+extern Reflection *find_refl(const RefList *list, signed int h, signed int k, signed int l);
+extern Reflection *next_found_refl(Reflection *refl);
+
+/* Get */
+extern double get_excitation_error(const Reflection *refl);
+extern void get_detector_pos(const Reflection *refl, double *fs, double *ss);
+extern double get_partiality(const Reflection *refl);
+extern void get_indices(const Reflection *refl,
+                        signed int *h, signed int *k, signed int *l);
+extern void get_symmetric_indices(const Reflection *refl,
+                                  signed int *hs, signed int *ks,
+                                  signed int *ls);
+extern double get_intensity(const Reflection *refl);
+extern void get_partial(const Reflection *refl, double *r1, double *r2,
+                        double *p, int *clamp_low, int *clamp_high);
+extern int get_scalable(const Reflection *refl);
+extern int get_refinable(const Reflection *refl);
+extern int get_redundancy(const Reflection *refl);
+extern double get_temp1(const Reflection *refl);
+extern double get_temp2(const Reflection *refl);
+extern double get_esd_intensity(const Reflection *refl);
+extern double get_phase(const Reflection *refl, int *have_phase);
+
+/* Set */
+extern void copy_data(Reflection *to, const Reflection *from);
+extern void set_detector_pos(Reflection *refl, double exerr,
+                             double fs, double ss);
+extern void set_partial(Reflection *refl, double r1, double r2, double p,
+                        double clamp_low, double clamp_high);
+extern void set_int(Reflection *refl, double intensity);
+extern void set_scalable(Reflection *refl, int scalable);
+extern void set_refinable(Reflection *refl, int refinable);
+extern void set_redundancy(Reflection *refl, int red);
+extern void set_temp1(Reflection *refl, double temp);
+extern void set_temp2(Reflection *refl, double temp);
+extern void set_esd_intensity(Reflection *refl, double esd);
+extern void set_ph(Reflection *refl, double phase);
+extern void set_symmetric_indices(Reflection *refl,
+                                  signed int hs, signed int ks, signed int ls);
+
+/* Insertion */
+extern Reflection *add_refl(RefList *list,
+                            signed int h, signed int k, signed int l);
+extern Reflection *add_refl_to_list(Reflection *refl, RefList *list);
+
+/* Iteration */
+extern Reflection *first_refl(RefList *list, RefListIterator **piter);
+extern Reflection *next_refl(Reflection *refl, RefListIterator *iter);
+
+/* Misc */
+extern int num_reflections(RefList *list);
+extern int tree_depth(RefList *list);
+extern void lock_reflection(Reflection *refl);
+extern void unlock_reflection(Reflection *refl);
+
+#endif	/* REFLIST_H */
diff --git a/libcrystfel/src/thread-pool.c b/libcrystfel/src/thread-pool.c
new file mode 100644
index 00000000..0ae26e17
--- /dev/null
+++ b/libcrystfel/src/thread-pool.c
@@ -0,0 +1,283 @@
+/*
+ * thread-pool.c
+ *
+ * A thread pool implementation
+ *
+ * (c) 2006-2011 Thomas White <taw@physics.org>
+ *
+ * Part of CrystFEL - crystallography with a FEL
+ *
+ */
+
+
+#ifdef HAVE_CONFIG_H
+#include <config.h>
+#endif
+
+
+#include <stdarg.h>
+#include <stdlib.h>
+#include <stdio.h>
+#include <string.h>
+#include <unistd.h>
+#include <pthread.h>
+#include <assert.h>
+
+#ifdef HAVE_CPU_AFFINITY
+#include <sched.h>
+#endif
+
+
+#include "utils.h"
+
+
+/**
+ * SECTION:thread-pool
+ * @short_description: The thread pool
+ * @title: The thread pool
+ * @section_id:
+ * @see_also:
+ * @include: "thread-pool.h"
+ * @Image:
+ *
+ * The thread pool helps when running many tasks in parallel.  It takes care of
+ * starting and stopping threads, and presents a relatively simple interface to
+ * the individual programs.
+ */
+
+/* ------------------------------ CPU affinity ------------------------------ */
+
+#ifdef HAVE_CPU_AFFINITY
+
+static void set_affinity(int n, int cpu_num, int cpu_groupsize, int cpu_offset)
+{
+	cpu_set_t c;
+	int group;
+	int n_cpu_groups;
+	int i;
+
+	if ( cpu_num == 0 ) return;
+
+	CPU_ZERO(&c);
+
+	/* Work out which group this thread belongs to */
+	group = (n / cpu_groupsize) + cpu_offset;
+
+	/* Work out which CPUs should be used for this group */
+	n_cpu_groups = cpu_num / cpu_groupsize;
+	group = group % n_cpu_groups;
+
+	/* Set flags */
+	for ( i=0; i<cpu_groupsize; i++ ) {
+
+		int cpu = cpu_groupsize*group + i;
+
+		CPU_SET(cpu, &c);
+
+	}
+
+	if ( sched_setaffinity(0, sizeof(cpu_set_t), &c) ) {
+
+		/* Cannot use ERROR() just yet */
+		fprintf(stderr, "%i: Failed to set CPU affinity.\n", n);
+
+	}
+}
+
+#else /* HAVE_CPU_AFFINITY */
+
+static void set_affinity(int n, int cpu_num, int cpu_groupsize, int cpu_offset)
+{
+	/* Do absolutely nothing */
+}
+
+#endif /* HAVE_CPU_AFFINITY */
+
+
+/* --------------------------- Status label stuff --------------------------- */
+
+static int use_status_labels = 0;
+static pthread_key_t status_label_key;
+pthread_mutex_t stderr_lock = PTHREAD_MUTEX_INITIALIZER;
+
+struct worker_args
+{
+	struct task_queue_range *tqr;
+	struct task_queue *tq;
+	int id;
+	int cpu_num;
+	int cpu_groupsize;
+	int cpu_offset;
+};
+
+
+signed int get_status_label()
+{
+	int *cookie;
+
+	if ( !use_status_labels ) {
+		return -1;
+	}
+
+	cookie = pthread_getspecific(status_label_key);
+	return *cookie;
+}
+
+
+struct task_queue
+{
+	pthread_mutex_t  lock;
+
+	int              n_started;
+	int              n_completed;
+	int              max;
+
+	void *(*get_task)(void *);
+	void (*finalise)(void *, void *);
+	void *queue_args;
+	void (*work)(void *, int);
+};
+
+
+static void *task_worker(void *pargsv)
+{
+	struct worker_args *w = pargsv;
+	struct task_queue *q = w->tq;
+	int *cookie;
+
+	set_affinity(w->id, w->cpu_num, w->cpu_groupsize, w->cpu_offset);
+
+	cookie = malloc(sizeof(int));
+	*cookie = w->id;
+	pthread_setspecific(status_label_key, cookie);
+
+	free(w);
+
+	do {
+
+		void *task;
+		int cookie;
+
+		/* Get a task */
+		pthread_mutex_lock(&q->lock);
+		if ( (q->max) && (q->n_started >= q->max) ) {
+			pthread_mutex_unlock(&q->lock);
+			break;
+		}
+		task = q->get_task(q->queue_args);
+
+		/* No more tasks? */
+		if ( task == NULL ) {
+			pthread_mutex_unlock(&q->lock);
+			break;
+		}
+
+		q->n_started++;
+		pthread_mutex_unlock(&q->lock);
+
+		cookie = *(int *)pthread_getspecific(status_label_key);
+		q->work(task, cookie);
+
+		/* Update totals etc */
+		pthread_mutex_lock(&q->lock);
+		q->n_completed++;
+		if ( q->finalise ) {
+			q->finalise(q->queue_args, task);
+		}
+		pthread_mutex_unlock(&q->lock);
+
+	} while ( 1 );
+
+	free(cookie);
+
+	return NULL;
+}
+
+
+/**
+ * run_threads:
+ * @n_threads: The number of threads to run in parallel
+ * @work: The function to be called to do the work
+ * @get_task: The function which will determine the next unassigned task
+ * @final: The function which will be called to clean up after a task
+ * @queue_args: A pointer to any data required to determine the next task
+ * @max: Stop calling get_task after starting this number of jobs
+ * @cpu_num: The number of CPUs in the system
+ * @cpu_groupsize: The group size into which the CPUs are grouped
+ * @cpu_offset: The CPU group number at which to start pinning threads
+ *
+ * 'get_task' will be called every time a worker is idle.  It returns either
+ * NULL, indicating that no further work is available, or a pointer which will
+ * be passed to 'work'.
+ *
+ * 'final' will be called once per image, and will be given both queue_args
+ * and the last task pointer.
+ *
+ * 'get_task' and 'final' will be called only under lock, and so do NOT need to
+ * be re-entrant or otherwise thread safe.  'work', of course, needs to be
+ * thread safe.
+ *
+ * Work will stop after 'max' tasks have been processed whether get_task
+ * returned NULL or not.  If "max" is zero, all tasks will be processed.
+ *
+ * Returns: The number of tasks completed.
+ **/
+int run_threads(int n_threads, TPWorkFunc work,
+                TPGetTaskFunc get_task, TPFinalFunc final,
+                void *queue_args, int max,
+                int cpu_num, int cpu_groupsize, int cpu_offset)
+{
+	pthread_t *workers;
+	int i;
+	struct task_queue q;
+
+	pthread_key_create(&status_label_key, NULL);
+
+	workers = malloc(n_threads * sizeof(pthread_t));
+
+	pthread_mutex_init(&q.lock, NULL);
+	q.work = work;
+	q.get_task = get_task;
+	q.finalise = final;
+	q.queue_args = queue_args;
+	q.n_started = 0;
+	q.n_completed = 0;
+	q.max = max;
+
+	/* Now it's safe to start using the status labels */
+	if ( n_threads > 1 ) use_status_labels = 1;
+
+	/* Start threads */
+	for ( i=0; i<n_threads; i++ ) {
+
+		struct worker_args *w;
+
+		w = malloc(sizeof(struct worker_args));
+
+		w->tq = &q;
+		w->tqr = NULL;
+		w->id = i;
+		w->cpu_num = cpu_num;
+		w->cpu_groupsize = cpu_groupsize;
+		w->cpu_offset = cpu_offset;
+
+		if ( pthread_create(&workers[i], NULL, task_worker, w) ) {
+			/* Not ERROR() here */
+			fprintf(stderr, "Couldn't start thread %i\n", i);
+			n_threads = i;
+			break;
+		}
+
+	}
+
+	/* Join threads */
+	for ( i=0; i<n_threads; i++ ) {
+		pthread_join(workers[i], NULL);
+	}
+
+	use_status_labels = 0;
+
+	free(workers);
+
+	return q.n_completed;
+}
diff --git a/libcrystfel/src/thread-pool.h b/libcrystfel/src/thread-pool.h
new file mode 100644
index 00000000..a99a7ade
--- /dev/null
+++ b/libcrystfel/src/thread-pool.h
@@ -0,0 +1,71 @@
+/*
+ * thread-pool.h
+ *
+ * A thread pool implementation
+ *
+ * (c) 2006-2010 Thomas White <taw@physics.org>
+ *
+ * Part of CrystFEL - crystallography with a FEL
+ *
+ */
+
+
+#ifndef THREAD_POOL_H
+#define THREAD_POOL_H
+
+#ifdef HAVE_CONFIG_H
+#include <config.h>
+#endif
+
+
+#include <pthread.h>
+
+extern pthread_mutex_t stderr_lock;
+extern signed int get_status_label(void);
+
+
+/**
+ * TPGetTaskFunc:
+ * @qargs: The queue_args pointer which was given to run_threads().
+ * Returns: A pointer which will be passed to the worker function.
+ *
+ * This function is called, non-reentrantly, to get a new work item to give to
+ * your work function.  The stuff you need to generate the new work item should
+ * have been stored in @qargs which was passed to run_threads().
+ *
+ **/
+typedef void *(*TPGetTaskFunc)(void *qargs);
+
+
+/**
+ * TPWorkFunc:
+ * @work: The queue_args pointer which was given to run_threads().
+ * @cookie: A small integral number which is guaranteed to be unique among all
+ * currently running threads.
+ *
+ * This function is called, reentrantly, for each work item.
+ *
+ **/
+typedef void (*TPWorkFunc)(void *work, int cookie);
+
+
+/**
+ * TPFinalFunc:
+ * @qargs: The queue_args pointer which was given to run_threads().
+ * @work: The pointer which was returned by your get_task function.
+ *
+ * This function is called, non-reentrantly, after each work item has been
+ * completed.  A typical use might be to update some counters inside @qargs
+ * according to fields withing @work which were filled by your 'work' function.
+ *
+ **/
+typedef void (*TPFinalFunc)(void *qargs, void *work);
+
+
+extern int run_threads(int n_threads, TPWorkFunc work,
+                       TPGetTaskFunc get_task, TPFinalFunc final,
+                       void *queue_args, int max,
+                       int cpu_num, int cpu_groupsize, int cpu_offset);
+
+
+#endif	/* THREAD_POOL_H */
diff --git a/libcrystfel/src/utils.c b/libcrystfel/src/utils.c
new file mode 100644
index 00000000..6b29627f
--- /dev/null
+++ b/libcrystfel/src/utils.c
@@ -0,0 +1,488 @@
+/*
+ * utils.c
+ *
+ * Utility stuff
+ *
+ * (c) 2006-2010 Thomas White <taw@physics.org>
+ *
+ * Part of CrystFEL - crystallography with a FEL
+ *
+ */
+
+#include <libgen.h>
+#include <math.h>
+#include <string.h>
+#include <stdio.h>
+#include <unistd.h>
+#include <sys/types.h>
+#include <sys/stat.h>
+#include <gsl/gsl_matrix.h>
+#include <gsl/gsl_vector.h>
+#include <gsl/gsl_blas.h>
+
+#include "utils.h"
+#include "image.h"
+
+
+/**
+ * SECTION:utils
+ * @short_description: Miscellaneous utilities
+ * @title: Utilities
+ * @section_id:
+ * @see_also:
+ * @include: "utils.h"
+ * @Image:
+ *
+ * Wibble
+ */
+
+/**
+ * show_matrix_eqn:
+ * @M: A matrix
+ * @v: A vector
+ * @r: The number of elements in @v and the side length of @M.
+ *
+ * Displays a matrix equation of the form @M.a = @v.
+ *
+ * @M must be square.
+ **/
+void show_matrix_eqn(gsl_matrix *M, gsl_vector *v, int r)
+{
+	int i, j;
+
+	for ( i=0; i<r; i++ ) {
+		STATUS("[ ");
+		for ( j=0; j<r; j++ ) {
+			STATUS("%+9.3e ", gsl_matrix_get(M, i, j));
+		}
+		STATUS("][ a%2i ] = [ %+9.3e ]\n", i, gsl_vector_get(v, i));
+	}
+}
+
+
+size_t notrail(char *s)
+{
+	size_t i;
+	size_t munched = 0;
+
+	for ( i=strlen(s)-1; i>=0; i-- ) {
+		if ( (s[i] == ' ') || (s[i] == '\t') ) {
+			s[i] = '\0';
+			munched++;
+		} else {
+			return munched;
+		}
+	}
+
+	return munched;
+}
+
+
+void chomp(char *s)
+{
+	size_t i;
+
+	if ( !s ) return;
+
+	for ( i=0; i<strlen(s); i++ ) {
+		if ( (s[i] == '\n') || (s[i] == '\r') ) {
+			s[i] = '\0';
+			return;
+		}
+	}
+}
+
+
+void progress_bar(int val, int total, const char *text)
+{
+	double frac;
+	int n, i;
+	char s[1024];
+	const int width = 50;
+
+	if ( total == 0 ) return;
+
+	if ( !isatty(STDERR_FILENO) ) return;
+	if ( tcgetpgrp(STDERR_FILENO) != getpgrp() ) return;
+
+	frac = (double)val/total;
+	n = (int)(frac*width);
+
+	for ( i=0; i<n; i++ ) s[i] = '=';
+	for ( i=n; i<width; i++ ) s[i] = '.';
+	s[width] = '\0';
+
+	pthread_mutex_lock(&stderr_lock);
+	fprintf(stderr, "\r%s: |%s|", text, s);
+	if ( val == total ) fprintf(stderr, "\n");
+	pthread_mutex_unlock(&stderr_lock);
+
+	fflush(stdout);
+}
+
+
+double random_flat(double max)
+{
+	return max * (double)random()/RAND_MAX;
+}
+
+
+double flat_noise(double expected, double width)
+{
+	double noise = random_flat(2.0*width);
+	return expected+noise-width;
+}
+
+
+double gaussian_noise(double expected, double stddev)
+{
+	double x1, x2, noise;
+
+	/* A uniformly distributed random number between 0 and 1 */
+	x1 = ((double)random()/RAND_MAX);
+	x2 = ((double)random()/RAND_MAX);
+
+	noise = sqrt(-2.0*log(x1)) * cos(2.0*M_PI*x2);
+
+	return expected + noise*stddev;
+}
+
+
+static int fake_poisson_noise(double expected)
+{
+	double rf = gaussian_noise(expected, sqrt(expected));
+	return (int)rf;
+}
+
+
+int poisson_noise(double expected)
+{
+	double L;
+	int k = 0;
+	double p = 1.0;
+
+	/* For large values of the mean, we get big problems with arithmetic.
+	 * In such cases, fall back on a Gaussian with the right variance. */
+	if ( expected > 100.0 ) return fake_poisson_noise(expected);
+
+	L = exp(-expected);
+
+	do {
+
+		double r;
+
+		k++;
+		r = (double)random()/(double)RAND_MAX;
+		p *= r;
+
+	} while ( p > L );
+
+	return k - 1;
+}
+
+
+/**
+ * SECTION:quaternion
+ * @short_description: Simple quaternion handling
+ * @title: Quaternion
+ * @section_id:
+ * @see_also:
+ * @include: "utils.h"
+ * @Image:
+ *
+ * There is a simple quaternion structure in CrystFEL.  At the moment, it is
+ * only used when simulating patterns, as an argument to cell_rotate() to
+ * orient the unit cell.
+ */
+
+/**
+ * quaternion_modulus:
+ * @q: A %quaternion
+ *
+ * If a quaternion represents a pure rotation, its modulus should be unity.
+ *
+ * Returns: the modulus of the given quaternion.
+ **/
+double quaternion_modulus(struct quaternion q)
+{
+	return sqrt(q.w*q.w + q.x*q.x + q.y*q.y + q.z*q.z);
+}
+
+
+/**
+ * normalise_quaternion:
+ * @q: A %quaternion
+ *
+ * Rescales the quaternion such that its modulus is unity.
+ *
+ * Returns: the normalised version of @q
+ **/
+struct quaternion normalise_quaternion(struct quaternion q)
+{
+	double mod;
+	struct quaternion r;
+
+	mod = quaternion_modulus(q);
+
+	r.w = q.w / mod;
+	r.x = q.x / mod;
+	r.y = q.y / mod;
+	r.z = q.z / mod;
+
+	return r;
+}
+
+
+/**
+ * random_quaternion:
+ *
+ * Returns: a randomly generated, normalised, quaternion.
+ **/
+struct quaternion random_quaternion()
+{
+	struct quaternion q;
+
+	q.w = 2.0*(double)random()/RAND_MAX - 1.0;
+	q.x = 2.0*(double)random()/RAND_MAX - 1.0;
+	q.y = 2.0*(double)random()/RAND_MAX - 1.0;
+	q.z = 2.0*(double)random()/RAND_MAX - 1.0;
+	q = normalise_quaternion(q);
+
+	return q;
+}
+
+
+/**
+ * quaternion_valid:
+ * @q: A %quaternion
+ *
+ * Checks if the given quaternion is normalised.
+ *
+ * This function performs a nasty floating point comparison of the form
+ * <code>(modulus > 0.999) && (modulus < 1.001)</code>, and so should not be
+ * relied upon to spot anything other than the most obvious input error.
+ *
+ * Returns: 1 if the quaternion is normalised, 0 if not.
+ **/
+int quaternion_valid(struct quaternion q)
+{
+	double qmod;
+
+	qmod = quaternion_modulus(q);
+
+	/* Modulus = 1 to within some tolerance?
+	 * Nasty allowance for floating-point accuracy follows... */
+	if ( (qmod > 0.999) && (qmod < 1.001) ) return 1;
+
+	return 0;
+}
+
+
+/**
+ * quat_rot
+ * @q: A vector (in the form of an %rvec)
+ * @z: A %quaternion
+ *
+ * Rotates a vector according to a quaternion.
+ *
+ * Returns: A rotated version of @p.
+ **/
+struct rvec quat_rot(struct rvec q, struct quaternion z)
+{
+	struct rvec res;
+	double t01, t02, t03, t11, t12, t13, t22, t23, t33;
+
+	t01 = z.w*z.x;
+	t02 = z.w*z.y;
+	t03 = z.w*z.z;
+	t11 = z.x*z.x;
+	t12 = z.x*z.y;
+	t13 = z.x*z.z;
+	t22 = z.y*z.y;
+	t23 = z.y*z.z;
+	t33 = z.z*z.z;
+
+	res.u = (1.0 - 2.0 * (t22 + t33)) * q.u
+	            + (2.0 * (t12 + t03)) * q.v
+	            + (2.0 * (t13 - t02)) * q.w;
+
+	res.v =       (2.0 * (t12 - t03)) * q.u
+	      + (1.0 - 2.0 * (t11 + t33)) * q.v
+	            + (2.0 * (t01 + t23)) * q.w;
+
+	res.w =       (2.0 * (t02 + t13)) * q.u
+	            + (2.0 * (t23 - t01)) * q.v
+	      + (1.0 - 2.0 * (t11 + t22)) * q.w;
+
+	return res;
+}
+
+
+/* Return non-zero if c is in delims */
+static int assplode_isdelim(const char c, const char *delims)
+{
+	size_t i;
+	for ( i=0; i<strlen(delims); i++ ) {
+		if ( c == delims[i] ) return 1;
+	}
+	return 0;
+}
+
+
+static int assplode_extract(char ***pbits, int n, size_t n_captured,
+                            size_t start, const char *a)
+{
+	char **bits = *pbits;
+	bits = realloc(bits, sizeof(char *)*(n+1));
+	bits[n] = malloc(n_captured+1);
+	memcpy(bits[n], a+start, n_captured);
+	bits[n][n_captured] = '\0';
+	n++;
+	*pbits = bits;
+	return n;
+}
+
+
+/* Split the string 'a' using 'delims' as a zero-terminated list of
+ *  deliminators.
+ * Store each segment in bits[0...n] where n is the number of segments and is
+ *  the return value.  pbits = &bits
+ * Each segment needs to be freed with free() when finished with.
+ * The array of bits also needs to be freed with free() when finished with,
+ *  unless n=0 in which case bits==NULL
+ */
+int assplode(const char *a, const char *delims, char ***pbits,
+             AssplodeFlag flags)
+{
+	size_t i, start, n_captured;
+	int n, last_was_delim;
+	char **bits;
+
+	n = 0;
+	i = 0;
+	n_captured = 0;
+	start = 0;
+	last_was_delim = 0;
+	bits = NULL;
+	while ( i < strlen(a) ) {
+
+		if ( assplode_isdelim(a[i], delims) ) {
+
+			if ( n_captured > 0 ) {
+				/* This is a deliminator after a sequence of
+				 * non-deliminator chars */
+				n = assplode_extract(&bits, n, n_captured,
+				                     start, a);
+			}
+
+			n_captured = 0;
+			if ( (flags & ASSPLODE_DUPS) && last_was_delim ) {
+				n = assplode_extract(&bits, n, 0, start, a);
+			}
+			last_was_delim = 1;
+
+		} else {
+
+			if ( n_captured == 0 ) {
+				/* No characters currently found, so this is the
+				 * start */
+				start = i;
+			}
+			n_captured++;
+			last_was_delim = 0;
+
+		}
+
+		i++;
+
+	}
+	/* Left over characters at the end? */
+	if ( n_captured > 0 ) {
+		n = assplode_extract(&bits, n, n_captured, start, a);
+	}
+
+	*pbits = bits;
+	return n;
+}
+
+
+char *check_prefix(char *prefix)
+{
+	int r;
+	struct stat statbuf;
+	char *new;
+	size_t len;
+
+	/* Is "prefix" a directory? */
+	r = stat(prefix, &statbuf);
+	if ( r != 0 ) {
+		/* "prefix" probably doesn't exist.  This is fine - assume
+		 * the user knows what they're doing, and that "prefix"
+		 * suffixed with the actual filename will produce something
+		 * sensible. */
+		return prefix;
+	}
+
+	if ( !S_ISDIR(statbuf.st_mode) ) {
+		/* Also fine, as above. */
+		return prefix;
+	}
+
+	/* Does the prefix end in a slash? */
+	if ( prefix[strlen(prefix)-1] == '/' ) {
+		/* This looks sensible. */
+		return prefix;
+	}
+
+	STATUS("Your prefix ('%s') is a directory, but doesn't end"
+	       " with a slash.  I'm going to add it for you.\n", prefix);
+	STATUS("If this isn't what you want, run with --no-check-prefix.\n");
+	len = strlen(prefix)+2;
+	new = malloc(len);
+	snprintf(new, len, "%s/", prefix);
+	free(prefix);
+	return new;
+}
+
+
+char *safe_basename(const char *in)
+{
+	int i;
+	char *cpy;
+	char *res;
+
+	cpy = strdup(in);
+
+	/* Get rid of any trailing slashes */
+	for ( i=strlen(cpy)-1; i>0; i-- ) {
+		if ( cpy[i] == '/' ) {
+			cpy[i] = '\0';
+		} else {
+			break;
+		}
+	}
+
+	/* Find the base name */
+	for ( i=strlen(cpy)-1; i>0; i-- ) {
+		if ( cpy[i] == '/' ) {
+			i++;
+			break;
+		}
+	}
+
+	res = strdup(cpy+i);
+	/* If we didn't find a previous slash, i==0 so res==cpy */
+
+	free(cpy);
+
+	return res;
+}
+
+
+#ifdef GSL_FUDGE
+/* Force the linker to bring in CBLAS to make GSL happy */
+void utils_fudge_gslcblas()
+{
+        STATUS("%p\n", cblas_sgemm);
+}
+#endif
diff --git a/libcrystfel/src/utils.h b/libcrystfel/src/utils.h
new file mode 100644
index 00000000..1179a57f
--- /dev/null
+++ b/libcrystfel/src/utils.h
@@ -0,0 +1,236 @@
+/*
+ * utils.h
+ *
+ * Utility stuff
+ *
+ * (c) 2006-2010 Thomas White <taw@physics.org>
+ *
+ * Part of CrystFEL - crystallography with a FEL
+ *
+ */
+
+#ifndef UTILS_H
+#define UTILS_H
+
+#ifdef HAVE_CONFIG_H
+#include <config.h>
+#endif
+
+#include <math.h>
+#include <complex.h>
+#include <string.h>
+#include <stdlib.h>
+#include <pthread.h>
+#include <gsl/gsl_matrix.h>
+#include <gsl/gsl_vector.h>
+
+
+#include "thread-pool.h"
+
+
+/* -------------------------- Fundamental constants  ------------------------ */
+
+/* Electron charge in C */
+#define ELECTRON_CHARGE (1.6021773e-19)
+
+/* Planck's constant (Js) */
+#define PLANCK (6.62606896e-34)
+
+/* Speed of light in vacuo (m/s) */
+#define C_VACUO (299792458)
+
+/* Thomson scattering length (m) */
+#define THOMSON_LENGTH (2.81794e-15)
+
+
+/* ------------------------------ Quaternions ------------------------------- */
+
+/**
+ * quaternion:
+ *
+ * <programlisting>
+ * struct quaternion
+ * {
+ *    double w
+ *    double x
+ *    double y
+ *    double z
+ * };
+ * </programlisting>
+ *
+ * A structure representing a quaternion.
+ *
+ **/
+struct quaternion;
+
+struct quaternion {
+	double w;
+	double x;
+	double y;
+	double z;
+};
+
+extern struct quaternion normalise_quaternion(struct quaternion q);
+extern double quaternion_modulus(struct quaternion q);
+extern struct quaternion random_quaternion(void);
+extern int quaternion_valid(struct quaternion q);
+extern struct rvec quat_rot(struct rvec q, struct quaternion z);
+
+
+/* --------------------------- Useful functions ----------------------------- */
+
+extern void show_matrix_eqn(gsl_matrix *M, gsl_vector *v, int r);
+extern size_t notrail(char *s);
+extern void chomp(char *s);
+
+typedef enum {
+	ASSPLODE_NONE	= 0,
+	ASSPLODE_DUPS	= 1<<0
+} AssplodeFlag;
+extern int assplode(const char *a, const char *delims, char ***pbits,
+                    AssplodeFlag flags);
+
+extern void progress_bar(int val, int total, const char *text);
+extern double random_flat(double max);
+extern double flat_noise(double expected, double width);
+extern double gaussian_noise(double expected, double stddev);
+extern int poisson_noise(double expected);
+
+/* Keep these ones inline, to avoid function call overhead */
+static inline struct quaternion invalid_quaternion(void)
+{
+	struct quaternion quat;
+	quat.w = 0.0;
+	quat.x = 0.0;
+	quat.y = 0.0;
+	quat.z = 0.0;
+	return quat;
+}
+
+static inline double distance(double x1, double y1, double x2, double y2)
+{
+	return sqrt((x2-x1)*(x2-x1) + (y2-y1)*(y2-y1));
+}
+
+static inline double modulus(double x, double y, double z)
+{
+	return sqrt(x*x + y*y + z*z);
+}
+
+static inline double modulus_squared(double x, double y, double z) {
+	return x*x + y*y + z*z;
+}
+
+static inline double distance3d(double x1, double y1, double z1,
+                                double x2, double y2, double z2)
+{
+	return modulus(x1-x2, y1-y2, z1-z2);
+}
+
+/* Answer in radians */
+static inline double angle_between(double x1, double y1, double z1,
+                                   double x2, double y2, double z2)
+{
+	double mod1 = modulus(x1, y1, z1);
+	double mod2 = modulus(x2, y2, z2);
+	double cosine = (x1*x2 + y1*y2 + z1*z2) / (mod1*mod2);
+
+	/* Fix domain if outside due to rounding */
+	if ( cosine > 1.0 ) cosine = 1.0;
+	if ( cosine < -1.0 ) cosine = -1.0;
+
+	return acos(cosine);
+}
+
+static inline int within_tolerance(double a, double b, double percent)
+{
+	double tol = fabs(a) * (percent/100.0);
+	if ( fabs(b-a) < tol ) return 1;
+	return 0;
+}
+
+
+/* ----------------------------- Useful macros ------------------------------ */
+
+#define rad2deg(a) ((a)*180/M_PI)
+#define deg2rad(a) ((a)*M_PI/180)
+
+#define is_odd(a) ((a)%2==1)
+
+#define biggest(a,b) ((a>b) ? (a) : (b))
+#define smallest(a,b) ((a<b) ? (a) : (b))
+
+
+/* Photon energy (J) to wavelength (m) */
+#define ph_en_to_lambda(a) ((PLANCK*C_VACUO)/(a))
+
+/* Photon wavelength (m) to energy (J) */
+#define ph_lambda_to_en(a) ((PLANCK*C_VACUO)/(a))
+
+/* eV to Joules */
+#define eV_to_J(a) ((a)*ELECTRON_CHARGE)
+
+/* Joules to eV */
+#define J_to_eV(a) ((a)/ELECTRON_CHARGE)
+
+#define UNUSED __attribute__((unused))
+
+/* -------------------- Indexed lists for specified types ------------------- */
+
+#include "defs.h"
+
+#define LIST_SIZE (IDIM*IDIM*IDIM)
+
+/* Create functions for storing reflection intensities indexed as h,k,l */
+#define LABEL(x) x##_intensity
+#define TYPE double
+#include "list_tmp.h"
+
+/* CAs above, but for phase values */
+#define LABEL(x) x##_phase
+#define TYPE double
+#include "list_tmp.h"
+
+/* As above, but for (unsigned) integer counts */
+#define LABEL(x) x##_count
+#define TYPE unsigned int
+#include "list_tmp.h"
+
+/* As above, but for simple flags */
+#define LABEL(x) x##_flag
+#define TYPE unsigned char
+#include "list_tmp.h"
+
+
+/* ------------------------------ Message macros ---------------------------- */
+
+extern pthread_mutex_t stderr_lock;
+
+#define ERROR(...) { \
+                      int error_print_val = get_status_label(); \
+                      pthread_mutex_lock(&stderr_lock); \
+                      if ( error_print_val >= 0 ) { \
+                         fprintf(stderr, "%3i: ", error_print_val); \
+                      } \
+                      fprintf(stderr, __VA_ARGS__); \
+                      pthread_mutex_unlock(&stderr_lock); \
+                   }
+
+#define STATUS(...) { \
+                       int status_print_val = get_status_label(); \
+                       pthread_mutex_lock(&stderr_lock); \
+                       if ( status_print_val >= 0 ) { \
+                          fprintf(stderr, "%3i: ", status_print_val); \
+                       } \
+                       fprintf(stderr, __VA_ARGS__); \
+                       pthread_mutex_unlock(&stderr_lock); \
+                    }
+
+
+/* ------------------------------ File handling ----------------------------- */
+
+extern char *check_prefix(char *prefix);
+extern char *safe_basename(const char *in);
+
+
+#endif	/* UTILS_H */