Merge branch 'tom/symmetry'

29 files changed, 2208 insertions, 639 deletions

diff --git a/Makefile.am b/Makefile.am
index 44e707ea..f0adc640 100644
--- a/Makefile.am
+++ b/Makefile.am
@@ -7,12 +7,13 @@ bin_PROGRAMS = src/pattern_sim src/process_hkl src/get_hkl src/indexamajig \
                src/partialator src/check_hkl src/partial_sim
 
 noinst_PROGRAMS = tests/list_check tests/integration_check \
-                  tests/pr_gradient_check tests/symmetry_check
+                  tests/pr_gradient_check tests/symmetry_check \
+                  tests/centering_check tests/transformation_check
 
 TESTS = tests/list_check tests/first_merge_check tests/second_merge_check \
         tests/third_merge_check tests/fourth_merge_check \
         tests/integration_check tests/pr_gradient_check \
-        tests/symmetry_check
+        tests/symmetry_check tests/centering_check tests/transformation_check
 
 EXTRA_DIST += tests/first_merge_check tests/second_merge_check \
               tests/third_merge_check tests/fourth_merge_check
@@ -80,10 +81,16 @@ tests_integration_check_SOURCES = tests/integration_check.c
 
 tests_symmetry_check_SOURCES = tests/symmetry_check.c
 
-
 tests_pr_gradient_check_SOURCES = tests/pr_gradient_check.c \
                                   src/post-refinement.c
 
+tests_centering_check_SOURCES = tests/centering_check.c libcrystfel/src/cell.c \
+                                libcrystfel/src/cell-utils.c
+
+tests_transformation_check_SOURCES = tests/transformation_check.c \
+                                     libcrystfel/src/cell.c \
+                                     libcrystfel/src/cell-utils.c
+
 INCLUDES = -I$(top_srcdir)/libcrystfel/src -I$(top_srcdir)/data
 
 EXTRA_DIST += src/dw-hdfsee.h src/hdfsee.h src/render_hkl.h \
diff --git a/doc/reference/CrystFEL-docs.sgml b/doc/reference/CrystFEL-docs.sgml
index 89b4832e..e1da1784 100644
--- a/doc/reference/CrystFEL-docs.sgml
+++ b/doc/reference/CrystFEL-docs.sgml
@@ -34,6 +34,7 @@
   <chapter>
     <title>Unit cells</title>
         <xi:include href="xml/unitcell.xml"/>
+        <xi:include href="xml/cell-utils.xml"/>
   </chapter>
 
   <chapter>
diff --git a/doc/reference/CrystFEL-sections.txt b/doc/reference/CrystFEL-sections.txt
index 284740a1..3ccc1511 100644
--- a/doc/reference/CrystFEL-sections.txt
+++ b/doc/reference/CrystFEL-sections.txt
@@ -88,12 +88,28 @@ cell_set_pointgroup
 cell_set_reciprocal
 cell_set_spacegroup
 <SUBSECTION>
+tfn_combine
+tfn_identity
+tfn_inverse
+tfn_print
+tfn_vector
+tfn_free
+</SECTION>
+
+<SECTION>
+<FILE>cell-utils</FILE>
 cell_rotate
 rotate_cell
 cell_print
 resolution
 match_cell
 match_cell_ab
+cell_is_sensible
+validate_cell
+uncenter_cell
+bravais_lattice
+right_handed
+str_lattice
 </SECTION>
 
 <SECTION>
diff --git a/libcrystfel/Makefile.am b/libcrystfel/Makefile.am
index d55f6e74..28d483ac 100644
--- a/libcrystfel/Makefile.am
+++ b/libcrystfel/Makefile.am
@@ -7,7 +7,8 @@ libcrystfel_la_SOURCES = src/reflist.c src/utils.c src/cell.c src/detector.c \
                          src/beam-parameters.c src/geometry.c src/statistics.c \
                          src/symmetry.c src/stream.c src/peaks.c \
                          src/reflist-utils.c src/filters.c \
-                         src/render.c src/index.c src/dirax.c src/mosflm.c
+                         src/render.c src/index.c src/dirax.c src/mosflm.c \
+                         src/cell-utils.c
 
 if HAVE_FFTW
 libcrystfel_la_SOURCES += src/reax.c
@@ -22,7 +23,7 @@ libcrystfel_la_include_HEADERS = src/beam-parameters.h src/hdf5-file.h \
                                  src/geometry.h src/peaks.h src/stream.h \
                                  src/render.h src/index.h src/image.h \
                                  src/filters.h src/dirax.h src/mosflm.h \
-                                 src/index-priv.h src/reax.h
+                                 src/index-priv.h src/reax.h src/cell-utils.h
 
 INCLUDES = "-I$(top_srcdir)/data"
 AM_CPPFLAGS = -DDATADIR=\""$(datadir)"\" -I$(top_builddir)/lib -Wall
diff --git a/libcrystfel/src/cell-utils.c b/libcrystfel/src/cell-utils.c
new file mode 100644
index 00000000..3aaae90e
--- /dev/null
+++ b/libcrystfel/src/cell-utils.c
@@ -0,0 +1,1205 @@
+/*
+ * cell-utils.c
+ *
+ * Unit Cell utility functions
+ *
+ * Copyright © 2012 Deutsches Elektronen-Synchrotron DESY,
+ *                  a research centre of the Helmholtz Association.
+ * Copyright © 2012 Lorenzo Galli
+ *
+ * Authors:
+ *   2009-2012 Thomas White <taw@physics.org>
+ *   2012      Lorenzo Galli
+ *
+ * This file is part of CrystFEL.
+ *
+ * CrystFEL is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU General Public License as published by
+ * the Free Software Foundation, either version 3 of the License, or
+ * (at your option) any later version.
+ *
+ * CrystFEL is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU General Public License
+ * along with CrystFEL.  If not, see <http://www.gnu.org/licenses/>.
+ *
+ */
+
+#ifdef HAVE_CONFIG_H
+#include <config.h>
+#endif
+
+#include <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 <assert.h>
+
+#include "cell.h"
+#include "cell-utils.h"
+#include "utils.h"
+#include "image.h"
+
+
+/**
+ * SECTION:cell-utils
+ * @short_description: Unit cell utilities
+ * @title: Unit cell utilities
+ * @section_id:
+ * @see_also:
+ * @include: "cell-utils.h"
+ * @Image:
+ *
+ * There are some utility functions associated with the core %UnitCell.
+ **/
+
+
+/* Weighting factor of lengths relative to angles */
+#define LWEIGHT (10.0e-9)
+
+
+/**
+ * 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;
+}
+
+
+const char *str_lattice(LatticeType l)
+{
+	switch ( l )
+	{
+		case L_TRICLINIC :    return "triclinic";
+		case L_MONOCLINIC :   return "monoclinic";
+		case L_ORTHORHOMBIC : return "orthorhombic";
+		case L_TETRAGONAL :   return "tetragonal";
+		case L_RHOMBOHEDRAL : return "rhombohedral";
+		case L_HEXAGONAL :    return "hexagonal";
+		case L_CUBIC :        return "cubic";
+	}
+
+	return "unknown lattice";
+}
+
+
+int right_handed(UnitCell *cell)
+{
+	double asx, asy, asz;
+	double bsx, bsy, bsz;
+	double csx, csy, csz;
+	struct rvec aCb;
+	double aCb_dot_c;
+	int rh_reciprocal;
+	int rh_direct;
+
+	if ( cell_get_reciprocal(cell, &asx, &asy, &asz,
+	                               &bsx, &bsy, &bsz,
+	                               &csx, &csy, &csz) ) {
+		ERROR("Couldn't get reciprocal cell.\n");
+		return 0;
+	}
+
+	/* "a" cross "b" */
+	aCb.u = asy*bsz - asz*bsy;
+	aCb.v = - (asx*bsz - asz*bsx);
+	aCb.w = asx*bsy - asy*bsx;
+
+	/* "a cross b" dot "c" */
+	aCb_dot_c = aCb.u*csx + aCb.v*csy + aCb.w*csz;
+
+	rh_reciprocal = aCb_dot_c > 0.0;
+
+	if ( cell_get_cartesian(cell, &asx, &asy, &asz,
+	                              &bsx, &bsy, &bsz,
+	                              &csx, &csy, &csz) ) {
+		ERROR("Couldn't get direct cell.\n");
+		return 0;
+	}
+
+	/* "a" cross "b" */
+	aCb.u = asy*bsz - asz*bsy;
+	aCb.v = - (asx*bsz - asz*bsx);
+	aCb.w = asx*bsy - asy*bsx;
+
+	/* "a cross b" dot "c" */
+	aCb_dot_c = aCb.u*csx + aCb.v*csy + aCb.w*csz;
+
+	rh_direct = aCb_dot_c > 0.0;
+
+	assert(rh_reciprocal == rh_direct);
+
+	return rh_direct;
+}
+
+
+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;
+	LatticeType lt;
+	char cen;
+
+	lt = cell_get_lattice_type(cell);
+	cen = cell_get_centering(cell);
+
+	STATUS("%s %c", str_lattice(lt), cen);
+
+	if ( (lt==L_MONOCLINIC) || (lt==L_TETRAGONAL) || ( lt==L_HEXAGONAL)
+	  || ( (lt==L_ORTHORHOMBIC) && (cen=='A') )
+	  || ( (lt==L_ORTHORHOMBIC) && (cen=='B') )
+	  || ( (lt==L_ORTHORHOMBIC) && (cen=='C') ) )
+	{
+		STATUS(", unique axis %c", cell_get_unique_axis(cell));
+	}
+
+	if ( right_handed(cell) ) {
+		STATUS(", right handed");
+	} else {
+		STATUS(", left handed");
+	}
+
+	STATUS(", point group '%s'.\n", cell_get_pointgroup(cell));
+
+	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("Cell representation is %s.\n", cell_rep(cell));
+}
+
+
+int bravais_lattice(UnitCell *cell)
+{
+	LatticeType lattice = cell_get_lattice_type(cell);
+	char centering = cell_get_centering(cell);
+	char ua = cell_get_unique_axis(cell);
+
+	switch ( centering )
+	{
+		case 'P' :
+		return 1;
+
+		case 'A' :
+		case 'B' :
+		case 'C' :
+		if ( lattice == L_MONOCLINIC ) {
+			if ( (ua=='a') && (centering=='A') ) return 1;
+			if ( (ua=='b') && (centering=='B') ) return 1;
+			if ( (ua=='c') && (centering=='C') ) return 1;
+		} else if ( lattice == L_ORTHORHOMBIC) {
+			return 1;
+		}
+		return 0;
+
+		case 'I' :
+		if ( (lattice == L_ORTHORHOMBIC)
+		  || (lattice == L_TETRAGONAL)
+		  || (lattice == L_CUBIC) )
+		{
+			return 1;
+		}
+		return 0;
+
+		case 'F' :
+		if ( (lattice == L_ORTHORHOMBIC) || (lattice == L_CUBIC) ) {
+			return 1;
+		}
+		return 0;
+
+		case 'H' :
+		/* "Hexagonal H" is not a Bravais lattice, but rather something
+		 * invented by the PDB to make life difficult for programmers.
+		 * Accepting it as Bravais seems to be the least painful way to
+		 * handle it correctly. Yuk. */
+		if ( ua != 'c' ) return 0;
+		if ( lattice == L_HEXAGONAL ) return 1;
+		return 0;
+
+		case 'R' :
+		if ( lattice == L_RHOMBOHEDRAL ) return 1;
+		return 0;
+
+		default :
+		return 0;
+	}
+}
+
+
+static UnitCellTransformation *uncentering_transformation(UnitCell *in,
+                                                          char *new_centering,
+                                                          LatticeType *new_latt)
+{
+	UnitCellTransformation *t;
+	const double OT = 1.0/3.0;
+	const double TT = 2.0/3.0;
+	const double H = 0.5;
+	LatticeType lt;
+	char ua, cen;
+
+	lt = cell_get_lattice_type(in);
+	ua = cell_get_unique_axis(in);
+	cen = cell_get_centering(in);
+
+	t = tfn_identity();
+	if ( t == NULL ) return NULL;
+
+	if ( ua == 'a' ) {
+		tfn_combine(t, tfn_vector(0,0,1),
+		               tfn_vector(0,1,0),
+		               tfn_vector(-1,0,0));
+	}
+
+	if ( ua == 'b' ) {
+		tfn_combine(t, tfn_vector(1,0,0),
+		               tfn_vector(0,0,1),
+		               tfn_vector(0,-1,0));
+	}
+
+	switch ( cen ) {
+
+		case 'P' :
+		*new_latt = lt;
+		*new_centering = 'P';
+		break;
+
+		case 'R' :
+		*new_latt = L_RHOMBOHEDRAL;
+		*new_centering = 'R';
+		break;
+
+		case 'I' :
+		tfn_combine(t, tfn_vector(-H,H,H),
+		               tfn_vector(H,-H,H),
+		               tfn_vector(H,H,-H));
+		if ( lt == L_CUBIC ) {
+			*new_latt = L_RHOMBOHEDRAL;
+			*new_centering = 'R';
+		} else {
+			/* Tetragonal or orthorhombic */
+			*new_latt = L_TRICLINIC;
+			*new_centering = 'P';
+		}
+		break;
+
+		case 'F' :
+		tfn_combine(t, tfn_vector(0,H,H),
+		               tfn_vector(H,0,H),
+		               tfn_vector(H,H,0));
+		if ( lt == L_CUBIC ) {
+			*new_latt = L_RHOMBOHEDRAL;
+			*new_centering = 'R';
+		} else {
+			assert(lt == L_ORTHORHOMBIC);
+			*new_latt = L_TRICLINIC;
+			*new_centering = 'P';
+		}
+		break;
+
+		case 'A' :
+		case 'B' :
+		case 'C' :
+		tfn_combine(t, tfn_vector(H,H,0),
+		               tfn_vector(-H,H,0),
+		               tfn_vector(0,0,1));
+		*new_latt = L_MONOCLINIC;
+		*new_centering = 'P';
+		break;
+
+		case 'H' :
+		/* Obverse setting */
+		tfn_combine(t, tfn_vector(TT,OT,OT),
+		               tfn_vector(-OT,OT,OT),
+		               tfn_vector(-OT,-TT,OT));
+		assert(lt == L_HEXAGONAL);
+		*new_latt = L_RHOMBOHEDRAL;
+		*new_centering = 'R';
+		break;
+
+		default :
+		ERROR("Invalid centering '%c'\n", cell_get_centering(in));
+		return NULL;
+
+	}
+
+	/* Reverse the axis permutation, but only if this was not an H->R
+	 * transformation */
+	if ( !((cen=='H') && (*new_latt == L_RHOMBOHEDRAL)) ) {
+		if ( ua == 'a' ) {
+			tfn_combine(t, tfn_vector(0,0,-1),
+				       tfn_vector(0,1,0),
+				       tfn_vector(1,0,0));
+		}
+
+		if ( ua == 'b' ) {
+			tfn_combine(t, tfn_vector(1,0,0),
+				       tfn_vector(0,0,-1),
+				       tfn_vector(0,1,0));
+		}
+	}
+
+	return t;
+}
+
+
+/**
+ * uncenter_cell:
+ * @in: A %UnitCell
+ * @t: Location at which to store the transformation which was used.
+ *
+ * Turns any cell into a primitive one, e.g. for comparison purposes.  The
+ * transformation which was used is stored at @t, which can be NULL if the
+ * transformation is not required.
+ *
+ * Returns: a primitive version of @in in a conventional (unique axis c)
+ * setting.
+ *
+ */
+UnitCell *uncenter_cell(UnitCell *in, UnitCellTransformation **t)
+{
+	UnitCellTransformation *tt;
+	char new_centering;
+	LatticeType new_latt;
+	UnitCell *out;
+
+	if ( !bravais_lattice(in) ) {
+		ERROR("Cannot uncenter: not a Bravais lattice.\n");
+		cell_print(in);
+		return NULL;
+	}
+
+	tt = uncentering_transformation(in, &new_centering, &new_latt);
+	if ( tt == NULL ) return NULL;
+
+	if ( t != NULL ) *t = tt;
+
+	out = cell_transform(in, tt);
+	if ( out == NULL ) return NULL;
+
+	cell_set_lattice_type(out, new_latt);
+	cell_set_centering(out, new_centering);
+
+	return out;
+}
+
+
+#define MAX_CAND (1024)
+
+static int right_handed_vec(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_in, UnitCell *template_in, int verbose,
+                     const float *tols, 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 angtol = deg2rad(tols[3]);
+	UnitCell *cell;
+	UnitCell *template;
+	UnitCellTransformation *uncentering;
+	UnitCell *new_cell_trans;
+
+	/* "Un-center" the template unit cell to make the comparison easier */
+	template = uncenter_cell(template_in, &uncentering);
+
+	/* The candidate cell is also uncentered, because it might be centered
+	 * if it came from (e.g.) MOSFLM */
+	cell = uncenter_cell(cell_in, NULL);
+
+	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,
+				                       tols[i]) )
+				{
+					continue;
+				}
+
+				if ( ncand[i] == MAX_CAND ) {
+					ERROR("Too many cell candidates - ");
+					ERROR("consider tightening the unit ");
+					ERROR("cell tolerances.\n");
+				} else {
+
+					double fom;
+
+					fom = fabs(lengths[i] - tlen);
+
+					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 = fom;
+
+					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_vec(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]);
+
+	/* Reverse the de-centering transformation */
+	new_cell_trans = cell_transform_inverse(new_cell, uncentering);
+	cell_free(new_cell);
+	cell_set_lattice_type(new_cell, cell_get_lattice_type(template_in));
+	cell_set_centering(new_cell, cell_get_centering(template_in));
+	cell_set_unique_axis(new_cell, cell_get_unique_axis(template_in));
+
+	return new_cell_trans;
+}
+
+
+UnitCell *match_cell_ab(UnitCell *cell_in, UnitCell *template_in)
+{
+	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;
+	float ltl = 5.0;     /* percent */
+	int have_real_a;
+	int have_real_b;
+	int have_real_c;
+	UnitCell *cell;
+	UnitCell *template;
+	UnitCellTransformation *to_given_cell;
+	UnitCell *new_cell;
+	UnitCell *new_cell_trans;
+
+	/* "Un-center" the template unit cell to make the comparison easier */
+	template = uncenter_cell(template_in, &to_given_cell);
+
+	/* The candidate cell is also uncentered, because it might be centered
+	 * if it came from (e.g.) MOSFLM */
+	cell = uncenter_cell(cell_in, NULL);
+
+	/* 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);
+
+	/* 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_vec(real_a, real_b, real_c) ) {
+		real_c.u = -real_c.u;
+		real_c.v = -real_c.v;
+		real_c.w = -real_c.w;
+	}
+
+	new_cell = cell_new_from_direct_axes(real_a, real_b, real_c);
+
+	 /* Reverse the de-centering transformation */
+	new_cell_trans = cell_transform_inverse(new_cell, to_given_cell);
+	cell_free(new_cell);
+	cell_set_lattice_type(new_cell, cell_get_lattice_type(template_in));
+	cell_set_centering(new_cell, cell_get_centering(template_in));
+	cell_set_unique_axis(new_cell, cell_get_unique_axis(template_in));
+
+	return new_cell_trans;
+}
+
+
+/* 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 determine_lattice(UnitCell *cell,
+                              const char *as, const char *bs, const char *cs,
+                              const char *als, const char *bes, const char *gas)
+{
+	int n_right;
+
+	/* Rhombohedral or cubic? */
+	if ( (strcmp(as, bs) == 0) && (strcmp(as, cs) == 0) ) {
+
+		if ( (strcmp(als, "  90.00") == 0)
+		  && (strcmp(bes, "  90.00") == 0)
+		  && (strcmp(gas, "  90.00") == 0) )
+		{
+			/* Cubic.  Unique axis irrelevant. */
+			cell_set_lattice_type(cell, L_CUBIC);
+			return;
+		}
+
+		if ( (strcmp(als, bes) == 0) && (strcmp(als, gas) == 0) ) {
+			/* Rhombohedral.  Unique axis irrelevant. */
+			cell_set_lattice_type(cell, L_RHOMBOHEDRAL);
+			return;
+		}
+
+	}
+
+	if ( (strcmp(als, "  90.00") == 0)
+	  && (strcmp(bes, "  90.00") == 0)
+	  && (strcmp(gas, "  90.00") == 0) )
+	{
+		if ( strcmp(bs, cs) == 0 ) {
+			/* Tetragonal, unique axis a */
+			cell_set_lattice_type(cell, L_TETRAGONAL);
+			cell_set_unique_axis(cell, 'a');
+			return;
+		}
+
+		if ( strcmp(as, cs) == 0 ) {
+			/* Tetragonal, unique axis b */
+			cell_set_lattice_type(cell, L_TETRAGONAL);
+			cell_set_unique_axis(cell, 'b');
+			return;
+		}
+
+		if ( strcmp(as, bs) == 0 ) {
+			/* Tetragonal, unique axis c */
+			cell_set_lattice_type(cell, L_TETRAGONAL);
+			cell_set_unique_axis(cell, 'c');
+			return;
+		}
+
+		/* Orthorhombic.  Unique axis irrelevant, but point group
+		 * can have different orientations. */
+		cell_set_lattice_type(cell, L_ORTHORHOMBIC);
+		return;
+	}
+
+	n_right = 0;
+	if ( strcmp(als, "  90.00") == 0 ) n_right++;
+	if ( strcmp(bes, "  90.00") == 0 ) n_right++;
+	if ( strcmp(gas, "  90.00") == 0 ) n_right++;
+
+	/* Hexgonal or monoclinic? */
+	if ( n_right == 2 ) {
+
+		if ( (strcmp(als, " 120.00") == 0)
+		  && (strcmp(bs, cs) == 0) )
+		{
+			/* Hexagonal, unique axis a */
+			cell_set_lattice_type(cell, L_HEXAGONAL);
+			cell_set_unique_axis(cell, 'a');
+			return;
+		}
+
+		if ( (strcmp(bes, " 120.00") == 0)
+		  && (strcmp(as, cs) == 0) )
+		{
+			/* Hexagonal, unique axis b */
+			cell_set_lattice_type(cell, L_HEXAGONAL);
+			cell_set_unique_axis(cell, 'b');
+			return;
+		}
+
+		if ( (strcmp(gas, " 120.00") == 0)
+		  && (strcmp(as, bs) == 0) )
+		{
+			/* Hexagonal, unique axis c */
+			cell_set_lattice_type(cell, L_HEXAGONAL);
+			cell_set_unique_axis(cell, 'c');
+			return;
+		}
+
+		if ( strcmp(als, "  90.00") != 0 ) {
+			/* Monoclinic, unique axis a */
+			cell_set_lattice_type(cell, L_MONOCLINIC);
+			cell_set_unique_axis(cell, 'a');
+			return;
+		}
+
+		if ( strcmp(bes, "  90.00") != 0 ) {
+			/* Monoclinic, unique axis b */
+			cell_set_lattice_type(cell, L_MONOCLINIC);
+			cell_set_unique_axis(cell, 'b');
+			return;
+		}
+
+		if ( strcmp(gas, "  90.00") != 0 ) {
+			/* Monoclinic, unique axis c */
+			cell_set_lattice_type(cell, L_MONOCLINIC);
+			cell_set_unique_axis(cell, 'c');
+			return;
+		}
+	}
+
+	/* Triclinic, unique axis irrelevant. */
+	cell_set_lattice_type(cell, L_TRICLINIC);
+}
+
+
+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];
+			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));
+
+			determine_lattice(cell, as, bs, cs, als, bes, gas);
+
+			if ( strlen(line) > 65 ) {
+				cell_set_centering(cell, line[55]);
+			} else {
+				cell_set_pointgroup(cell, "1");
+				ERROR("CRYST1 line without centering.\n");
+			}
+
+			break;  /* Done */
+		}
+
+	} while ( rval != NULL );
+
+	fclose(fh);
+
+	validate_cell(cell);
+
+	return cell;
+}
+
+
+/* Force the linker to bring in CBLAS to make GSL happy */
+void cell_fudge_gslcblas()
+{
+        STATUS("%p\n", cblas_sgemm);
+}
+
+
+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;
+}
+
+
+/**
+ * validate_cell:
+ * @cell: A %UnitCell to validate
+ *
+ * Perform some checks for crystallographic validity @cell, such as that the
+ * lattice is a conventional Bravais lattice.
+ * Warnings are printied if any of the checks are failed.
+ *
+ * Returns: true if cell is invalid.
+ *
+ */
+int validate_cell(UnitCell *cell)
+{
+	int err = 0;
+	char cen, ua;
+
+	if ( !cell_is_sensible(cell) ) {
+		ERROR("Warning: Unit cell parameters are not sensible.\n");
+		err = 1;
+	}
+
+	if ( !bravais_lattice(cell) ) {
+		ERROR("Warning: Unit cell is not a conventional Bravais"
+		      " lattice.\n");
+		err = 1;
+	}
+
+	if ( !right_handed(cell) ) {
+		ERROR("Warning: Unit cell is not right handed.\n");
+		err = 1;
+	}
+
+	cen = cell_get_centering(cell);
+	ua = cell_get_unique_axis(cell);
+	if ( (cen == 'A') && (ua != 'a') ) {
+		ERROR("Warning: centering doesn't match unique axis.\n");
+		err = 1;
+	}
+	if ( (cen == 'B') && (ua != 'b') ) {
+		ERROR("Warning: centering doesn't match unique axis.\n");
+		err = 1;
+	}
+	if ( (cen == 'C') && (ua != 'c') ) {
+		ERROR("Warning: centering doesn't match unique axis.\n");
+		err = 1;
+	}
+
+	return err;
+}
+
+
+/**
+ * forbidden_reflection:
+ * @cell: A %UnitCell
+ * @h: h index to check
+ * @k: k index to check
+ * @l: l index to check
+ *
+ * Returns: true if this reflection is forbidden.
+ *
+ */
+int forbidden_reflection(UnitCell *cell,
+                         signed int h, signed int k, signed int l)
+{
+	char cen;
+
+	cen = cell_get_centering(cell);
+
+	/* Reflection conditions here must match the transformation matrices
+	 * in uncentering_transformation().  tests/centering_check verifies
+	 * this (amongst other things). */
+
+	if ( cen == 'P' ) return 0;
+	if ( cen == 'R' ) return 0;
+
+	if ( cen == 'A' ) return (k+l) % 2;
+	if ( cen == 'B' ) return (h+l) % 2;
+	if ( cen == 'C' ) return (h+k) % 2;
+
+	if ( cen == 'I' ) return (h+k+l) % 2;
+	if ( cen == 'F' ) return ((h+k) % 2) || ((h+l) % 2) || ((k+l) % 2);
+
+	/* Obverse setting */
+	if ( cen == 'H' ) return (-h+k+l) % 3;
+
+	return 0;
+}
diff --git a/libcrystfel/src/cell-utils.h b/libcrystfel/src/cell-utils.h
new file mode 100644
index 00000000..f92ab22d
--- /dev/null
+++ b/libcrystfel/src/cell-utils.h
@@ -0,0 +1,69 @@
+/*
+ * cell-utils.h
+ *
+ * Unit Cell utility functions
+ *
+ * Copyright © 2012 Deutsches Elektronen-Synchrotron DESY,
+ *                  a research centre of the Helmholtz Association.
+ * Copyright © 2012 Lorenzo Galli
+ *
+ * Authors:
+ *   2009-2012 Thomas White <taw@physics.org>
+ *   2012      Lorenzo Galli
+ *
+ * This file is part of CrystFEL.
+ *
+ * CrystFEL is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU General Public License as published by
+ * the Free Software Foundation, either version 3 of the License, or
+ * (at your option) any later version.
+ *
+ * CrystFEL is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU General Public License
+ * along with CrystFEL.  If not, see <http://www.gnu.org/licenses/>.
+ *
+ */
+
+#ifndef CELL_UTILS_H
+#define CELL_UTILS_H
+
+#ifdef HAVE_CONFIG_H
+#include <config.h>
+#endif
+
+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 *tempcell, int verbose,
+                            const float *ltl, int reduce);
+
+extern UnitCell *match_cell_ab(UnitCell *cell, UnitCell *tempcell);
+
+extern UnitCell *load_cell_from_pdb(const char *filename);
+
+extern int cell_is_sensible(UnitCell *cell);
+
+extern int validate_cell(UnitCell *cell);
+
+extern UnitCell *uncenter_cell(UnitCell *in, UnitCellTransformation **tr);
+
+extern int bravais_lattice(UnitCell *cell);
+
+extern int right_handed(UnitCell *cell);
+
+extern const char *str_lattice(LatticeType l);
+
+extern int forbidden_reflection(UnitCell *cell,
+                                signed int h, signed int k, signed int l);
+
+#endif	/* CELL_UTILS_H */
diff --git a/libcrystfel/src/cell.c b/libcrystfel/src/cell.c
index 3d2a28bf..e7c9dace 100644
--- a/libcrystfel/src/cell.c
+++ b/libcrystfel/src/cell.c
@@ -1,7 +1,7 @@
 /*
  * cell.c
  *
- * Unit Cell Calculations
+ * A class representing a unit cell
  *
  * Copyright © 2012 Deutsches Elektronen-Synchrotron DESY,
  *                  a research centre of the Helmholtz Association.
@@ -60,10 +60,6 @@
  */
 
 
-
-/* Weighting factor of lengths relative to angles */
-#define LWEIGHT (10.0e-9)
-
 typedef enum {
 	CELL_REP_CRYST,
 	CELL_REP_CART,
@@ -92,8 +88,10 @@ struct _unitcell {
 	double ays;	double bys;	double cys;
 	double azs;	double bzs;	double czs;
 
-	char *pointgroup;
-	char *spacegroup;
+	char        *pointgroup;
+	LatticeType  lattice_type;
+	char         centering;
+	char         unique_axis;
 };
 
 
@@ -125,7 +123,9 @@ UnitCell *cell_new()
 	cell->rep = CELL_REP_CRYST;
 
 	cell->pointgroup = strdup("1");
-	cell->spacegroup = strdup("P 1");
+	cell->lattice_type = L_TRICLINIC;
+	cell->centering = 'P';
+	cell->unique_axis = 'c';
 
 	return cell;
 }
@@ -142,7 +142,6 @@ void cell_free(UnitCell *cell)
 {
 	if ( cell == NULL ) return;
 	free(cell->pointgroup);
-	free(cell->spacegroup);
 	free(cell);
 }
 
@@ -261,7 +260,9 @@ UnitCell *cell_new_from_cell(UnitCell *orig)
 	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);
+	cell_set_lattice_type(new, orig->lattice_type);
+	cell_set_centering(new, orig->centering);
+	cell_set_unique_axis(new, orig->unique_axis);
 
 	return new;
 }
@@ -282,17 +283,28 @@ void cell_set_reciprocal(UnitCell *cell,
 }
 
 
-void cell_set_spacegroup(UnitCell *cell, const char *sym)
+void cell_set_pointgroup(UnitCell *cell, const char *sym)
 {
-	free(cell->spacegroup);
-	cell->spacegroup = strdup(sym);
+	free(cell->pointgroup);
+	cell->pointgroup = strdup(sym);
 }
 
 
-void cell_set_pointgroup(UnitCell *cell, const char *sym)
+void cell_set_centering(UnitCell *cell, char centering)
 {
-	free(cell->pointgroup);
-	cell->pointgroup = strdup(sym);
+	cell->centering = centering;
+}
+
+
+void cell_set_lattice_type(UnitCell *cell, LatticeType lattice_type)
+{
+	cell->lattice_type = lattice_type;
+}
+
+
+void cell_set_unique_axis(UnitCell *cell, char unique_axis)
+{
+	cell->unique_axis = unique_axis;
 }
 
 
@@ -555,18 +567,25 @@ const char *cell_get_pointgroup(UnitCell *cell)
 }
 
 
-const char *cell_get_spacegroup(UnitCell *cell)
+char cell_get_centering(UnitCell *cell)
 {
-	return cell->spacegroup;
+	return cell->centering;
 }
 
 
+LatticeType cell_get_lattice_type(UnitCell *cell)
+{
+	return cell->lattice_type;
+}
 
 
+char cell_get_unique_axis(UnitCell *cell)
+{
+	return cell->unique_axis;
+}
 
-/********************************* Utilities **********************************/
 
-static const char *cell_rep(UnitCell *cell)
+const char *cell_rep(UnitCell *cell)
 {
 	switch ( cell->rep ) {
 
@@ -585,614 +604,307 @@ static const char *cell_rep(UnitCell *cell)
 }
 
 
+struct _unitcelltransformation
+{
+	gsl_matrix *m;
+};
+
+
 /**
- * cell_rotate:
- * @in: A %UnitCell to rotate
- * @quat: A %quaternion
+ * tfn_inverse:
+ * @t: A %UnitCellTransformation.
  *
- * Rotate a %UnitCell using a %quaternion.
+ * Calculates the inverse of @t.  That is, if you apply cell_transform() to a
+ * %UnitCell using @t, and then apply cell_transform() to the result using
+ * tfn_inverse(@t) instead of @t, you will recover the same lattice vectors
+ * (but note that the lattice type, centering and unique axis information will
+ * be lost).
  *
- * Returns: a newly allocated rotated copy of @in.
+ * Returns: The inverse of @t.
  *
  */
-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)
+UnitCellTransformation *tfn_inverse(UnitCellTransformation *t)
 {
-	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;
-}
+	int s;
+	gsl_matrix *m;
+	gsl_matrix *inv;
+	gsl_permutation *perm;
+	UnitCellTransformation *out;
 
+	m = gsl_matrix_alloc(3, 3);
+	if ( m == NULL ) return NULL;
 
-/* Attempt to make 'cell' fit into 'template' somehow */
-UnitCell *match_cell(UnitCell *cell, UnitCell *template, int verbose,
-                     const float *tols, 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 angtol = deg2rad(tols[3]);
-
-	if ( cell_get_reciprocal(template, &asx, &asy, &asz,
-	                         &bsx, &bsy, &bsz,
-	                         &csx, &csy, &csz) ) {
-		ERROR("Couldn't get reciprocal cell for template.\n");
+	out = tfn_identity();
+	if ( out == NULL ) {
+		gsl_matrix_free(m);
 		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);
+	gsl_matrix_memcpy(m, t->m);
 
-	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");
+	perm = gsl_permutation_alloc(m->size1);
+	if ( perm == NULL ) {
+		ERROR("Couldn't allocate permutation\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,
-				                       tols[i]) )
-				{
-					continue;
-				}
-
-				if ( ncand[i] == MAX_CAND ) {
-					ERROR("Too many cell candidates - ");
-					ERROR("consider tightening the unit ");
-					ERROR("cell tolerances.\n");
-				} else {
-
-					double fom;
-
-					fom = fabs(lengths[i] - tlen);
-
-					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 = fom;
-
-					ncand[i]++;
-
-				}
-
-			}
-
-		}
-		}
-		}
-
-	}
-	}
-	}
-
-	if ( verbose ) {
-		STATUS("Candidates: %i %i %i\n", ncand[0], ncand[1], ncand[2]);
+	inv = gsl_matrix_alloc(m->size1, m->size2);
+	if ( inv == NULL ) {
+		ERROR("Couldn't allocate inverse\n");
+		gsl_permutation_free(perm);
+		return NULL;
 	}
-
-	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;
-			}
-
-		}
-
+	if ( gsl_linalg_LU_decomp(m, perm, &s) ) {
+		ERROR("Couldn't decompose matrix\n");
+		gsl_permutation_free(perm);
+		return NULL;
 	}
+	if ( gsl_linalg_LU_invert(m, perm, inv)  ) {
+		ERROR("Couldn't invert matrix\n");
+		gsl_permutation_free(perm);
+		return NULL;
 	}
+	gsl_permutation_free(perm);
 
-	free(cand[0]);
-	free(cand[1]);
-	free(cand[2]);
-
-	return new_cell;
+	gsl_matrix_free(out->m);
+	gsl_matrix_free(m);
+	out->m = inv;
+	return out;
 }
 
 
-
-UnitCell *match_cell_ab(UnitCell *cell, UnitCell *template)
+/**
+ * cell_transform:
+ * @cell: A %UnitCell.
+ * @t: A %UnitCellTransformation.
+ *
+ * Applies @t to @cell.  Note that the lattice type, centering and unique axis
+ * information will not be preserved.
+ *
+ * Returns: Transformed copy of @cell.
+ *
+ */
+UnitCell *cell_transform(UnitCell *cell, UnitCellTransformation *t)
 {
+	UnitCell *out;
 	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;
-	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);
-
-	/* 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;
-		}
-	}
+	gsl_matrix *m;
+	gsl_matrix *a;
 
-	/* Have we matched both a and b? */
-	if ( !(have_real_a && have_real_b) ) return NULL;
+	if ( t == NULL ) 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;
-		}
-	}
+	out = cell_new_from_cell(cell);
+	if ( out == NULL ) return NULL;
 
-	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;
-	}
+	cell_get_cartesian(out, &ax, &ay, &az,
+	                        &bx, &by, &bz,
+	                        &cx, &cy, &cz);
 
-	/* 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;
+	m = gsl_matrix_alloc(3,3);
+	a = gsl_matrix_calloc(3,3);
+	if ( (m == NULL) || (a == NULL) ) {
+		cell_free(out);
+		return NULL;
 	}
 
-	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);
+	gsl_matrix_set(m, 0, 0, ax);
+	gsl_matrix_set(m, 0, 1, ay);
+	gsl_matrix_set(m, 0, 2, az);
+	gsl_matrix_set(m, 1, 0, bx);
+	gsl_matrix_set(m, 1, 1, by);
+	gsl_matrix_set(m, 1, 2, bz);
+	gsl_matrix_set(m, 2, 0, cx);
+	gsl_matrix_set(m, 2, 1, cy);
+	gsl_matrix_set(m, 2, 2, cz);
 
-	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) );
+	gsl_blas_dgemm(CblasNoTrans, CblasNoTrans, 1.0, t->m, m, 1.0, a);
 
-	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));
+	cell_set_cartesian(out, gsl_matrix_get(a, 0, 0),
+	                        gsl_matrix_get(a, 0, 1),
+	                        gsl_matrix_get(a, 0, 2),
+	                        gsl_matrix_get(a, 1, 0),
+	                        gsl_matrix_get(a, 1, 1),
+	                        gsl_matrix_get(a, 1, 2),
+	                        gsl_matrix_get(a, 2, 0),
+	                        gsl_matrix_get(a, 2, 1),
+	                        gsl_matrix_get(a, 2, 2));
 
-	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);
+	gsl_matrix_free(a);
+	gsl_matrix_free(m);
 
-	return oneoverd / 2;
+	return out;
 }
 
 
-static void cell_set_pointgroup_from_pdb(UnitCell *cell, const char *sym)
+/**
+ * cell_transform_inverse:
+ * @cell: A %UnitCell.
+ * @t: A %UnitCellTransformation.
+ *
+ * Applies the inverse of @t to @cell.
+ *
+ * Returns: Transformed copy of @cell.
+ *
+ */
+UnitCell *cell_transform_inverse(UnitCell *cell, UnitCellTransformation *t)
 {
-	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);
-	}
+	UnitCellTransformation *inv;
+	UnitCell *out;
+
+	inv = tfn_inverse(t);
+	out = cell_transform(cell, inv);
+	tfn_free(inv);
+	return out;
 }
 
 
-UnitCell *load_cell_from_pdb(const char *filename)
+/**
+ * tfn_identity:
+ *
+ * Returns: A %UnitCellTransformation corresponding to an identity operation.
+ *
+ */
+UnitCellTransformation *tfn_identity()
 {
-	FILE *fh;
-	char *rval;
-	UnitCell *cell = NULL;
+	UnitCellTransformation *tfn;
+
+	tfn = calloc(1, sizeof(UnitCellTransformation));
+	if ( tfn == NULL ) return NULL;
 
-	fh = fopen(filename, "r");
-	if ( fh == NULL ) {
-		ERROR("Couldn't open '%s'\n", filename);
+	tfn->m = gsl_matrix_alloc(3, 3);
+	if ( tfn->m == NULL ) {
+		free(tfn);
 		return NULL;
 	}
 
-	do {
+	gsl_matrix_set_identity(tfn->m);
 
-		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;
-			}
+	return tfn;
+}
 
-			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");
-			}
+/**
+ * tfn_combine:
+ * @t: A %UnitCellTransformation
+ * @na: Pointer to three doubles representing naa, nab, nac
+ * @nb: Pointer to three doubles representing nba, nbb, nbc
+ * @nc: Pointer to three doubles representing nca, ncb, ncc
+ *
+ * Updates @t such that it represents its previous transformation followed by
+ * a new transformation, corresponding to letting a = naa*a + nab*b + nac*c.
+ * Likewise, a = nba*a + nbb*b + nbc*c and c = nca*a + ncb*b + ncc*c.
+ *
+ */
+void tfn_combine(UnitCellTransformation *t, double *na, double *nb, double *nc)
+{
+	gsl_matrix *a;
+	gsl_matrix *n;
 
-			break;  /* Done */
-		}
+	n = gsl_matrix_alloc(3, 3);
+	a = gsl_matrix_calloc(3, 3);
+	if ( (n == NULL) || (a == NULL) ) {
+		return;
+	}
+	gsl_matrix_set(n, 0, 0, na[0]);
+	gsl_matrix_set(n, 0, 1, na[1]);
+	gsl_matrix_set(n, 0, 2, na[2]);
+	gsl_matrix_set(n, 1, 0, nb[0]);
+	gsl_matrix_set(n, 1, 1, nb[1]);
+	gsl_matrix_set(n, 1, 2, nb[2]);
+	gsl_matrix_set(n, 2, 0, nc[0]);
+	gsl_matrix_set(n, 2, 1, nc[1]);
+	gsl_matrix_set(n, 2, 2, nc[2]);
 
-	} while ( rval != NULL );
+	free(na);
+	free(nb);
+	free(nc);
 
-	fclose(fh);
+	gsl_blas_dgemm(CblasNoTrans, CblasNoTrans, 1.0, n, t->m, 1.0, a);
 
-	return cell;
+	gsl_matrix_free(t->m);
+	t->m = a;
+	gsl_matrix_free(n);
 }
 
 
-/* Force the linker to bring in CBLAS to make GSL happy */
-void cell_fudge_gslcblas()
+/**
+ * tfn_vector:
+ * @a: Amount of "a" to include in new vector
+ * @b: Amount of "b" to include in new vector
+ * @c: Amount of "c" to include in new vector
+ *
+ * This is a convenience function to use when sending vectors to tfn_combine():
+ *	tfn_combine(tfn, tfn_vector(1,0,0),
+ *	                 tfn_vector(0,2,0),
+ *	                 tfn_vector(0,0,1));
+ *
+ */
+double *tfn_vector(double a, double b, double c)
 {
-        STATUS("%p\n", cblas_sgemm);
+	double *vec = malloc(3*sizeof(double));
+	if ( vec == NULL ) return NULL;
+	vec[0] = a;  vec[1] = b;  vec[2] = c;
+	return vec;
 }
 
 
-UnitCell *rotate_cell(UnitCell *in, double omega, double phi, double rot)
+/**
+ * tfn_print:
+ * @t: A %UnitCellTransformation
+ *
+ * Prints information about @t to stderr.
+ *
+ */
+void tfn_print(UnitCellTransformation *t)
 {
-	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);
+	gsl_permutation *perm;
+	gsl_matrix *lu;
+	int s;
 
-	return out;
+	STATUS("New a = %+.2fa %+.2fb %+.2fc\n", gsl_matrix_get(t->m, 0, 0),
+	                                         gsl_matrix_get(t->m, 0, 1),
+	                                         gsl_matrix_get(t->m, 0, 2));
+	STATUS("New b = %+.2fa %+.2fb %+.2fc\n", gsl_matrix_get(t->m, 1, 0),
+	                                         gsl_matrix_get(t->m, 1, 1),
+	                                         gsl_matrix_get(t->m, 1, 2));
+	STATUS("New c = %+.2fa %+.2fb %+.2fc\n", gsl_matrix_get(t->m, 2, 0),
+	                                         gsl_matrix_get(t->m, 2, 1),
+	                                         gsl_matrix_get(t->m, 2, 2));
+	lu = gsl_matrix_alloc(3, 3);
+	if ( lu == NULL ) {
+		ERROR("Couldn't allocate LU decomposition.\n");
+		return;
+	}
+
+	gsl_matrix_memcpy(lu, t->m);
+
+	perm = gsl_permutation_alloc(t->m->size1);
+	if ( perm == NULL ) {
+		ERROR("Couldn't allocate permutation.\n");
+		gsl_matrix_free(lu);
+		return;
+	}
+	if ( gsl_linalg_LU_decomp(lu, perm, &s) ) {
+		ERROR("LU decomposition failed.\n");
+		gsl_permutation_free(perm);
+		gsl_matrix_free(lu);
+		return;
+	}
+
+	STATUS("Transformation determinant = %.2f\n", gsl_linalg_LU_det(lu, s));
 }
 
 
-int cell_is_sensible(UnitCell *cell)
+/**
+ * tfn_free:
+ * @t: A %UnitCellTransformation
+ *
+ * Frees all resources associated with @t.
+ *
+ */
+void tfn_free(UnitCellTransformation *t)
 {
-	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;
+	gsl_matrix_free(t->m);
+	free(t);
 }
diff --git a/libcrystfel/src/cell.h b/libcrystfel/src/cell.h
index bd2719dd..c7b8f8d6 100644
--- a/libcrystfel/src/cell.h
+++ b/libcrystfel/src/cell.h
@@ -1,7 +1,7 @@
 /*
  * cell.h
  *
- * Unit Cell Calculations
+ * A class representing a unit cell
  *
  * Copyright © 2012 Deutsches Elektronen-Synchrotron DESY,
  *                  a research centre of the Helmholtz Association.
@@ -48,6 +48,17 @@ struct rvec
 	double   w;
 };
 
+typedef enum
+{
+	L_TRICLINIC,
+	L_MONOCLINIC,
+	L_ORTHORHOMBIC,
+	L_TETRAGONAL,
+	L_RHOMBOHEDRAL,
+	L_HEXAGONAL,
+	L_CUBIC
+} LatticeType;
+
 
 /**
  * UnitCell:
@@ -57,6 +68,15 @@ struct rvec
  **/
 typedef struct _unitcell UnitCell;
 
+
+/**
+ * UnitCellTransformation:
+ *
+ * This opaque data structure represents a tranformation of a unit cell, such
+ * as a rotation or a centering operation.
+ **/
+typedef struct _unitcelltransformation UnitCellTransformation;
+
 extern UnitCell *cell_new(void);
 extern UnitCell *cell_new_from_cell(UnitCell *orig);
 extern void cell_free(UnitCell *cell);
@@ -82,7 +102,6 @@ extern void cell_set_parameters(UnitCell *cell, double a, double b, double c,
 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);
 
 
@@ -106,24 +125,27 @@ extern void cell_set_reciprocal(UnitCell *cell,
 
 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 LatticeType cell_get_lattice_type(UnitCell *cell);
+extern void cell_set_lattice_type(UnitCell *cell, LatticeType lattice_type);
 
-extern void cell_print(UnitCell *cell);
+extern char cell_get_centering(UnitCell *cell);
+extern void cell_set_centering(UnitCell *cell, char centering);
 
-extern UnitCell *match_cell(UnitCell *cell, UnitCell *tempcell, int verbose,
-                            const float *ltl, int reduce);
+extern char cell_get_unique_axis(UnitCell *cell);
+extern void cell_set_unique_axis(UnitCell *cell, char unique_axis);
 
-extern UnitCell *match_cell_ab(UnitCell *cell, UnitCell *tempcell);
+extern const char *cell_rep(UnitCell *cell);
 
-extern UnitCell *load_cell_from_pdb(const char *filename);
+extern UnitCell *cell_transform(UnitCell *cell, UnitCellTransformation *t);
+extern UnitCell *cell_transform_inverse(UnitCell *cell,
+                                        UnitCellTransformation *t);
 
-extern int cell_is_sensible(UnitCell *cell);
+extern UnitCellTransformation *tfn_identity(void);
+extern void tfn_combine(UnitCellTransformation *t,
+                        double *na, double *nb, double *nc);
+extern void tfn_print(UnitCellTransformation *t);
+extern UnitCellTransformation *tfn_inverse(UnitCellTransformation *t);
+extern double *tfn_vector(double a, double b, double c);
+extern void tfn_free(UnitCellTransformation *t);
 
 #endif	/* CELL_H */
diff --git a/libcrystfel/src/geometry.c b/libcrystfel/src/geometry.c
index 24669aef..091b4fed 100644
--- a/libcrystfel/src/geometry.c
+++ b/libcrystfel/src/geometry.c
@@ -36,6 +36,7 @@
 
 #include "utils.h"
 #include "cell.h"
+#include "cell-utils.h"
 #include "image.h"
 #include "peaks.h"
 #include "beam-parameters.h"
diff --git a/libcrystfel/src/index.c b/libcrystfel/src/index.c
index 102a5312..7d912902 100644
--- a/libcrystfel/src/index.c
+++ b/libcrystfel/src/index.c
@@ -49,6 +49,7 @@
 #include "index-priv.h"
 #include "reax.h"
 #include "geometry.h"
+#include "cell-utils.h"
 
 
 /* Base class constructor for unspecialised indexing private data */
diff --git a/libcrystfel/src/mosflm.c b/libcrystfel/src/mosflm.c
index ed118aa4..63166919 100644
--- a/libcrystfel/src/mosflm.c
+++ b/libcrystfel/src/mosflm.c
@@ -272,14 +272,74 @@ static void mosflm_sendline(const char *line, struct mosflm_data *mosflm)
 }
 
 
+/* Turn what we know about the unit cell into something which we can give to
+ * MOSFLM to make it give us only indexing results compatible with the cell. */
+static const char *spacegroup_for_lattice(UnitCell *cell)
+{
+	LatticeType latt;
+	char centering;
+	char ua;
+	char *g = NULL;
+	char *result;
+
+	latt = cell_get_lattice_type(cell);
+	centering = cell_get_centering(cell);
+	ua = cell_get_unique_axis(cell);
+
+	switch ( latt )
+	{
+		case L_TRICLINIC :
+		g = "1";
+		break;
+
+		case L_MONOCLINIC :
+		/* "2 1 1", "1 2 1" or "1 1 2" depending on unique axis */
+		if ( ua == 'a' ) g = "2 1 1";
+		if ( ua == 'b' ) g = "1 2 1";
+		if ( ua == 'c' ) g = "1 1 2";
+		break;
+
+		case L_ORTHORHOMBIC :
+		g = "2 2 2";
+		break;
+
+		case L_TETRAGONAL :
+		/* "4 1 1", "1 4 1" or "1 1 4" depending on unique axis */
+		if ( ua == 'a' ) g = "4 1 1";
+		if ( ua == 'b' ) g = "1 4 1";
+		if ( ua == 'c' ) g = "1 1 4";
+		break;
+
+		case L_RHOMBOHEDRAL :
+		g = "3";
+		break;
+
+		case L_HEXAGONAL :
+		/* "6 1 1", "1 6 1" or "1 1 6" depending on unique axis */
+		if ( ua == 'a' ) g = "6 1 1";
+		if ( ua == 'b' ) g = "1 6 1";
+		if ( ua == 'c' ) g = "6";
+		break;
+
+		case L_CUBIC :
+		g = "2 3";
+		break;
+	}
+	assert(g != NULL);
+
+	result = malloc(32);
+	if ( result == NULL ) return NULL;
+
+	snprintf(result, 31, "%c%s", centering, g);
+
+	return result;
+}
+
+
 static void mosflm_send_next(struct image *image, struct mosflm_data *mosflm)
 {
 	char tmp[256];
-	char symm[32];
-	const char *sg;
 	double wavelength;
-	double a, b, c, alpha, beta, gamma;
-	int i, j;
 
 	switch ( mosflm->step ) {
 
@@ -291,29 +351,8 @@ static void mosflm_send_next(struct image *image, struct mosflm_data *mosflm)
 
 		case 2 :
 		if ( mosflm->target_cell != NULL ) {
-			cell_get_parameters(mosflm->target_cell, &a, &b, &c,
-				            &alpha, &beta, &gamma);
-			snprintf(tmp, 255,
-			         "CELL %6.2f %6.2f %6.2f %6.2f %6.2f %6.2f\n",
-		                 a*1e10, b*1e10, c*1e10,
-		                 rad2deg(alpha), rad2deg(beta), rad2deg(gamma));
-			mosflm_sendline(tmp, mosflm);
-		} else {
-			mosflm_sendline("# Do nothing\n", mosflm);
-		}
-		break;
-
-		case 3 :
-		if ( mosflm->target_cell != NULL ) {
-			sg = cell_get_spacegroup(mosflm->target_cell);
-			/* Remove white space from space group */
-			j = 0;
-			for ( i=0; i<strlen(sg); i++ ) {
-				if (sg[i] != ' ') {
-					symm[j++] = sg[i];
-				}
-			}
-			symm[j] = '\0';
+			const char *symm;
+			symm = spacegroup_for_lattice(mosflm->target_cell);
 			snprintf(tmp, 255, "SYMM %s\n", symm);
 			mosflm_sendline(tmp, mosflm);
 		} else {
@@ -321,31 +360,31 @@ static void mosflm_send_next(struct image *image, struct mosflm_data *mosflm)
 		}
 		break;
 
-		case 4 :
+		case 3 :
 		mosflm_sendline("DISTANCE 67.8\n", mosflm);
 		break;
 
-		case 5 :
+		case 4 :
 		mosflm_sendline("BEAM 0.0 0.0\n", mosflm);
 		break;
 
-		case 6 :
+		case 5 :
 		wavelength = image->lambda*1e10;
 		snprintf(tmp, 255, "WAVELENGTH %10.5f\n", wavelength);
 		mosflm_sendline(tmp, mosflm);
 		break;
 
-		case 7 :
+		case 6 :
 		snprintf(tmp, 255, "NEWMAT %s\n", mosflm->newmatfile);
 		mosflm_sendline(tmp, mosflm);
 		break;
 
-		case 8 :
+		case 7 :
 		snprintf(tmp, 255, "IMAGE %s phi 0 0\n", mosflm->imagefile);
 		mosflm_sendline(tmp, mosflm);
 		break;
 
-		case 9 :
+		case 8 :
 		snprintf(tmp, 255, "AUTOINDEX DPS FILE %s"
 		                   " IMAGE 1 MAXCELL 1000 REFINE\n",
 		         mosflm->sptfile);
@@ -357,7 +396,7 @@ static void mosflm_send_next(struct image *image, struct mosflm_data *mosflm)
 		mosflm_sendline(tmp, mosflm);
 		break;
 
-		case 10 :
+		case 9 :
 		mosflm_sendline("GO\n", mosflm);
 		mosflm->finished_ok = 1;
 		break;
diff --git a/libcrystfel/src/peaks.c b/libcrystfel/src/peaks.c
index b87eb56e..f7f6c650 100644
--- a/libcrystfel/src/peaks.c
+++ b/libcrystfel/src/peaks.c
@@ -51,6 +51,7 @@
 #include "filters.h"
 #include "reflist-utils.h"
 #include "beam-parameters.h"
+#include "cell-utils.h"
 
 
 /* Degree of polarisation of X-ray beam */
diff --git a/libcrystfel/src/reax.c b/libcrystfel/src/reax.c
index 5cfa908a..7529d12f 100644
--- a/libcrystfel/src/reax.c
+++ b/libcrystfel/src/reax.c
@@ -41,12 +41,12 @@
 #include <gsl/gsl_vector.h>
 #include <gsl/gsl_linalg.h>
 #include <gsl/gsl_eigen.h>
-#include <gsl/gsl_blas.h>
 
 #include "image.h"
 #include "utils.h"
 #include "peaks.h"
 #include "cell.h"
+#include "cell-utils.h"
 #include "index.h"
 #include "index-priv.h"
 
@@ -744,7 +744,7 @@ static double max_feature_resolution(ImageFeatureList *flist)
 }
 
 
-static int right_handed(struct rvec a, struct rvec b, struct rvec c)
+static int right_handed_vec(struct rvec a, struct rvec b, struct rvec c)
 {
 	struct rvec aCb;
 	double aCb_dot_c;
@@ -957,7 +957,7 @@ static void assemble_cells_from_candidates(struct image *image,
 		bi.u = vj.x;  bi.v = vj.y;  bi.w = vj.z;
 		ci.u = vk.x;  ci.v = vk.y;  ci.w = vk.z;
 
-		if ( !right_handed(ai, bi, ci) ) continue;
+		if ( !right_handed_vec(ai, bi, ci) ) continue;
 
 		/* We have three vectors with the right angles */
 		cnew = cell_new_from_direct_axes(ai, bi, ci);
diff --git a/libcrystfel/src/reflist-utils.c b/libcrystfel/src/reflist-utils.c
index 2e6715f2..f2292929 100644
--- a/libcrystfel/src/reflist-utils.c
+++ b/libcrystfel/src/reflist-utils.c
@@ -35,6 +35,7 @@
 
 #include "reflist.h"
 #include "cell.h"
+#include "cell-utils.h"
 #include "utils.h"
 #include "reflist-utils.h"
 #include "symmetry.h"
diff --git a/libcrystfel/src/symmetry.c b/libcrystfel/src/symmetry.c
index 0352e441..52a9aae6 100644
--- a/libcrystfel/src/symmetry.c
+++ b/libcrystfel/src/symmetry.c
@@ -68,7 +68,6 @@ struct _symoplist
 	int n_ops;
 	int max_ops;
 	char *name;
-	int *divisors;
 	int num_equivs;
 };
 
@@ -84,7 +83,6 @@ struct _symopmask
 static void alloc_ops(SymOpList *ops)
 {
 	ops->ops = realloc(ops->ops, ops->max_ops*sizeof(struct sym_op));
-	ops->divisors = realloc(ops->divisors, ops->max_ops*sizeof(int));
 }
 
 
@@ -127,7 +125,6 @@ static SymOpList *new_symoplist()
 	new->max_ops = 16;
 	new->n_ops = 0;
 	new->ops = NULL;
-	new->divisors = NULL;
 	new->name = NULL;
 	new->num_equivs = 1;
 	alloc_ops(new);
@@ -259,6 +256,7 @@ int num_equivs(const SymOpList *ops, const SymOpMask *m)
 static signed int *v(signed int h, signed int k, signed int i, signed int l)
 {
 	signed int *vec = malloc(3*sizeof(signed int));
+	if ( vec == NULL ) return NULL;
 	/* Convert back to 3-index form now */
 	vec[0] = h-i;  vec[1] = k-i;  vec[2] = l;
 	return vec;
diff --git a/src/check_hkl.c b/src/check_hkl.c
index 1dfe91a0..84cfdccf 100644
--- a/src/check_hkl.c
+++ b/src/check_hkl.c
@@ -43,6 +43,7 @@
 #include "symmetry.h"
 #include "reflist.h"
 #include "reflist-utils.h"
+#include "cell-utils.h"
 
 
 /* Number of bins for plot of resolution shells */
diff --git a/src/compare_hkl.c b/src/compare_hkl.c
index b6a525fc..2b0e5f90 100644
--- a/src/compare_hkl.c
+++ b/src/compare_hkl.c
@@ -44,6 +44,7 @@
 #include "statistics.h"
 #include "symmetry.h"
 #include "reflist-utils.h"
+#include "cell-utils.h"
 
 
 /* Number of bins for plot of resolution shells */
diff --git a/src/hrs-scaling.c b/src/hrs-scaling.c
index 3ae0d563..399d5f90 100644
--- a/src/hrs-scaling.c
+++ b/src/hrs-scaling.c
@@ -37,7 +37,6 @@
 #include <gsl/gsl_vector.h>
 #include <gsl/gsl_linalg.h>
 #include <gsl/gsl_eigen.h>
-#include <gsl/gsl_blas.h>
 
 #include "image.h"
 #include "peaks.h"
diff --git a/src/indexamajig.c b/src/indexamajig.c
index 074be20a..4b3dcb5d 100644
--- a/src/indexamajig.c
+++ b/src/indexamajig.c
@@ -62,6 +62,7 @@
 #include "geometry.h"
 #include "stream.h"
 #include "reflist-utils.h"
+#include "cell-utils.h"
 
 #include "im-sandbox.h"
 
diff --git a/src/partial_sim.c b/src/partial_sim.c
index 62c72ceb..e9a7003e 100644
--- a/src/partial_sim.c
+++ b/src/partial_sim.c
@@ -40,14 +40,15 @@
 #include <assert.h>
 #include <pthread.h>
 
-#include <utils.h>
-#include <reflist-utils.h>
-#include <symmetry.h>
-#include <beam-parameters.h>
-#include <detector.h>
-#include <geometry.h>
-#include <stream.h>
-#include <thread-pool.h>
+#include "image.h"
+#include "utils.h"
+#include "reflist-utils.h"
+#include "symmetry.h"
+#include "beam-parameters.h"
+#include "geometry.h"
+#include "stream.h"
+#include "thread-pool.h"
+#include "cell-utils.h"
 
 /* Number of bins for partiality graph */
 #define NBINS 50
diff --git a/src/pattern_sim.c b/src/pattern_sim.c
index f9831e79..06a1f5e5 100644
--- a/src/pattern_sim.c
+++ b/src/pattern_sim.c
@@ -42,6 +42,7 @@
 #include "diffraction.h"
 #include "diffraction-gpu.h"
 #include "cell.h"
+#include "cell-utils.h"
 #include "utils.h"
 #include "hdf5-file.h"
 #include "detector.h"
diff --git a/src/post-refinement.c b/src/post-refinement.c
index 78c47d07..7410d931 100644
--- a/src/post-refinement.c
+++ b/src/post-refinement.c
@@ -37,7 +37,6 @@
 #include <gsl/gsl_vector.h>
 #include <gsl/gsl_linalg.h>
 #include <gsl/gsl_eigen.h>
-#include <gsl/gsl_blas.h>
 
 #include "image.h"
 #include "post-refinement.h"
@@ -45,6 +44,7 @@
 #include "symmetry.h"
 #include "geometry.h"
 #include "cell.h"
+#include "cell-utils.h"
 
 
 /* Maximum number of iterations of NLSq to do for each image per macrocycle. */
diff --git a/src/powder_plot.c b/src/powder_plot.c
index dd7f3f29..7917ccd0 100644
--- a/src/powder_plot.c
+++ b/src/powder_plot.c
@@ -51,6 +51,7 @@
 #include "beam-parameters.h"
 #include "reflist-utils.h"
 #include "symmetry.h"
+#include "cell-utils.h"
 
 
 struct bin_stats {
diff --git a/src/render_hkl.c b/src/render_hkl.c
index 9e7a6dda..7818b304 100644
--- a/src/render_hkl.c
+++ b/src/render_hkl.c
@@ -48,6 +48,7 @@
 #include "render_hkl.h"
 #include "reflist.h"
 #include "reflist-utils.h"
+#include "cell-utils.h"
 
 
 #define KEY_FILENAME "key.pdf"
@@ -157,7 +158,7 @@ static void draw_circles(double xh, double xk, double xl,
                          cairo_t *dctx, int wght, double boost, int colscale,
                          UnitCell *cell, double radius, double theta,
                          double as, double bs, double cx, double cy,
-                         double scale, double max_val)
+                         double scale, double max_val, signed int zone)
 {
 	Reflection *refl;
 	RefListIterator *iter;
@@ -192,7 +193,7 @@ static void draw_circles(double xh, double xk, double xl,
 			get_equiv(sym, m, i, ha, ka, la, &h, &k, &l);
 
 			/* Is the reflection in the zone? */
-			if ( h*zh + k*zk + l*zl != 0 ) continue;
+			if ( h*zh + k*zk + l*zl != zone ) continue;
 
 			xi = (h*xh + k*xk + l*xl) / nx;
 			yi = (h*yh + k*yk + l*yl) / ny;
@@ -299,7 +300,7 @@ static void render_za(UnitCell *cell, RefList *list,
                       int colscale,
                       signed int xh, signed int xk, signed int xl,
                       signed int yh, signed int yk, signed int yl,
-                      const char *outfile, double scale_top)
+                      const char *outfile, double scale_top, signed int zone)
 {
 	cairo_surface_t *surface;
 	cairo_t *dctx;
@@ -415,7 +416,7 @@ static void render_za(UnitCell *cell, RefList *list,
 	draw_circles(xh, xk, xl, yh, yk, yl, zh, zk, zl,
 	             list, sym, dctx, wght, boost, colscale, cell,
 	             max_r, theta, as, bs, cx, cy, scale,
-	             max_val);
+	             max_val, zone);
 
 	/* Centre marker */
 	cairo_arc(dctx, (double)cx,
@@ -609,6 +610,7 @@ int main(int argc, char *argv[])
 	char *outfile = NULL;
 	double scale_top = -1.0;
 	char *endptr;
+	long int zone = 0;
 
 	/* Long options */
 	const struct option longopts[] = {
@@ -625,6 +627,7 @@ int main(int argc, char *argv[])
 		{"counts",             0, &config_sqrt,        1},
 		{"colour-key",         0, &config_colkey,      1},
 		{"scale-top",          1, NULL,                2},
+		{"zone",               1, NULL,                3},
 		{0, 0, NULL, 0}
 	};
 
@@ -685,6 +688,17 @@ int main(int argc, char *argv[])
 			}
 			break;
 
+			case 3 :
+			errno = 0;
+			zone = strtol(optarg, &endptr, 10);
+			if ( !( (optarg[0] != '\0') && (endptr[0] == '\0') )
+			   || (errno != 0) )
+			{
+				ERROR("Invalid zone number ('%s')\n", optarg);
+				return 1;
+			}
+			break;
+
 			case 0 :
 			break;
 
@@ -804,7 +818,7 @@ int main(int argc, char *argv[])
 	}
 
 	render_za(cell, list, boost, sym, wght, colscale,
-	          rh, rk, rl, dh, dk, dl, outfile, scale_top);
+	          rh, rk, rl, dh, dk, dl, outfile, scale_top, zone);
 
 	free(pdb);
 	free_symoplist(sym);
diff --git a/tests/.gitignore b/tests/.gitignore
index c2a45f5a..81fd558f 100644
--- a/tests/.gitignore
+++ b/tests/.gitignore
@@ -5,4 +5,6 @@ gpu_sim_check
 integration_check
 pr_gradient_check
 symmetry_check
+centering_check
+transformation_check
 .dirstamp
diff --git a/tests/centering_check.c b/tests/centering_check.c
new file mode 100644
index 00000000..e17915d3
--- /dev/null
+++ b/tests/centering_check.c
@@ -0,0 +1,296 @@
+/*
+ * centering_check.c
+ *
+ * Check that centering of cells works
+ *
+ * Copyright © 2012 Thomas White <taw@physics.org>
+ *
+ * This file is part of CrystFEL.
+ *
+ * CrystFEL is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU General Public License as published by
+ * the Free Software Foundation, either version 3 of the License, or
+ * (at your option) any later version.
+ *
+ * CrystFEL is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU General Public License
+ * along with CrystFEL.  If not, see <http://www.gnu.org/licenses/>.
+ *
+ */
+
+#ifdef HAVE_CONFIG_H
+#include <config.h>
+#endif
+
+
+#include <stdlib.h>
+#include <stdio.h>
+#include <stdarg.h>
+#include <fenv.h>
+
+#include <cell.h>
+#include <cell-utils.h>
+
+
+static int check_cell(UnitCell *cell, const char *text)
+{
+	int err = 0;
+
+	err = validate_cell(cell);
+
+	if ( err ) {
+		ERROR("%s cell:\n", text);
+		cell_print(cell);
+	}
+
+	return err;
+}
+
+
+static int check_centering(double a, double b, double c,
+                           double al, double be, double ga,
+                           LatticeType latt, char cen, char ua)
+{
+	UnitCell *cell, *cref;
+	UnitCell *n;
+	UnitCellTransformation *t;
+	int fail = 0;
+	int i;
+	double asx, asy, asz;
+	double bsx, bsy, bsz;
+	double csx, csy, csz;
+	double ax, ay, az;
+	double bx, by, bz;
+	double cx, cy, cz;
+
+	STATUS("   --------------->  "
+	       "Checking %s %c (ua %c) %5.2e %5.2e %5.2e %5.2f %5.2f %5.2f\n",
+	       str_lattice(latt), cen, ua, a, b, c, al, be, ga);
+
+	cref = cell_new_from_parameters(a, b, c,
+	                                deg2rad(al), deg2rad(be), deg2rad(ga));
+	cell_set_lattice_type(cref, latt);
+	cell_set_centering(cref, cen);
+	cell_set_unique_axis(cref, ua);
+
+	cell = cell_rotate(cref, random_quaternion());
+	if ( cell == NULL ) return 1;
+	cell_free(cref);
+
+	check_cell(cell, "Input");
+	n = uncenter_cell(cell, &t);
+	if ( n != NULL ) {
+		STATUS("Transformation was:\n");
+		tfn_print(t);
+		if ( check_cell(n, "Output") ) fail = 1;
+		if ( !fail ) cell_print(n);
+	} else {
+		fail = 1;
+	}
+
+	cell_get_reciprocal(cell, &asx, &asy, &asz,
+	                          &bsx, &bsy, &bsz,
+	                          &csx, &csy, &csz);
+	cell_get_cartesian(n, &ax, &ay, &az,
+	                      &bx, &by, &bz,
+	                      &cx, &cy, &cz);
+
+	fesetround(1);  /* Round towards nearest */
+	for ( i=0; i<100; i++ ) {
+
+		signed int h, k, l;
+		double x, y, z;
+		double nh, nk, nl;
+		double dh, dk, dl;
+		int f = 0;
+
+		do {
+
+			h = flat_noise(0, 30);
+			k = flat_noise(0, 30);
+			l = flat_noise(0, 30);
+
+		} while ( forbidden_reflection(cell, h, k, l) );
+
+		x = h*asx + k*bsx + l*csx;
+		y = h*asy + k*bsy + l*csy;
+		z = h*asz + k*bsz + l*csz;
+
+		nh = x*ax + y*ay + z*az;
+		nk = x*bx + y*by + z*bz;
+		nl = x*cx + y*cy + z*cz;
+
+		dh = nh - lrint(nh);
+		dk = nk - lrint(nk);
+		dl = nl - lrint(nl);
+		if ( fabs(dh) > 0.1 ) f++;
+		if ( fabs(dk) > 0.1 ) f++;
+		if ( fabs(dl) > 0.1 ) f++;
+
+		if ( f ) {
+			STATUS("Centered %3i %3i %3i -> "
+			       "Primitive %7.2f %7.2f %7.2f\n",
+			       h, k, l, nh, nk, nl);
+			fail = 1;
+		}
+
+	}
+
+	cell_get_reciprocal(n, &asx, &asy, &asz,
+	                       &bsx, &bsy, &bsz,
+	                       &csx, &csy, &csz);
+	cell_get_cartesian(cell, &ax, &ay, &az,
+	                         &bx, &by, &bz,
+	                         &cx, &cy, &cz);
+
+	for ( i=0; i<100; i++ ) {
+
+		signed int h, k, l;
+		double x, y, z;
+		double nh, nk, nl;
+		double dh, dk, dl;
+		int f = 0;
+		long int ih, ik, il;
+
+		h = flat_noise(0, 30);
+		k = flat_noise(0, 30);
+		l = flat_noise(0, 30);
+
+		x = h*asx + k*bsx + l*csx;
+		y = h*asy + k*bsy + l*csy;
+		z = h*asz + k*bsz + l*csz;
+
+		nh = x*ax + y*ay + z*az;
+		nk = x*bx + y*by + z*bz;
+		nl = x*cx + y*cy + z*cz;
+
+		dh = nh - lrint(nh);  dk = nk - lrint(nk);  dl = nl - lrint(nl);
+
+		if ( fabs(dh) > 0.1 ) f++;
+		if ( fabs(dk) > 0.1 ) f++;
+		if ( fabs(dl) > 0.1 ) f++;
+
+		ih = lrint(nh);  ik = lrint(nk);  il = lrint(nl);
+		if ( forbidden_reflection(cell, ih, ik, il) ) {
+			STATUS("Primitive %3i %3i %3i -> "
+			       "Centered %3li %3li %3li, "
+			       "which is forbidden\n",
+			       h, k, l, ih, ik, il);
+			fail = 1;
+		}
+
+		if ( f ) {
+			STATUS("Primitive %3i %3i %3i -> "
+			       "Centered %7.2f %7.2f %7.2f\n",
+			       h, k, l, nh, nk, nl);
+			fail = 1;
+		}
+
+	}
+
+	return fail;
+}
+
+
+int main(int argc, char *argv[])
+{
+	int fail = 0;
+
+	/* Triclinic P */
+	fail += check_centering(50e-10, 55e-10, 70e-10, 67.0, 70.0, 77.0,
+	                        L_TRICLINIC, 'P', '*');
+
+	/* Monoclinic P */
+	fail += check_centering(10e-10, 20e-10, 30e-10, 100.0, 90.0, 90.0,
+	                        L_MONOCLINIC, 'P', 'a');
+	fail += check_centering(10e-10, 20e-10, 30e-10, 90.0, 100.0, 90.0,
+	                        L_MONOCLINIC, 'P', 'b');
+	fail += check_centering(10e-10, 20e-10, 30e-10, 90.0, 90.0, 100.0,
+	                        L_MONOCLINIC, 'P', 'c');
+
+	/* Monoclinic A */
+	fail += check_centering(10e-10, 20e-10, 30e-10, 100.0, 90.0, 90.0,
+	                        L_MONOCLINIC, 'A', 'a');
+
+	/* Monoclinic B */
+	fail += check_centering(10e-10, 20e-10, 30e-10, 90.0, 100.0, 90.0,
+	                        L_MONOCLINIC, 'B', 'b');
+
+	/* Monoclinic C */
+	fail += check_centering(10e-10, 20e-10, 30e-10, 90.0, 90.0, 100.0,
+	                        L_MONOCLINIC, 'C', 'c');
+
+	/* Orthorhombic P */
+	fail += check_centering(10e-10, 20e-10, 30e-10, 90.0, 90.0, 90.0,
+	                        L_ORTHORHOMBIC, 'P', '*');
+
+	/* Orthorhombic A */
+	fail += check_centering(10e-10, 20e-10, 30e-10, 90.0, 90.0, 90.0,
+	                        L_ORTHORHOMBIC, 'A', 'a');
+
+	/* Orthorhombic B */
+	fail += check_centering(10e-10, 20e-10, 30e-10, 90.0, 90.0, 90.0,
+	                        L_ORTHORHOMBIC, 'B', 'b');
+
+	/* Orthorhombic C */
+	fail += check_centering(10e-10, 20e-10, 30e-10, 90.0, 90.0, 90.0,
+	                        L_ORTHORHOMBIC, 'C', 'c');
+
+	/* Orthorhombic I */
+	fail += check_centering(10e-10, 20e-10, 30e-10, 90.0, 90.0, 90.0,
+	                        L_ORTHORHOMBIC, 'I', '*');
+
+	/* Orthorhombic F */
+	fail += check_centering(10e-10, 20e-10, 30e-10, 90.0, 90.0, 90.0,
+	                        L_ORTHORHOMBIC, 'F', '*');
+
+	/* Tetragonal P */
+	fail += check_centering(10e-10, 20e-10, 30e-10, 90.0, 90.0, 90.0,
+	                        L_TETRAGONAL, 'P', 'a');
+	fail += check_centering(10e-10, 20e-10, 30e-10, 90.0, 90.0, 90.0,
+	                        L_TETRAGONAL, 'P', 'b');
+	fail += check_centering(10e-10, 20e-10, 30e-10, 90.0, 90.0, 90.0,
+	                        L_TETRAGONAL, 'P', 'c');
+
+	/* Tetragonal I */
+	fail += check_centering(10e-10, 20e-10, 30e-10, 90.0, 90.0, 90.0,
+	                        L_TETRAGONAL, 'I', 'a');
+	fail += check_centering(10e-10, 20e-10, 30e-10, 90.0, 90.0, 90.0,
+	                        L_TETRAGONAL, 'I', 'b');
+	fail += check_centering(10e-10, 20e-10, 30e-10, 90.0, 90.0, 90.0,
+	                        L_TETRAGONAL, 'I', 'c');
+
+	/* Rhombohedral R */
+	fail += check_centering(10e-10, 10e-10, 10e-10, 60.0, 60.0, 60.0,
+	                        L_RHOMBOHEDRAL, 'R', '*');
+
+	/* Hexagonal P */
+	fail += check_centering(10e-10, 20e-10, 30e-10, 120.0, 90.0, 90.0,
+	                        L_HEXAGONAL, 'P', 'a');
+	fail += check_centering(10e-10, 20e-10, 30e-10, 90.0, 120.0, 90.0,
+	                        L_HEXAGONAL, 'P', 'b');
+	fail += check_centering(10e-10, 20e-10, 30e-10, 90.0, 90.0, 120.0,
+	                        L_HEXAGONAL, 'P', 'c');
+
+	/* Hexagonal H (PDB-speak for rhombohedral) */
+	fail += check_centering(20e-10, 20e-10, 40e-10, 90.0, 90.0, 120.0,
+	                        L_HEXAGONAL, 'H', 'c');
+
+	/* Cubic P */
+	fail += check_centering(30e-10, 30e-10, 30e-10, 90.0, 90.0, 90.0,
+	                        L_CUBIC, 'P', '*');
+
+	/* Cubic I */
+	fail += check_centering(30e-10, 30e-10, 30e-10, 90.0, 90.0, 90.0,
+	                        L_CUBIC, 'I', '*');
+
+	/* Cubic F */
+	fail += check_centering(30e-10, 30e-10, 30e-10, 90.0, 90.0, 90.0,
+	                        L_CUBIC, 'F', '*');
+
+	return fail;
+}
diff --git a/tests/gpu_sim_check.c b/tests/gpu_sim_check.c
index 5b6f6da7..5c45ea02 100644
--- a/tests/gpu_sim_check.c
+++ b/tests/gpu_sim_check.c
@@ -36,6 +36,7 @@
 #include <beam-parameters.h>
 #include <utils.h>
 #include <symmetry.h>
+#include <cell-utils.h>
 
 
 #ifdef HAVE_CLOCK_GETTIME
diff --git a/tests/pr_gradient_check.c b/tests/pr_gradient_check.c
index 6972b66a..b39dfd05 100644
--- a/tests/pr_gradient_check.c
+++ b/tests/pr_gradient_check.c
@@ -32,6 +32,7 @@
 
 #include <image.h>
 #include <cell.h>
+#include <cell-utils.h>
 #include <geometry.h>
 #include <reflist.h>
 #include "../src/post-refinement.h"
diff --git a/tests/transformation_check.c b/tests/transformation_check.c
new file mode 100644
index 00000000..7d25aa04
--- /dev/null
+++ b/tests/transformation_check.c
@@ -0,0 +1,176 @@
+/*
+ * transformation_check.c
+ *
+ * Check that unit cell transformations work
+ *
+ * Copyright © 2012 Thomas White <taw@physics.org>
+ *
+ * This file is part of CrystFEL.
+ *
+ * CrystFEL is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU General Public License as published by
+ * the Free Software Foundation, either version 3 of the License, or
+ * (at your option) any later version.
+ *
+ * CrystFEL is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU General Public License
+ * along with CrystFEL.  If not, see <http://www.gnu.org/licenses/>.
+ *
+ */
+
+#ifdef HAVE_CONFIG_H
+#include <config.h>
+#endif
+
+
+#include <stdlib.h>
+#include <stdio.h>
+#include <stdarg.h>
+
+#include <cell.h>
+#include <cell-utils.h>
+
+
+static int check_transformation(UnitCell *cell, UnitCellTransformation *tfn)
+{
+	UnitCell *cnew, *cback;
+	UnitCellTransformation *inv;
+	double a[9], b[9];
+	int i;
+	int fail = 0;
+
+	cnew = cell_transform(cell, tfn);
+
+	cback = cell_transform_inverse(cnew, tfn);
+
+	cell_get_cartesian(cell, &a[0], &a[1], &a[2],
+	                         &a[3], &a[4], &a[5],
+	                         &a[6], &a[7], &a[8]);
+	cell_get_cartesian(cback, &b[0], &b[1], &b[2],
+	                          &b[3], &b[4], &b[5],
+	                          &b[6], &b[7], &b[8]);
+	for ( i=0; i<9; i++ ) {
+		if ( !within_tolerance(a[i], b[i], 0.1) ) {
+			fail = 1;
+			STATUS("%e %e\n", a[i], b[i]);
+		}
+	}
+
+	if ( fail ) {
+		ERROR("Original cell not recovered after transformation:\n");
+		cell_print(cell);
+		tfn_print(tfn);
+		inv = tfn_inverse(tfn);
+		tfn_print(inv);
+		cell_print(cback);
+	}
+
+	return fail;
+}
+
+
+static int check_identity(UnitCell *cell, UnitCellTransformation *tfn)
+{
+	UnitCell *cnew;
+	double a[9], b[9];
+	int i;
+	int fail = 0;
+
+	cnew = cell_transform(cell, tfn);
+
+	cell_get_cartesian(cell, &a[0], &a[1], &a[2],
+	                         &a[3], &a[4], &a[5],
+	                         &a[6], &a[7], &a[8]);
+	cell_get_cartesian(cnew, &b[0], &b[1], &b[2],
+	                         &b[3], &b[4], &b[5],
+	                         &b[6], &b[7], &b[8]);
+	for ( i=0; i<9; i++ ) {
+		if ( !within_tolerance(a[i], b[i], 0.1) ) {
+			fail = 1;
+			STATUS("%e %e\n", a[i], b[i]);
+		}
+	}
+
+	if ( fail ) {
+		ERROR("Original cell not recovered after transformation:\n");
+		cell_print(cell);
+		tfn_print(tfn);
+		cell_print(cnew);
+	}
+
+	return fail;
+}
+
+
+int main(int argc, char *argv[])
+{
+	int fail = 0;
+	UnitCell *cell, *cref;
+	UnitCellTransformation *tfn;
+
+	cref = cell_new_from_parameters(50e-10, 55e-10, 70e-10,
+	                                deg2rad(67.0),
+	                                deg2rad(70.0),
+	                                deg2rad(77.0));
+	if ( cref == NULL ) return 1;
+
+	cell = cell_rotate(cref, random_quaternion());
+	if ( cell == NULL ) return 1;
+	cell_free(cref);
+
+	/* Permutation of axes */
+	tfn = tfn_identity();
+	if ( tfn == NULL ) return 1;
+	tfn_combine(tfn, tfn_vector(0,1,0),
+	                 tfn_vector(0,0,1),
+	                 tfn_vector(1,0,0));
+	fail += check_transformation(cell, tfn);
+	tfn_free(tfn);
+
+	/* Doubling of cell in one direction */
+	tfn = tfn_identity();
+	if ( tfn == NULL ) return 1;
+	tfn_combine(tfn, tfn_vector(2,0,0),
+	                 tfn_vector(0,1,0),
+	                 tfn_vector(0,0,1));
+	fail += check_transformation(cell, tfn);
+	tfn_free(tfn);
+
+	/* Diagonal supercell */
+	tfn = tfn_identity();
+	if ( tfn == NULL ) return 1;
+	tfn_combine(tfn, tfn_vector(1,1,0),
+	                 tfn_vector(0,1,0),
+	                 tfn_vector(0,0,1));
+	fail += check_transformation(cell, tfn);
+	tfn_free(tfn);
+
+	/* Crazy */
+	tfn = tfn_identity();
+	if ( tfn == NULL ) return 1;
+	tfn_combine(tfn, tfn_vector(1,1,0),
+	                 tfn_vector(0,1,0),
+	                 tfn_vector(0,1,1));
+	fail += check_transformation(cell, tfn);
+	tfn_free(tfn);
+
+	/* Identity in two parts */
+	tfn = tfn_identity();
+	if ( tfn == NULL ) return 1;
+	tfn_combine(tfn, tfn_vector(0,0,1),
+	                 tfn_vector(0,1,0),
+	                 tfn_vector(-1,0,0));
+	tfn_combine(tfn,  tfn_vector(0,0,-1),
+	                  tfn_vector(0,1,0),
+	                  tfn_vector(1,0,0));
+	fail += check_identity(cell, tfn);
+	tfn_free(tfn);
+
+	cell_free(cell);
+
+	return fail;
+}