Merge branch 'tom/transformations'

16 files changed, 2156 insertions, 680 deletions

diff --git a/libcrystfel/CMakeLists.txt b/libcrystfel/CMakeLists.txt
index 247b925e..014d2a5a 100644
--- a/libcrystfel/CMakeLists.txt
+++ b/libcrystfel/CMakeLists.txt
@@ -7,6 +7,8 @@ find_package(PINKINDEXER)
 find_package(NBP)
 find_package(FDIP)
 find_package(ZLIB REQUIRED)
+find_package(FLEX REQUIRED)
+find_package(BISON REQUIRED)
 pkg_search_module(FFTW fftw3)
 
 set(HAVE_CURSES ${CURSES_FOUND})
@@ -30,6 +32,12 @@ endif()
 
 configure_file(config.h.cmake.in config.h)
 
+bison_target(symopp src/symop.y ${CMAKE_CURRENT_BINARY_DIR}/symop-parse.c COMPILE_FLAGS --report=all)
+flex_target(symopl src/symop.l ${CMAKE_CURRENT_BINARY_DIR}/symop-lex.c
+	DEFINES_FILE ${CMAKE_CURRENT_BINARY_DIR}/symop-lex.h)
+add_flex_bison_dependency(symopl symopp)
+include_directories(${PROJECT_SOURCE_DIR}/src)
+
 set(LIBCRYSTFEL_SOURCES
     src/reflist.c
     src/utils.c
@@ -62,6 +70,9 @@ set(LIBCRYSTFEL_SOURCES
     src/peakfinder8.c
     src/taketwo.c
     src/xgandalf.c
+    src/rational.c
+    ${BISON_symopp_OUTPUTS}
+    ${FLEX_symopl_OUTPUTS}
 )
 
 if (HAVE_FFTW)
@@ -101,6 +112,7 @@ set(LIBCRYSTFEL_HEADERS
     src/peakfinder8.h
     src/taketwo.h
     src/xgandalf.h
+    src/rational.h
 )
 
 add_library(${PROJECT_NAME} SHARED
diff --git a/libcrystfel/src/cell-utils.c b/libcrystfel/src/cell-utils.c
index 7b1984bb..a36ebe1c 100644
--- a/libcrystfel/src/cell-utils.c
+++ b/libcrystfel/src/cell-utils.c
@@ -3,13 +3,13 @@
  *
  * Unit Cell utility functions
  *
- * Copyright © 2012-2017 Deutsches Elektronen-Synchrotron DESY,
+ * Copyright © 2012-2019 Deutsches Elektronen-Synchrotron DESY,
  *                       a research centre of the Helmholtz Association.
  * Copyright © 2012 Lorenzo Galli
  *
  * Authors:
- *   2009-2012,2014-2017 Thomas White <taw@physics.org>
- *   2012                Lorenzo Galli
+ *   2009-2019 Thomas White <taw@physics.org>
+ *   2012      Lorenzo Galli
  *
  * This file is part of CrystFEL.
  *
@@ -219,11 +219,6 @@ int right_handed(UnitCell *cell)
 
 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;
 
@@ -251,12 +246,31 @@ void cell_print(UnitCell *cell)
 	}
 
 	if ( cell_has_parameters(cell) ) {
+
+		double a, b, c, alpha, beta, gamma;
 		cell_get_parameters(cell, &a, &b, &c, &alpha, &beta, &gamma);
 
 		STATUS("a      b      c            alpha   beta  gamma\n");
 		STATUS("%6.2f %6.2f %6.2f A    %6.2f %6.2f %6.2f deg\n",
 		       a*1e10, b*1e10, c*1e10,
 		       rad2deg(alpha), rad2deg(beta), rad2deg(gamma));
+	} else {
+		STATUS("Unit cell parameters are not specified.\n");
+	}
+}
+
+
+void cell_print_full(UnitCell *cell)
+{
+
+	cell_print(cell);
+
+	if ( cell_has_parameters(cell) ) {
+
+		double asx, asy, asz;
+		double bsx, bsy, bsz;
+		double csx, csy, csz;
+		double ax, ay, az, bx, by, bz, cx, cy, cz;
 
 		cell_get_cartesian(cell, &ax, &ay, &az, &bx, &by, &bz, &cx, &cy, &cz);
 
@@ -283,8 +297,6 @@ void cell_print(UnitCell *cell)
 
 		STATUS("Cell representation is %s.\n", cell_rep(cell));
 
-	} else {
-		STATUS("Unit cell parameters are not specified.\n");
 	}
 }
 
@@ -349,203 +361,218 @@ int bravais_lattice(UnitCell *cell)
 }
 
 
-static UnitCellTransformation *uncentering_transformation(UnitCell *in,
-                                                          char *new_centering,
-                                                          LatticeType *new_latt)
+static RationalMatrix *create_rtnl_mtx(signed int a1, signed int a2,
+                                       signed int b1, signed int b2,
+                                       signed int c1, signed int c2,
+                                       signed int d1, signed int d2,
+                                       signed int e1, signed int e2,
+                                       signed int f1, signed int f2,
+                                       signed int g1, signed int g2,
+                                       signed int h1, signed int h2,
+                                       signed int i1, signed int i2)
+{
+	RationalMatrix *m = rtnl_mtx_new(3, 3);
+	if ( m == NULL ) return NULL;
+	rtnl_mtx_set(m, 0, 0, rtnl(a1, a2));
+	rtnl_mtx_set(m, 0, 1, rtnl(b1, b2));
+	rtnl_mtx_set(m, 0, 2, rtnl(c1, c2));
+	rtnl_mtx_set(m, 1, 0, rtnl(d1, d2));
+	rtnl_mtx_set(m, 1, 1, rtnl(e1, e2));
+	rtnl_mtx_set(m, 1, 2, rtnl(f1, f2));
+	rtnl_mtx_set(m, 2, 0, rtnl(g1, g2));
+	rtnl_mtx_set(m, 2, 1, rtnl(h1, h2));
+	rtnl_mtx_set(m, 2, 2, rtnl(i1, i2));
+	return m;
+}
+
+
+/* Given a centered cell @in, return the integer transformation matrix which
+ * turns a primitive cell into @in. Set new_centering and new_latt to the
+ * centering and lattice type of the primitive cell (usually aP, sometimes rR,
+ * rarely mP).  Store the inverse matrix at pCi */
+static IntegerMatrix *centering_transformation(UnitCell *in,
+                                               char *new_centering,
+                                               LatticeType *new_latt,
+                                               char *new_ua,
+                                               RationalMatrix **pCi)
 {
-	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;
+	IntegerMatrix *C = NULL;
+	RationalMatrix *Ci = NULL;
 
 	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,1,0),
-		               tfn_vector(0,0,1),
-		               tfn_vector(1,0,0));
-		if ( lt == L_MONOCLINIC ) {
-			assert(cen != 'A');
-			switch ( cen ) {
-				case 'B' : cen = 'A'; break;
-				case 'C' : cen = 'B'; break;
-				case 'I' : cen = 'I'; break;
-			}
-		}
-	}
-
-	if ( ua == 'b' ) {
-		tfn_combine(t, tfn_vector(0,0,1),
-		               tfn_vector(1,0,0),
-		               tfn_vector(0,1,0));
-		if ( lt == L_MONOCLINIC ) {
-			assert(cen != 'B');
-			switch ( cen ) {
-				case 'C' : cen = 'A'; break;
-				case 'A' : cen = 'B'; break;
-				case 'I' : cen = 'I'; break;
-			}
-		}
-	}
-
-	switch ( cen ) {
+	/* Write the matrices exactly as they appear in ITA Table 5.1.3.1.
+	 * C is "P", and Ci is "Q=P^-1".  Vice-versa if the transformation
+	 * should go the opposite way to what's written in the first column. */
 
-		case 'P' :
+	if ( (cen=='P') || (cen=='R') ) {
+		*new_centering = cen;
 		*new_latt = lt;
-		*new_centering = 'P';
-		break;
-
-		case 'R' :
-		*new_latt = L_RHOMBOHEDRAL;
-		*new_centering = 'R';
-		break;
+		*new_ua = ua;
+		C = intmat_identity(3);
+		Ci = rtnl_mtx_identity(3);
+	}
 
-		case 'I' :
-		tfn_combine(t, tfn_vector(-H,H,H),
-		               tfn_vector(H,-H,H),
-		               tfn_vector(H,H,-H));
+	if ( cen == 'I' ) {
+		C = intmat_create_3x3(0, 1, 1,
+		                      1, 0, 1,
+		                      1, 1, 0);
+		Ci = create_rtnl_mtx(-1,2,  1,2,  1,2,
+		                      1,2, -1,2,  1,2,
+		                      1,2,  1,2, -1,2);
 		if ( lt == L_CUBIC ) {
 			*new_latt = L_RHOMBOHEDRAL;
 			*new_centering = 'R';
+			*new_ua = '*';
 		} else {
-			/* Tetragonal or orthorhombic */
 			*new_latt = L_TRICLINIC;
 			*new_centering = 'P';
+			*new_ua = '*';
 		}
-		break;
+	}
 
-		case 'F' :
-		tfn_combine(t, tfn_vector(0,H,H),
-		               tfn_vector(H,0,H),
-		               tfn_vector(H,H,0));
+	if ( cen == 'F' ) {
+		C = intmat_create_3x3(-1,  1,  1,
+		                       1, -1,  1,
+		                       1,  1, -1);
+		Ci = create_rtnl_mtx( 0,1,  1,2,  1,2,
+		                      1,2,  0,1,  1,2,
+		                      1,2,  1,2,  0,1);
 		if ( lt == L_CUBIC ) {
 			*new_latt = L_RHOMBOHEDRAL;
 			*new_centering = 'R';
+			*new_ua = '*';
 		} else {
-			assert(lt == L_ORTHORHOMBIC);
 			*new_latt = L_TRICLINIC;
 			*new_centering = 'P';
+			*new_ua = '*';
 		}
-		break;
+	}
 
-		case 'A' :
-		tfn_combine(t, tfn_vector( 1, 0, 0),
-		               tfn_vector( 0, H, H),
-		               tfn_vector( 0,-H, H));
+	if ( (lt == L_HEXAGONAL) && (cen == 'H') && (ua == 'c') ) {
+		/* Obverse setting */
+		C = intmat_create_3x3( 1,  0,  1,
+		                      -1,  1,  1,
+		                       0, -1,  1);
+		Ci = create_rtnl_mtx( 2,3, -1,3, -1,3,
+		                      1,3,  1,3, -2,3,
+		                      1,3,  1,3,  1,3);
+		assert(lt == L_HEXAGONAL);
+		assert(ua == 'c');
+		*new_latt = L_RHOMBOHEDRAL;
+		*new_centering = 'R';
+		*new_ua = '*';
+	}
+
+	if ( cen == 'A' ) {
+		C = intmat_create_3x3( 1,  0,  0,
+		                       0,  1,  1,
+		                       0, -1,  1);
+		Ci = create_rtnl_mtx( 1,1,  0,1,  0,1,
+		                      0,1,  1,2, -1,2,
+		                      0,1,  1,2,  1,2);
 		if ( lt == L_ORTHORHOMBIC ) {
 			*new_latt = L_MONOCLINIC;
+			*new_centering = 'P';
+			*new_ua = 'a';
 		} else {
 			*new_latt = L_TRICLINIC;
+			*new_centering = 'P';
+			*new_ua = '*';
 		}
-		*new_centering = 'P';
-		break;
+	}
 
-		case 'B' :
-		tfn_combine(t, tfn_vector( H, 0, H),
-		               tfn_vector( 0, 1, 0),
-		               tfn_vector(-H, 0, H));
+	if ( cen == 'B' ) {
+		C = intmat_create_3x3( 1,  0,  1,
+		                       0,  1,  0,
+		                      -1,  0,  1);
+		Ci = create_rtnl_mtx( 1,2,  0,1, -1,2,
+		                      0,1,  1,1,  0,1,
+		                      1,2,  0,1,  1,2);
 		if ( lt == L_ORTHORHOMBIC ) {
 			*new_latt = L_MONOCLINIC;
+			*new_centering = 'P';
+			*new_ua = 'b';
 		} else {
 			*new_latt = L_TRICLINIC;
+			*new_centering = 'P';
+			*new_ua = '*';
 		}
-		*new_centering = 'P';
-		break;
+	}
 
-		case 'C' :
-		tfn_combine(t, tfn_vector( H, H, 0),
-		               tfn_vector(-H, H, 0),
-		               tfn_vector( 0, 0, 1));
+	if ( cen == 'C' ) {
+		C = intmat_create_3x3( 1,  1,  0,
+		                      -1,  1,  0,
+		                       0,  0,  1);
+		Ci = create_rtnl_mtx( 1,2, -1,2,  0,1,
+		                      1,2,  1,2,  0,1,
+		                      0,1,  0,1,  1,1);
 		if ( lt == L_ORTHORHOMBIC ) {
 			*new_latt = L_MONOCLINIC;
+			*new_centering = 'P';
+			*new_ua = 'c';
 		} else {
 			*new_latt = L_TRICLINIC;
-		}
-		*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(1,0,0),
-				       tfn_vector(0,1,0));
-		}
-
-		if ( ua == 'b' ) {
-			tfn_combine(t, tfn_vector(0,1,0),
-				       tfn_vector(0,0,1),
-				       tfn_vector(1,0,0));
+			*new_centering = 'P';
+			*new_ua = '*';
 		}
 	}
 
-	return t;
+	*pCi = Ci;
+	return C;
 }
 
 
 /**
  * uncenter_cell:
  * @in: A %UnitCell
- * @t: Location at which to store the transformation which was used.
+ * @C: Location at which to store the centering transformation
+ * @Ci: Location at which to store the inverse centering transformation
  *
- * 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.
+ * Turns any cell into a primitive one, e.g. for comparison purposes.
  *
- * Returns: a primitive version of @in in a conventional (unique axis c)
- * setting.
+ * The transformation which was used is stored at @Ci. The centering
+ * transformation, which is the transformation you should apply if you want to
+ * get back the original cell, will be stored at @C.  Either or both of these
+ * can be NULL if you don't need that information.
+ *
+ * Returns: a primitive version of @in.
  *
  */
-UnitCell *uncenter_cell(UnitCell *in, UnitCellTransformation **t)
+UnitCell *uncenter_cell(UnitCell *in, IntegerMatrix **pC, RationalMatrix **pCi)
 {
-	UnitCellTransformation *tt;
+	IntegerMatrix *C;
+	RationalMatrix *Ci;
 	char new_centering;
 	LatticeType new_latt;
+	char new_ua;
 	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;
+	C = centering_transformation(in, &new_centering, &new_latt,
+	                             &new_ua, &Ci);
+	if ( C == NULL ) return NULL;
 
-	out = cell_transform(in, tt);
+	out = cell_transform_rational(in, Ci);
 	if ( out == NULL ) return NULL;
 
 	cell_set_lattice_type(out, new_latt);
 	cell_set_centering(out, new_centering);
+	cell_set_unique_axis(out, new_ua);
+
+	if ( pC != NULL ) {
+		*pC = C;
+	} else {
+		intmat_free(C);
+	}
 
-	if ( t != NULL ) {
-		*t = tt;
+	if ( pCi != NULL ) {
+		*pCi = Ci;
 	} else {
-		tfn_free(tt);
+		rtnl_mtx_free(Ci);
 	}
 
 	return out;
@@ -609,16 +636,16 @@ UnitCell *match_cell(UnitCell *cell_in, UnitCell *template_in, int verbose,
 	float angtol = deg2rad(tols[3]);
 	UnitCell *cell;
 	UnitCell *template;
-	UnitCellTransformation *uncentering;
+	IntegerMatrix *centering;
 	UnitCell *new_cell_trans;
 
 	/* "Un-center" the template unit cell to make the comparison easier */
-	template = uncenter_cell(template_in, &uncentering);
+	template = uncenter_cell(template_in, &centering, NULL);
 	if ( template == NULL ) return NULL;
 
 	/* The candidate cell is also uncentered, because it might be centered
 	 * if it came from (e.g.) MOSFLM */
-	cell = uncenter_cell(cell_in, NULL);
+	cell = uncenter_cell(cell_in, NULL, NULL);
 	if ( cell == NULL ) return NULL;
 
 	if ( cell_get_reciprocal(template, &asx, &asy, &asz,
@@ -627,7 +654,7 @@ UnitCell *match_cell(UnitCell *cell_in, UnitCell *template_in, int verbose,
 		ERROR("Couldn't get reciprocal cell for template.\n");
 		cell_free(template);
 		cell_free(cell);
-		tfn_free(uncentering);
+		intmat_free(centering);
 		return NULL;
 	}
 
@@ -649,7 +676,7 @@ UnitCell *match_cell(UnitCell *cell_in, UnitCell *template_in, int verbose,
 		ERROR("Couldn't get reciprocal cell.\n");
 		cell_free(template);
 		cell_free(cell);
-		tfn_free(uncentering);
+		intmat_free(centering);
 		return NULL;
 	}
 
@@ -811,7 +838,7 @@ UnitCell *match_cell(UnitCell *cell_in, UnitCell *template_in, int verbose,
 	/* Reverse the de-centering transformation */
 	if ( new_cell != NULL ) {
 
-		new_cell_trans = cell_transform_inverse(new_cell, uncentering);
+		new_cell_trans = cell_transform_intmat(new_cell, centering);
 		cell_free(new_cell);
 		cell_set_lattice_type(new_cell_trans,
 		                      cell_get_lattice_type(template_in));
@@ -821,13 +848,13 @@ UnitCell *match_cell(UnitCell *cell_in, UnitCell *template_in, int verbose,
 		                     cell_get_unique_axis(template_in));
 
 		cell_free(template);
-		tfn_free(uncentering);
+		intmat_free(centering);
 
 		return new_cell_trans;
 
 	} else {
 		cell_free(template);
-		tfn_free(uncentering);
+		intmat_free(centering);
 		return NULL;
 	}
 }
@@ -850,16 +877,16 @@ UnitCell *match_cell_ab(UnitCell *cell_in, UnitCell *template_in)
 	int have_real_c;
 	UnitCell *cell;
 	UnitCell *template;
-	UnitCellTransformation *to_given_cell;
+	IntegerMatrix *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);
+	template = uncenter_cell(template_in, &to_given_cell, NULL);
 
 	/* The candidate cell is also uncentered, because it might be centered
 	 * if it came from (e.g.) MOSFLM */
-	cell = uncenter_cell(cell_in, NULL);
+	cell = uncenter_cell(cell_in, NULL, NULL);
 
 	/* Get the lengths to match */
 	if ( cell_get_cartesian(template, &ax, &ay, &az,
@@ -940,7 +967,7 @@ UnitCell *match_cell_ab(UnitCell *cell_in, UnitCell *template_in)
 	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);
+	new_cell_trans = cell_transform_intmat_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));
@@ -1443,49 +1470,6 @@ void cell_fudge_gslcblas()
 }
 
 
-UnitCell *transform_cell_gsl(UnitCell *in, gsl_matrix *m)
-{
-	gsl_matrix *c;
-	double asx, asy, asz;
-	double bsx, bsy, bsz;
-	double csx, csy, csz;
-	gsl_matrix *res;
-	UnitCell *out;
-
-	cell_get_reciprocal(in, &asx, &asy, &asz, &bsx, &bsy,
-	                        &bsz, &csx, &csy, &csz);
-
-	c = gsl_matrix_alloc(3, 3);
-	gsl_matrix_set(c, 0, 0, asx);
-	gsl_matrix_set(c, 1, 0, asy);
-	gsl_matrix_set(c, 2, 0, asz);
-	gsl_matrix_set(c, 0, 1, bsx);
-	gsl_matrix_set(c, 1, 1, bsy);
-	gsl_matrix_set(c, 2, 1, bsz);
-	gsl_matrix_set(c, 0, 2, csx);
-	gsl_matrix_set(c, 1, 2, csy);
-	gsl_matrix_set(c, 2, 2, csz);
-
-	res = gsl_matrix_calloc(3, 3);
-	gsl_blas_dgemm(CblasNoTrans, CblasNoTrans, 1.0, m, c, 0.0, res);
-
-	out = cell_new_from_cell(in);
-	cell_set_reciprocal(out, gsl_matrix_get(res, 0, 0),
-	                         gsl_matrix_get(res, 1, 0),
-	                         gsl_matrix_get(res, 2, 0),
-	                         gsl_matrix_get(res, 0, 1),
-	                         gsl_matrix_get(res, 1, 1),
-	                         gsl_matrix_get(res, 2, 1),
-	                         gsl_matrix_get(res, 0, 2),
-	                         gsl_matrix_get(res, 1, 2),
-	                         gsl_matrix_get(res, 2, 2));
-
-	gsl_matrix_free(res);
-	gsl_matrix_free(c);
-	return out;
-}
-
-
 /**
  * rotate_cell:
  * @in: A %UnitCell to rotate
@@ -1586,7 +1570,9 @@ int cell_is_sensible(UnitCell *cell)
  * lattice is a conventional Bravais lattice.
  * Warnings are printied if any of the checks are failed.
  *
- * Returns: true if cell is invalid.
+ * Returns: zero if the cell is fine, 1 if it is unconventional but otherwise
+ *  OK (e.g. left-handed or not a Bravais lattice), and 2 if there is a serious
+ *  problem such as the parameters being physically impossible.
  *
  */
 int validate_cell(UnitCell *cell)
@@ -1596,7 +1582,7 @@ int validate_cell(UnitCell *cell)
 
 	if ( cell_has_parameters(cell) && !cell_is_sensible(cell) ) {
 		ERROR("WARNING: Unit cell parameters are not sensible.\n");
-		err = 1;
+		err = 2;
 	}
 
 	if ( !bravais_lattice(cell) ) {
@@ -1620,7 +1606,7 @@ int validate_cell(UnitCell *cell)
 		  || ((cen == 'C') && (ua == 'c')) ) {
 			ERROR("WARNING: A, B or C centering matches unique"
 			      " axis.\n");
-			err = 1;
+			err = 2;
 		}
 	}
 
@@ -1646,7 +1632,7 @@ int forbidden_reflection(UnitCell *cell,
 	cen = cell_get_centering(cell);
 
 	/* Reflection conditions here must match the transformation matrices
-	 * in uncentering_transformation().  tests/centering_check verifies
+	 * in centering_transformation().  tests/centering_check verifies
 	 * this (amongst other things). */
 
 	if ( cen == 'P' ) return 0;
@@ -1694,6 +1680,48 @@ double cell_get_volume(UnitCell *cell)
 }
 
 
+/**
+ * compare_cell_parameters:
+ * @a: A %UnitCell
+ * @b: Another %UnitCell
+ * @ltl: Maximum allowable fractional difference in axis lengths
+ * @atl: Maximum allowable difference in reciprocal angles (in radians)
+ *
+ * Compare the two unit cells.  If the real space parameters match to within
+ * fractional difference @ltl, and the inter-axial angles match within @atl,
+ * and the centering matches, this function returns 1.  Otherwise 0.
+ *
+ * This function considers the cell parameters and centering, but ignores the
+ * orientation of the cell.  If you want to compare the orientation as well,
+ * use compare_cell_parameters_and_orientation() instead.
+ *
+ * Returns: non-zero if the cells match.
+ *
+ */
+int compare_cell_parameters(UnitCell *cell1, UnitCell *cell2,
+                            float ltl, float atl)
+{
+	double a1, b1, c1, al1, be1, ga1;
+	double a2, b2, c2, al2, be2, ga2;
+
+	/* Centering must match: we don't arbitrarte primitive vs centered,
+	 * different cell choices etc */
+	if ( cell_get_centering(cell1) != cell_get_centering(cell2) ) return 0;
+
+	cell_get_parameters(cell1, &a1, &b1, &c1, &al1, &be1, &ga1);
+	cell_get_parameters(cell2, &a2, &b2, &c2, &al2, &be2, &ga2);
+
+	if ( !within_tolerance(a1, a2, ltl*100.0) ) return 0;
+	if ( !within_tolerance(b1, b2, ltl*100.0) ) return 0;
+	if ( !within_tolerance(c1, c2, ltl*100.0) ) return 0;
+	if ( fabs(al1-al2) > atl ) return 0;
+	if ( fabs(be1-be2) > atl ) return 0;
+	if ( fabs(ga1-ga2) > atl ) return 0;
+
+	return 1;
+}
+
+
 static double moduli_check(double ax, double ay, double az,
                            double bx, double by, double bz)
 {
@@ -1703,28 +1731,42 @@ static double moduli_check(double ax, double ay, double az,
 }
 
 
-static int cells_are_similar(UnitCell *cell1, UnitCell *cell2,
-                             const double ltl, const double atl)
+/**
+ * compare_cell_parameters_and_orientation:
+ * @a: A %UnitCell
+ * @b: Another %UnitCell
+ * @ltl: Maximum allowable fractional difference in reciprocal axis lengths
+ * @atl: Maximum allowable difference in reciprocal angles (in radians)
+ *
+ * Compare the two unit cells.  If the axes match in length (to within
+ * fractional difference @ltl) and the axes are aligned to within @atl radians,
+ * this function returns non-zero.
+ *
+ * This function compares the orientation of the cell as well as the parameters.
+ * If you just want to see if the parameters are the same, use
+ * compare_cell_parameters() instead.
+ *
+ * The cells @a and @b must have the same centering.  Otherwise, this function
+ * always returns zero.
+ *
+ * Returns: non-zero if the cells match.
+ *
+ */
+int compare_cell_parameters_and_orientation(UnitCell *cell1, UnitCell *cell2,
+                                            const double ltl, const double atl)
 {
 	double asx1, asy1, asz1, bsx1, bsy1, bsz1, csx1, csy1, csz1;
 	double asx2, asy2, asz2, bsx2, bsy2, bsz2, csx2, csy2, csz2;
-	UnitCell *pcell1, *pcell2;
-
-	/* Compare primitive cells, not centered */
-	pcell1 = uncenter_cell(cell1, NULL);
-	pcell2 = uncenter_cell(cell2, NULL);
-
-	cell_get_reciprocal(pcell1, &asx1, &asy1, &asz1,
-	                            &bsx1, &bsy1, &bsz1,
-	                            &csx1, &csy1, &csz1);
 
-	cell_get_reciprocal(pcell2, &asx2, &asy2, &asz2,
-	                            &bsx2, &bsy2, &bsz2,
-	                            &csx2, &csy2, &csz2);
+	if ( cell_get_centering(cell1) != cell_get_centering(cell2) ) return 0;
 
+	cell_get_cartesian(cell1, &asx1, &asy1, &asz1,
+	                          &bsx1, &bsy1, &bsz1,
+	                          &csx1, &csy1, &csz1);
 
-	cell_free(pcell1);
-	cell_free(pcell2);
+	cell_get_cartesian(cell2, &asx2, &asy2, &asz2,
+	                          &bsx2, &bsy2, &bsz2,
+	                          &csx2, &csy2, &csz2);
 
 	if ( angle_between(asx1, asy1, asz1, asx2, asy2, asz2) > atl ) return 0;
 	if ( angle_between(bsx1, bsy1, bsz1, bsx2, bsy2, bsz2) > atl ) return 0;
@@ -1738,28 +1780,38 @@ static int cells_are_similar(UnitCell *cell1, UnitCell *cell2,
 }
 
 
-
 /**
- * compare_cells:
+ * compare_reindexed_cell_parameters_and_orientation:
  * @a: A %UnitCell
  * @b: Another %UnitCell
  * @ltl: Maximum allowable fractional difference in reciprocal axis lengths
  * @atl: Maximum allowable difference in reciprocal angles (in radians)
  * @pmb: Place to store pointer to matrix
  *
- * Compare the two units cells.  If they agree to within @ltl and @atl, using
- *  any change of axes, returns non-zero and stores the transformation to map @b
- *  onto @a.
+ * Compare the two unit cells.  If, using any permutation of the axes, the
+ * axes can be made to match in length (to within fractional difference @ltl)
+ * and the axes aligned to within @atl radians, this function returns non-zero
+ * and stores the transformation to map @b onto @a.
+ *
+ * Note that the orientations of the cells must match, not just the parameters.
+ * The comparison is done after reindexing using
+ * compare_cell_parameters_and_orientation().
+ *
+ * The cells @a and @b must have the same centering.  Otherwise, this function
+ * always returns zero.
  *
  * Returns: non-zero if the cells match.
  *
  */
-int compare_cells(UnitCell *a, UnitCell *b, double ltl, double atl,
-                  IntegerMatrix **pmb)
+int compare_reindexed_cell_parameters_and_orientation(UnitCell *a, UnitCell *b,
+                                                      double ltl, double atl,
+                                                      IntegerMatrix **pmb)
 {
 	IntegerMatrix *m;
 	int i[9];
 
+	if ( cell_get_centering(a) != cell_get_centering(b) ) return 0;
+
 	m = intmat_new(3, 3);
 
 	for ( i[0]=-1; i[0]<=+1; i[0]++ ) {
@@ -1772,7 +1824,6 @@ int compare_cells(UnitCell *a, UnitCell *b, double ltl, double atl,
 	for ( i[7]=-1; i[7]<=+1; i[7]++ ) {
 	for ( i[8]=-1; i[8]<=+1; i[8]++ ) {
 
-		UnitCellTransformation *tfn;
 		UnitCell *nc;
 		int j, k;
 		int l = 0;
@@ -1783,17 +1834,14 @@ int compare_cells(UnitCell *a, UnitCell *b, double ltl, double atl,
 
 		if ( intmat_det(m) != +1 ) continue;
 
-		tfn = tfn_from_intmat(m);
-		nc = cell_transform(b, tfn);
+		nc = cell_transform_intmat(b, m);
 
-		if ( cells_are_similar(a, nc, ltl, atl) ) {
+		if ( compare_cell_parameters_and_orientation(a, nc, ltl, atl) ) {
 			if ( pmb != NULL ) *pmb = m;
-			tfn_free(tfn);
 			cell_free(nc);
 			return 1;
 		}
 
-		tfn_free(tfn);
 		cell_free(nc);
 
 	}
@@ -1809,3 +1857,275 @@ int compare_cells(UnitCell *a, UnitCell *b, double ltl, double atl,
 	intmat_free(m);
 	return 0;
 }
+
+
+struct cand
+{
+	Rational abc[3];
+	double fom;
+};
+
+
+static int cmpcand(const void *av, const void *bv)
+{
+	const struct cand *a = av;
+	const struct cand *b = bv;
+	return a->fom > b->fom;
+}
+
+
+static Rational *find_candidates(double len, double *a, double *b, double *c,
+                                 double ltl, int *pncand)
+{
+	Rational *r;
+	struct cand *cands;
+	const int max_cand = 1024;
+	int ncand = 0;
+	Rational *rat;
+	int nrat;
+	int nrej = 0;
+	int ia, ib, ic;
+	int i;
+
+	cands = malloc(max_cand * sizeof(struct cand));
+	if ( cands == NULL ) return NULL;
+
+	rat = rtnl_list(-5, 5, 1, 4, &nrat);
+	if ( rat == NULL ) return NULL;
+
+	for ( ia=0; ia<nrat; ia++ ) {
+	for ( ib=0; ib<nrat; ib++ ) {
+	for ( ic=0; ic<nrat; ic++ ) {
+		double vec[3];
+		double abc[3];
+		double veclen;
+		abc[0] = rtnl_as_double(rat[ia]);
+		abc[1] = rtnl_as_double(rat[ib]);
+		abc[2] = rtnl_as_double(rat[ic]);
+		vec[0] = a[0]*abc[0] + b[0]*abc[1] + c[0]*abc[2];
+		vec[1] = a[1]*abc[0] + b[1]*abc[1] + c[1]*abc[2];
+		vec[2] = a[2]*abc[0] + b[2]*abc[1] + c[2]*abc[2];
+		veclen = modulus(vec[0], vec[1], vec[2]);
+		if ( within_tolerance(len, veclen, ltl*100.0) ) {
+			if ( ncand == max_cand ) {
+				nrej++;
+			} else {
+				cands[ncand].abc[0] = rat[ia];
+				cands[ncand].abc[1] = rat[ib];
+				cands[ncand].abc[2] = rat[ic];
+				cands[ncand].fom = fabs(veclen - len);
+				ncand++;
+			}
+		}
+	}
+	}
+	}
+
+	if ( nrej ) {
+		ERROR("WARNING: Too many vector candidates (%i rejected)\n", nrej);
+	}
+
+	/* Sort by difference from reference vector length */
+	qsort(cands, ncand, sizeof(struct cand), cmpcand);
+
+	r = malloc(ncand * 3 * sizeof(Rational));
+	if ( r == 0 ) return NULL;
+
+	for ( i=0; i<ncand; i++ ) {
+		r[3*i+0] = cands[i].abc[0];
+		r[3*i+1] = cands[i].abc[1];
+		r[3*i+2] = cands[i].abc[2];
+	}
+	free(cands);
+
+	*pncand = ncand;
+	return r;
+}
+
+
+static void g6_components(double *g6, double a, double b, double c,
+                                      double al, double be, double ga)
+{
+	g6[0] = a*a;
+	g6[1] = b*b;
+	g6[2] = c*c;
+	g6[3] = 2.0*b*c*cos(al);
+	g6[4] = 2.0*a*c*cos(be);
+	g6[5] = 2.0*a*b*cos(ga);
+}
+
+
+static double g6_distance(double a1, double b1, double c1,
+                          double al1, double be1, double ga1,
+                          double a2, double b2, double c2,
+                          double al2, double be2, double ga2)
+{
+	double g1[6], g2[6];
+	int i;
+	double total = 0.0;
+	g6_components(g1, a1, b1, c1, al1, be1, ga1);
+	g6_components(g2, a2, b2, c2, al2, be2, ga2);
+	for ( i=0; i<6; i++ ) {
+		total += (g1[i]-g2[i])*(g1[i]-g2[i]);
+	}
+	return sqrt(total);
+}
+
+
+/**
+ * compare_reindexed_cell_parameters:
+ * @cell_in: A %UnitCell
+ * @reference_in: Another %UnitCell
+ * @ltl: Maximum allowable fractional difference in direct-space axis lengths
+ * @atl: Maximum allowable difference in direct-space angles (in radians)
+ * @pmb: Place to store pointer to matrix
+ *
+ * Compare the @cell_in with @reference_in.  If @cell is a derivative lattice
+ * of @reference, within fractional axis length difference @ltl and absolute angle
+ * difference @atl (in radians), this function returns non-zero and stores the
+ * transformation which needs to be applied to @cell_in at @pmb.
+ *
+ * Note that the tolerances will be applied to the primitive unit cell.  If
+ * the reference cell is centered, a primitive unit cell will first be calculated.
+ *
+ * Subject to the tolerances, this function will find the transformation which
+ * gives the best match to the reference cell, using the Euclidian norm in
+ * G6 [see e.g. Andrews and Bernstein, Acta Cryst. A44 (1988) p1009].
+ *
+ * Only the cell parameters will be compared.  The relative orientations are
+ * irrelevant.
+ *
+ * Returns: non-zero if the cells match, zero for no match or error.
+ *
+ */
+int compare_reindexed_cell_parameters(UnitCell *cell_in, UnitCell *reference_in,
+                                      double ltl, double atl,
+                                      RationalMatrix **pmb)
+{
+	UnitCell *cell;
+	UnitCell *reference;
+	IntegerMatrix *CBint;
+	RationalMatrix *CiA;
+	RationalMatrix *CB;
+	RationalMatrix *m;
+	double a, b, c, al, be, ga;
+	double av[3], bv[3], cv[3];
+	Rational *cand_a;
+	Rational *cand_b;
+	Rational *cand_c;
+	int ncand_a, ncand_b, ncand_c;
+	int ia, ib;
+	RationalMatrix *MCiA;
+	RationalMatrix *CBMCiA;
+	double min_dist = +INFINITY;
+
+	/* Actually compare against primitive version of reference */
+	reference = uncenter_cell(reference_in, &CBint, NULL);
+	if ( reference == NULL ) return 0;
+	CB = rtnl_mtx_from_intmat(CBint);
+	intmat_free(CBint);
+
+	/* Actually compare primitive version of cell */
+	cell = uncenter_cell(cell_in, NULL, &CiA);
+	if ( cell == NULL ) return 0;
+
+	/* Get target parameters */
+	cell_get_parameters(reference, &a, &b, &c, &al, &be, &ga);
+	cell_get_cartesian(cell, &av[0], &av[1], &av[2],
+	                         &bv[0], &bv[1], &bv[2],
+	                         &cv[0], &cv[1], &cv[2]);
+
+	/* Find vectors in 'cell' with lengths close to a, b and c */
+	cand_a = find_candidates(a, av, bv, cv, ltl, &ncand_a);
+	cand_b = find_candidates(b, av, bv, cv, ltl, &ncand_b);
+	cand_c = find_candidates(c, av, bv, cv, ltl, &ncand_c);
+
+	m = rtnl_mtx_new(3, 3);
+	MCiA = rtnl_mtx_new(3, 3);
+	CBMCiA = rtnl_mtx_new(3, 3);
+	for ( ia=0; ia<ncand_a; ia++ ) {
+		for ( ib=0; ib<ncand_b; ib++ ) {
+
+			UnitCell *test;
+			double at, bt, ct, alt, bet, gat;
+			double dist;
+			int ic = 0;
+
+			/* Form the matrix using the first candidate for c */
+			rtnl_mtx_set(m, 0, 0, cand_a[3*ia+0]);
+			rtnl_mtx_set(m, 1, 0, cand_a[3*ia+1]);
+			rtnl_mtx_set(m, 2, 0, cand_a[3*ia+2]);
+			rtnl_mtx_set(m, 0, 1, cand_b[3*ib+0]);
+			rtnl_mtx_set(m, 1, 1, cand_b[3*ib+1]);
+			rtnl_mtx_set(m, 2, 1, cand_b[3*ib+2]);
+			rtnl_mtx_set(m, 0, 2, cand_c[3*ic+0]);
+			rtnl_mtx_set(m, 1, 2, cand_c[3*ic+1]);
+			rtnl_mtx_set(m, 2, 2, cand_c[3*ic+2]);
+
+			/* Check angle between a and b */
+			test = cell_transform_rational(cell, m);
+			cell_get_parameters(test, &at, &bt, &ct, &alt, &bet, &gat);
+			cell_free(test);
+			if ( fabs(gat - ga) > atl ) continue;
+
+			/* Gamma OK, now look for place for c axis */
+			for ( ic=0; ic<ncand_c; ic++ ) {
+
+				rtnl_mtx_set(m, 0, 0, cand_a[3*ia+0]);
+				rtnl_mtx_set(m, 1, 0, cand_a[3*ia+1]);
+				rtnl_mtx_set(m, 2, 0, cand_a[3*ia+2]);
+				rtnl_mtx_set(m, 0, 1, cand_b[3*ib+0]);
+				rtnl_mtx_set(m, 1, 1, cand_b[3*ib+1]);
+				rtnl_mtx_set(m, 2, 1, cand_b[3*ib+2]);
+				rtnl_mtx_set(m, 0, 2, cand_c[3*ic+0]);
+				rtnl_mtx_set(m, 1, 2, cand_c[3*ic+1]);
+				rtnl_mtx_set(m, 2, 2, cand_c[3*ic+2]);
+
+				if ( rtnl_cmp(rtnl_mtx_det(m),rtnl_zero()) == 0 ) continue;
+
+				test = cell_transform_rational(cell, m);
+				cell_get_parameters(test, &at, &bt, &ct, &alt, &bet, &gat);
+				if ( !right_handed(test) ) {
+					cell_free(test);
+					continue;
+				}
+				if ( fabs(alt - al) > atl ) {
+					cell_free(test);
+					continue;
+				}
+				if ( fabs(bet - be) > atl ) {
+					cell_free(test);
+					continue;
+				}
+
+				dist = g6_distance(at, bt, ct, alt, bet, gat,
+				                   a, b, c, al, be, ga);
+				if ( dist < min_dist ) {
+					min_dist = dist;
+					rtnl_mtx_mtxmult(m, CiA, MCiA);
+				}
+
+				cell_free(test);
+
+			}
+		}
+	}
+
+	rtnl_mtx_free(m);
+	free(cand_a);
+	free(cand_b);
+	free(cand_c);
+
+	if ( isinf(min_dist) ) {
+		rtnl_mtx_free(CBMCiA);
+		rtnl_mtx_free(MCiA);
+		*pmb = NULL;
+		return 0;
+	}
+
+	/* Solution found */
+	rtnl_mtx_mtxmult(CB, MCiA, CBMCiA);
+	rtnl_mtx_free(MCiA);
+	*pmb = CBMCiA;
+	return 1;
+}
diff --git a/libcrystfel/src/cell-utils.h b/libcrystfel/src/cell-utils.h
index 5e2b2825..a97f2bfb 100644
--- a/libcrystfel/src/cell-utils.h
+++ b/libcrystfel/src/cell-utils.h
@@ -3,13 +3,13 @@
  *
  * Unit Cell utility functions
  *
- * Copyright © 2012-2017 Deutsches Elektronen-Synchrotron DESY,
+ * Copyright © 2012-2019 Deutsches Elektronen-Synchrotron DESY,
  *                       a research centre of the Helmholtz Association.
  * Copyright © 2012 Lorenzo Galli
  *
  * Authors:
- *   2009-2013,2014,2017 Thomas White <taw@physics.org>
- *   2012                 Lorenzo Galli
+ *   2009-2018 Thomas White <taw@physics.org>
+ *   2012      Lorenzo Galli
  *
  * This file is part of CrystFEL.
  *
@@ -49,9 +49,9 @@ extern double resolution(UnitCell *cell,
 extern UnitCell *cell_rotate(UnitCell *in, struct quaternion quat);
 extern UnitCell *rotate_cell(UnitCell *in, double omega, double phi,
                              double rot);
-extern UnitCell *transform_cell_gsl(UnitCell *in, gsl_matrix *m);
 
 extern void cell_print(UnitCell *cell);
+extern void cell_print_full(UnitCell *cell);
 
 extern UnitCell *match_cell(UnitCell *cell, UnitCell *tempcell, int verbose,
                             const float *ltl, int reduce);
@@ -66,7 +66,8 @@ extern int cell_is_sensible(UnitCell *cell);
 
 extern int validate_cell(UnitCell *cell);
 
-extern UnitCell *uncenter_cell(UnitCell *in, UnitCellTransformation **t);
+extern UnitCell *uncenter_cell(UnitCell *in, IntegerMatrix **pC,
+                               RationalMatrix **pCi);
 
 extern int bravais_lattice(UnitCell *cell);
 
@@ -80,8 +81,24 @@ extern int forbidden_reflection(UnitCell *cell,
 
 extern double cell_get_volume(UnitCell *cell);
 
-extern int compare_cells(UnitCell *a, UnitCell *b, double ltl, double atl,
-                         IntegerMatrix **pmb);
+extern int compare_cell_parameters(UnitCell *cell1, UnitCell *cell2,
+                                   float ltl, float atl);
+
+
+extern int compare_cell_parameters_and_orientation(UnitCell *cell1,
+                                                   UnitCell *cell2,
+                                                   const double ltl,
+                                                   const double atl);
+
+extern int compare_reindexed_cell_parameters_and_orientation(UnitCell *a,
+                                                             UnitCell *b,
+                                                             double ltl,
+                                                             double atl,
+                                                             IntegerMatrix **pmb);
+
+extern int compare_reindexed_cell_parameters(UnitCell *cell, UnitCell *reference,
+                                             double ltl, double atl,
+                                             RationalMatrix **pmb);
 
 #ifdef __cplusplus
 }
diff --git a/libcrystfel/src/cell.c b/libcrystfel/src/cell.c
index cc18b49d..ce9d37bb 100644
--- a/libcrystfel/src/cell.c
+++ b/libcrystfel/src/cell.c
@@ -45,6 +45,8 @@
 #include "cell.h"
 #include "utils.h"
 #include "image.h"
+#include "integer_matrix.h"
+#include "rational.h"
 
 
 /**
@@ -95,6 +97,20 @@ struct _unitcell {
 	char         unique_axis;
 };
 
+typedef enum {
+	CMASK_P = 1<<0,
+	CMASK_A = 1<<1,
+	CMASK_B = 1<<2,
+	CMASK_C = 1<<3,
+	CMASK_I = 1<<4,
+	CMASK_F = 1<<5,
+	CMASK_H = 1<<6,
+	CMASK_R = 1<<7
+} CenteringMask;
+
+#define CMASK_ALL (CMASK_P | CMASK_A | CMASK_B | CMASK_C | CMASK_I \
+                     | CMASK_F | CMASK_H | CMASK_R)
+
 
 /************************** Setters and Constructors **************************/
 
@@ -352,13 +368,13 @@ static int cell_invert(double ax, double ay, double az,
 		return 1;
 	}
 	gsl_matrix_set(m, 0, 0, ax);
-	gsl_matrix_set(m, 0, 1, bx);
-	gsl_matrix_set(m, 0, 2, cx);
 	gsl_matrix_set(m, 1, 0, ay);
-	gsl_matrix_set(m, 1, 1, by);
-	gsl_matrix_set(m, 1, 2, cy);
 	gsl_matrix_set(m, 2, 0, az);
+	gsl_matrix_set(m, 0, 1, bx);
+	gsl_matrix_set(m, 1, 1, by);
 	gsl_matrix_set(m, 2, 1, bz);
+	gsl_matrix_set(m, 0, 2, cx);
+	gsl_matrix_set(m, 1, 2, cy);
 	gsl_matrix_set(m, 2, 2, cz);
 
 	/* Invert */
@@ -394,13 +410,13 @@ static int cell_invert(double ax, double ay, double az,
 	gsl_matrix_transpose(inv);
 
 	*asx = gsl_matrix_get(inv, 0, 0);
-	*bsx = gsl_matrix_get(inv, 0, 1);
-	*csx = gsl_matrix_get(inv, 0, 2);
 	*asy = gsl_matrix_get(inv, 1, 0);
-	*bsy = gsl_matrix_get(inv, 1, 1);
-	*csy = gsl_matrix_get(inv, 1, 2);
 	*asz = gsl_matrix_get(inv, 2, 0);
+	*bsx = gsl_matrix_get(inv, 0, 1);
+	*bsy = gsl_matrix_get(inv, 1, 1);
 	*bsz = gsl_matrix_get(inv, 2, 1);
+	*csx = gsl_matrix_get(inv, 0, 2);
+	*csy = gsl_matrix_get(inv, 1, 2);
 	*csz = gsl_matrix_get(inv, 2, 2);
 
 	gsl_matrix_free(inv);
@@ -600,333 +616,405 @@ const char *cell_rep(UnitCell *cell)
 }
 
 
-struct _unitcelltransformation
+UnitCell *cell_transform_gsl_direct(UnitCell *in, gsl_matrix *m)
 {
-	gsl_matrix *m;
-};
-
-
-/**
- * tfn_inverse:
- * @t: A %UnitCellTransformation.
- *
- * 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: The inverse of @t.
- *
- */
-UnitCellTransformation *tfn_inverse(UnitCellTransformation *t)
-{
-	int s;
-	gsl_matrix *m;
-	gsl_matrix *inv;
-	gsl_permutation *perm;
-	UnitCellTransformation *out;
-
-	m = gsl_matrix_alloc(3, 3);
-	if ( m == NULL ) return NULL;
-
-	out = tfn_identity();
-	if ( out == NULL ) {
-		gsl_matrix_free(m);
-		return NULL;
-	}
-
-	gsl_matrix_memcpy(m, t->m);
-
-	perm = gsl_permutation_alloc(m->size1);
-	if ( perm == NULL ) {
-		ERROR("Couldn't allocate permutation\n");
-		return NULL;
-	}
-	inv = gsl_matrix_alloc(m->size1, m->size2);
-	if ( inv == NULL ) {
-		ERROR("Couldn't allocate inverse\n");
-		gsl_permutation_free(perm);
-		return NULL;
-	}
-	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 transformation matrix\n");
-		gsl_permutation_free(perm);
-		return NULL;
-	}
-	gsl_permutation_free(perm);
+	gsl_matrix *c;
+	double asx, asy, asz;
+	double bsx, bsy, bsz;
+	double csx, csy, csz;
+	gsl_matrix *res;
+	UnitCell *out;
 
-	gsl_matrix_free(out->m);
-	gsl_matrix_free(m);
-	out->m = inv;
+	cell_get_cartesian(in, &asx, &asy, &asz, &bsx, &bsy,
+	                       &bsz, &csx, &csy, &csz);
+
+	c = gsl_matrix_alloc(3, 3);
+	gsl_matrix_set(c, 0, 0, asx);
+	gsl_matrix_set(c, 1, 0, asy);
+	gsl_matrix_set(c, 2, 0, asz);
+	gsl_matrix_set(c, 0, 1, bsx);
+	gsl_matrix_set(c, 1, 1, bsy);
+	gsl_matrix_set(c, 2, 1, bsz);
+	gsl_matrix_set(c, 0, 2, csx);
+	gsl_matrix_set(c, 1, 2, csy);
+	gsl_matrix_set(c, 2, 2, csz);
+
+	res = gsl_matrix_calloc(3, 3);
+	gsl_blas_dgemm(CblasNoTrans, CblasNoTrans, 1.0, c, m, 0.0, res);
+
+	out = cell_new_from_cell(in);
+	cell_set_cartesian(out, gsl_matrix_get(res, 0, 0),
+	                        gsl_matrix_get(res, 1, 0),
+	                        gsl_matrix_get(res, 2, 0),
+	                        gsl_matrix_get(res, 0, 1),
+	                        gsl_matrix_get(res, 1, 1),
+	                        gsl_matrix_get(res, 2, 1),
+	                        gsl_matrix_get(res, 0, 2),
+	                        gsl_matrix_get(res, 1, 2),
+	                        gsl_matrix_get(res, 2, 2));
+
+	gsl_matrix_free(res);
+	gsl_matrix_free(c);
 	return out;
 }
 
 
-/**
- * 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;
-	gsl_matrix *m;
-	gsl_matrix *a;
+static int centering_has_point(char cen, Rational *p)
+{
+	/* First, put the point into the range 0..1 */
+	while ( rtnl_cmp(p[0], rtnl_zero()) < 0 ) p[0] = rtnl_add(p[0], rtnl(1, 1));
+	while ( rtnl_cmp(p[1], rtnl_zero()) < 0 ) p[1] = rtnl_add(p[1], rtnl(1, 1));
+	while ( rtnl_cmp(p[2], rtnl_zero()) < 0 ) p[2] = rtnl_add(p[2], rtnl(1, 1));
+	while ( rtnl_cmp(p[0], rtnl(1, 1)) >= 0 ) p[0] = rtnl_sub(p[0], rtnl(1, 1));
+	while ( rtnl_cmp(p[1], rtnl(1, 1)) >= 0 ) p[1] = rtnl_sub(p[1], rtnl(1, 1));
+	while ( rtnl_cmp(p[2], rtnl(1, 1)) >= 0 ) p[2] = rtnl_sub(p[2], rtnl(1, 1));
+
+	/* 0,0,0 is present in all centerings */
+	if ( (rtnl_cmp(p[0], rtnl_zero()) == 0)
+	  && (rtnl_cmp(p[1], rtnl_zero()) == 0)
+	  && (rtnl_cmp(p[2], rtnl_zero()) == 0) ) return 1;
+
+	/* Only I has 1/2 , 1/2, 1/2 */
+	if ( (rtnl_cmp(p[0], rtnl(1,2)) == 0)
+	  && (rtnl_cmp(p[1], rtnl(1,2)) == 0)
+	  && (rtnl_cmp(p[2], rtnl(1,2)) == 0)
+	  && (cen == 'I') ) return 1;
+
+	/* A or F has 0 , 1/2, 1/2 */
+	if ( (rtnl_cmp(p[0], rtnl_zero()) == 0)
+	  && (rtnl_cmp(p[1], rtnl(1,2)) == 0)
+	  && (rtnl_cmp(p[2], rtnl(1,2)) == 0)
+	  && ((cen == 'A') || (cen == 'F')) ) return 1;
+
+	/* B or F has 1/2 , 0 , 1/2 */
+	if ( (rtnl_cmp(p[0], rtnl(1,2)) == 0)
+	  && (rtnl_cmp(p[1], rtnl_zero()) == 0)
+	  && (rtnl_cmp(p[2], rtnl(1,2)) == 0)
+	  && ((cen == 'B') || (cen == 'F')) ) return 1;
+
+	/* C or F has 1/2 , 1/2 , 0 */
+	if ( (rtnl_cmp(p[0], rtnl(1,2)) == 0)
+	  && (rtnl_cmp(p[1], rtnl(1,2)) == 0)
+	  && (rtnl_cmp(p[2], rtnl_zero()) == 0)
+	  && ((cen == 'C') || (cen == 'F')) ) return 1;
+
+	/* H has 2/3 , 1/3 , 1/3 */
+	if ( (rtnl_cmp(p[0], rtnl(2,3)) == 0)
+	  && (rtnl_cmp(p[1], rtnl(1,3)) == 0)
+	  && (rtnl_cmp(p[2], rtnl(1,3)) == 0)
+	  && (cen == 'H') ) return 1;
+
+	/* H has 1/3 , 2/3 , 2/3 */
+	if ( (rtnl_cmp(p[0], rtnl(1,3)) == 0)
+	  && (rtnl_cmp(p[1], rtnl(2,3)) == 0)
+	  && (rtnl_cmp(p[2], rtnl(2,3)) == 0)
+	  && (cen == 'H') ) return 1;
 
-	if ( t == NULL ) return NULL;
+	return 0;
+}
 
-	out = cell_new_from_cell(cell);
-	if ( out == NULL ) return NULL;
 
-	if ( cell_get_cartesian(out, &ax, &ay, &az,
-	                             &bx, &by, &bz,
-	                             &cx, &cy, &cz) ) return NULL;
+static void maybe_eliminate(CenteringMask c, CenteringMask *cmask, Rational *nc,
+                            char cen)
+{
+	/* Skip test if this centering isn't even a candidate */
+	if ( !(*cmask & c) ) return;
 
-	m = gsl_matrix_alloc(3,3);
-	a = gsl_matrix_calloc(3,3);
-	if ( (m == NULL) || (a == NULL) ) {
-		cell_free(out);
-		return NULL;
+	if ( !centering_has_point(cen, nc) ) {
+		*cmask |= c;
+		*cmask ^= c;
 	}
+}
 
-	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);
 
-	gsl_blas_dgemm(CblasNoTrans, CblasNoTrans, 1.0, t->m, m, 0.0, a);
+/* Check if the point x,y,z in the original cell matches any lattice point
+ * in the transformed cell */
+static void check_point_fwd(RationalMatrix *P, CenteringMask *cmask,
+                            Rational x, Rational y, Rational z)
+{
+	Rational c[3] = {x, y, z};
+	Rational nc[3];
 
-	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));
+	/* Transform the lattice point */
+	transform_fractional_coords_rtnl(P, c, nc);
 
-	gsl_matrix_free(a);
-	gsl_matrix_free(m);
-
-	return out;
+	/* Eliminate any centerings which don't include the transformed point */
+	maybe_eliminate(CMASK_P, cmask, nc, 'P');
+	maybe_eliminate(CMASK_R, cmask, nc, 'R');
+	maybe_eliminate(CMASK_A, cmask, nc, 'A');
+	maybe_eliminate(CMASK_B, cmask, nc, 'B');
+	maybe_eliminate(CMASK_C, cmask, nc, 'C');
+	maybe_eliminate(CMASK_I, cmask, nc, 'I');
+	maybe_eliminate(CMASK_F, cmask, nc, 'F');
+	maybe_eliminate(CMASK_H, cmask, nc, 'H');
 }
 
 
-/**
- * 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)
+/* Check if the point x,y,z in the transformed cell matches any lattice point
+ * in the original cell.  If not, eliminate "exclude" from "*mask". */
+static void check_point_bwd(RationalMatrix *P, CenteringMask *mask,
+                            char cen, CenteringMask exclude,
+                            Rational x, Rational y, Rational z)
 {
-	UnitCellTransformation *inv;
-	UnitCell *out;
+	Rational nc[3];
+	Rational c[3] = {x, y, z};
 
-	inv = tfn_inverse(t);
-	out = cell_transform(cell, inv);
-	tfn_free(inv);
-	return out;
+	transform_fractional_coords_rtnl_inverse(P, c, nc);
+
+	if ( !centering_has_point(cen, nc) ) {
+		*mask |= exclude;
+		*mask ^= exclude;  /* Unset bits */
+	}
 }
 
 
-/**
- * tfn_identity:
- *
- * Returns: A %UnitCellTransformation corresponding to an identity operation.
- *
- */
-UnitCellTransformation *tfn_identity()
+static char cmask_decode(CenteringMask mask)
 {
-	UnitCellTransformation *tfn;
+	char res[32];
 
-	tfn = calloc(1, sizeof(UnitCellTransformation));
-	if ( tfn == NULL ) return NULL;
+	res[0] = '\0';
 
-	tfn->m = gsl_matrix_alloc(3, 3);
-	if ( tfn->m == NULL ) {
-		free(tfn);
-		return NULL;
-	}
-
-	gsl_matrix_set_identity(tfn->m);
+	if ( mask & CMASK_H ) strcat(res, "H");
+	if ( mask & CMASK_F ) strcat(res, "F");
+	if ( mask & CMASK_I ) strcat(res, "I");
+	if ( mask & CMASK_A ) strcat(res, "A");
+	if ( mask & CMASK_B ) strcat(res, "B");
+	if ( mask & CMASK_C ) strcat(res, "C");
+	if ( mask & CMASK_P ) strcat(res, "P");
+	if ( mask & CMASK_R ) strcat(res, "R");
 
-	return tfn;
+	if ( strlen(res) == 0 ) return '?';
+	return res[0];
 }
 
 
+static char determine_centering(RationalMatrix *P, char cen)
+{
+	CenteringMask cmask = CMASK_ALL;
+
+	/* Check whether the current centering can provide all the lattice
+	 * points for the transformed cell.  Eliminate any centerings for which
+	 * it can't. */
+	check_point_bwd(P, &cmask, cen, CMASK_A | CMASK_F, rtnl_zero(), rtnl(1,2), rtnl(1,2));
+	check_point_bwd(P, &cmask, cen, CMASK_B | CMASK_F, rtnl(1,2), rtnl_zero(), rtnl(1,2));
+	check_point_bwd(P, &cmask, cen, CMASK_C | CMASK_F, rtnl(1,2), rtnl(1,2), rtnl_zero());
+	check_point_bwd(P, &cmask, cen, CMASK_I, rtnl(1,2), rtnl(1,2), rtnl(1,2));
+	check_point_bwd(P, &cmask, cen, CMASK_H, rtnl(2,3), rtnl(1,3), rtnl(1,3));
+	check_point_bwd(P, &cmask, cen, CMASK_H, rtnl(1,3), rtnl(2,3), rtnl(2,3));
+	check_point_bwd(P, &cmask, cen, CMASK_ALL, rtnl(1,1), rtnl(1,1), rtnl(1,1));
+
+	/* Check whether the current centering's lattice points will all
+	 * coincide with lattice points in the new centering.  Eliminate any
+	 * centerings for which they don't (they give "excess lattice points"). */
+	switch ( cen ) {
+
+		case 'P' :
+		case 'R' :
+		check_point_fwd(P, &cmask, rtnl(1,1), rtnl(1,1), rtnl(1,1));
+		break;
+
+		case 'A' :
+		check_point_fwd(P, &cmask, rtnl(1,1), rtnl(1,1), rtnl(1,1));
+		check_point_fwd(P, &cmask, rtnl_zero(), rtnl(1,2), rtnl(1,2));
+		break;
+
+		case 'B' :
+		check_point_fwd(P, &cmask, rtnl(1,1), rtnl(1,1), rtnl(1,1));
+		check_point_fwd(P, &cmask, rtnl(1,2), rtnl_zero(), rtnl(1,2));
+		break;
+
+		case 'C' :
+		check_point_fwd(P, &cmask, rtnl(1,1), rtnl(1,1), rtnl(1,1));
+		check_point_fwd(P, &cmask, rtnl(1,2), rtnl(1,2), rtnl_zero());
+		break;
+
+		case 'I' :
+		check_point_fwd(P, &cmask, rtnl(1,1), rtnl(1,1), rtnl(1,1));
+		check_point_fwd(P, &cmask, rtnl(1,2), rtnl(1,2), rtnl(1,2));
+		break;
+
+		case 'F' :
+		check_point_fwd(P, &cmask, rtnl(1,1), rtnl(1,1), rtnl(1,1));
+		check_point_fwd(P, &cmask, rtnl_zero(), rtnl(1,2), rtnl(1,2));
+		check_point_fwd(P, &cmask, rtnl(1,2), rtnl_zero(), rtnl(1,2));
+		check_point_fwd(P, &cmask, rtnl(1,2), rtnl(1,2), rtnl_zero());
+		break;
+
+		case 'H' :
+		check_point_fwd(P, &cmask, rtnl(1,1), rtnl(1,1), rtnl(1,1));
+		check_point_fwd(P, &cmask, rtnl(2,3), rtnl(1,3), rtnl(1,3));
+		check_point_fwd(P, &cmask, rtnl(1,3), rtnl(2,3), rtnl(2,3));
+		break;
 
-/**
- * tfn_from_intmat:
- * @m: An %IntegerMatrix
- *
- * Returns: A %UnitCellTransformation corresponding to @m.
- *
- */
-UnitCellTransformation *tfn_from_intmat(IntegerMatrix *m)
-{
-	UnitCellTransformation *tfn;
-	int i, j;
-
-	tfn = tfn_identity();
-	if ( tfn == NULL ) return NULL;
-
-	for ( i=0; i<3; i++ ) {
-	for ( j=0; j<3; j++ ) {
-		gsl_matrix_set(tfn->m, i, j, intmat_get(m, i, j));
-	}
 	}
 
-	return tfn;
+	return cmask_decode(cmask);
 }
 
 
 /**
- * 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
+ * cell_transform_rational:
+ * @cell: A %UnitCell.
+ * @t: A %RationalMatrix.
+ *
+ * Applies @t to @cell.
  *
- * 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.
+ * This function will determine the centering of the resulting unit cell,
+ * producing '?' if any lattice points cannot be accounted for.  Note that if
+ * there are 'excess' lattice points in the transformed cell, the centering
+ * will still be '?' even if the lattice points for another centering are
+ * all present.
+ *
+ * The lattice type will be set to triclinic, and the unique axis to '?'.
+ *
+ * Returns: Transformed copy of @cell.
  *
  */
-void tfn_combine(UnitCellTransformation *t, double *na, double *nb, double *nc)
+UnitCell *cell_transform_rational(UnitCell *cell, RationalMatrix *m)
 {
-	gsl_matrix *a;
-	gsl_matrix *n;
+	UnitCell *out;
+	gsl_matrix *tm;
+	char ncen;
+	int i, j;
+	Rational det;
 
-	n = gsl_matrix_alloc(3, 3);
-	a = gsl_matrix_calloc(3, 3);
-	if ( (n == NULL) || (a == NULL) ) {
-		return;
+	if ( m == NULL ) return NULL;
+
+	det = rtnl_mtx_det(m);
+	if ( rtnl_cmp(det, rtnl_zero()) == 0 ) return NULL;
+
+	tm = gsl_matrix_alloc(3,3);
+	if ( tm == NULL ) {
+		return NULL;
 	}
-	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]);
 
-	free(na);
-	free(nb);
-	free(nc);
+	for ( i=0; i<3; i++ ) {
+		for ( j=0; j<3; j++ ) {
+			gsl_matrix_set(tm, i, j,
+			               rtnl_as_double(rtnl_mtx_get(m, i, j)));
+		}
+	}
+
+	out = cell_transform_gsl_direct(cell, tm);
+	gsl_matrix_free(tm);
+
+	ncen = determine_centering(m, cell_get_centering(cell));
+	cell_set_centering(out, ncen);
 
-	gsl_blas_dgemm(CblasNoTrans, CblasNoTrans, 1.0, n, t->m, 0.0, a);
+	/* FIXME: Update unique axis, lattice type */
+	cell_set_lattice_type(out, L_TRICLINIC);
+	cell_set_unique_axis(out, '?');
 
-	gsl_matrix_free(t->m);
-	t->m = a;
-	gsl_matrix_free(n);
+	return out;
 }
 
 
 /**
- * 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
+ * cell_transform_intmat:
+ * @cell: A %UnitCell.
+ * @t: An %IntegerMatrix.
+ *
+ * Applies @t to @cell.
  *
- * 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));
+ * This is just a convenience function which turns @m into a %RationalMatrix
+ * and then calls cell_transform_rational().  See the documentation for that
+ * function for some important information.
+ *
+ * Returns: Transformed copy of @cell.
  *
  */
-double *tfn_vector(double a, double b, double c)
+UnitCell *cell_transform_intmat(UnitCell *cell, IntegerMatrix *m)
 {
-	double *vec = malloc(3*sizeof(double));
-	if ( vec == NULL ) return NULL;
-	vec[0] = a;  vec[1] = b;  vec[2] = c;
-	return vec;
+	UnitCell *ans;
+	RationalMatrix *mtx = rtnl_mtx_from_intmat(m);
+	ans = cell_transform_rational(cell, mtx);
+	rtnl_mtx_free(mtx);
+	return ans;
 }
 
 
 /**
- * tfn_print:
- * @t: A %UnitCellTransformation
+ * cell_transform_rational_inverse:
+ * @cell: A %UnitCell.
+ * @m: A %RationalMatrix
  *
- * Prints information about @t to stderr.
+ * Applies the inverse of @m to @cell.
+ *
+ * Returns: Transformed copy of @cell.
  *
  */
-void tfn_print(UnitCellTransformation *t)
+UnitCell *cell_transform_rational_inverse(UnitCell *cell, RationalMatrix *m)
 {
+	UnitCell *out;
+	gsl_matrix *tm;
+	gsl_matrix *inv;
 	gsl_permutation *perm;
-	gsl_matrix *lu;
 	int s;
+	int i, j;
+
+	if ( m == NULL ) return NULL;
 
-	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;
+	tm = gsl_matrix_alloc(3,3);
+	if ( tm == NULL ) {
+		return NULL;
 	}
 
-	gsl_matrix_memcpy(lu, t->m);
+	for ( i=0; i<3; i++ ) {
+		for ( j=0; j<3; j++ ) {
+			gsl_matrix_set(tm, i, j,
+			               rtnl_as_double(rtnl_mtx_get(m, i, j)));
+		}
+	}
 
-	perm = gsl_permutation_alloc(t->m->size1);
+	perm = gsl_permutation_alloc(3);
 	if ( perm == NULL ) {
-		ERROR("Couldn't allocate permutation.\n");
-		gsl_matrix_free(lu);
-		return;
+		ERROR("Couldn't allocate permutation\n");
+		return NULL;
 	}
-	if ( gsl_linalg_LU_decomp(lu, perm, &s) ) {
-		ERROR("LU decomposition failed.\n");
+	inv = gsl_matrix_alloc(3, 3);
+	if ( inv == NULL ) {
+		ERROR("Couldn't allocate inverse\n");
+		gsl_permutation_free(perm);
+		return NULL;
+	}
+	if ( gsl_linalg_LU_decomp(tm, perm, &s) ) {
+		ERROR("Couldn't decompose matrix\n");
 		gsl_permutation_free(perm);
-		gsl_matrix_free(lu);
-		return;
+		return NULL;
+	}
+	if ( gsl_linalg_LU_invert(tm, perm, inv)  ) {
+		ERROR("Couldn't invert transformation matrix\n");
+		gsl_permutation_free(perm);
+		return NULL;
 	}
+	gsl_permutation_free(perm);
+
+	out = cell_transform_gsl_direct(cell, inv);
+
+	/* FIXME: Update centering, unique axis, lattice type */
+
+	gsl_matrix_free(tm);
+	gsl_matrix_free(inv);
 
-	STATUS("Transformation determinant = %.2f\n", gsl_linalg_LU_det(lu, s));
+	return out;
 }
 
 
 /**
- * tfn_free:
- * @t: A %UnitCellTransformation
+ * cell_transform_intmat_inverse:
+ * @cell: A %UnitCell.
+ * @m: An %IntegerMatrix
  *
- * Frees all resources associated with @t.
+ * Applies the inverse of @m to @cell.
+ *
+ * Returns: Transformed copy of @cell.
  *
  */
-void tfn_free(UnitCellTransformation *t)
+UnitCell *cell_transform_intmat_inverse(UnitCell *cell, IntegerMatrix *m)
 {
-	gsl_matrix_free(t->m);
-	free(t);
+	UnitCell *ans;
+	RationalMatrix *mtx = rtnl_mtx_from_intmat(m);
+	ans = cell_transform_rational_inverse(cell, mtx);
+	rtnl_mtx_free(mtx);
+	return ans;
 }
diff --git a/libcrystfel/src/cell.h b/libcrystfel/src/cell.h
index 39e6a1ed..c868cd29 100644
--- a/libcrystfel/src/cell.h
+++ b/libcrystfel/src/cell.h
@@ -91,14 +91,6 @@ typedef enum
 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;
-
 #ifdef __cplusplus
 extern "C" {
 #endif
@@ -156,18 +148,13 @@ extern void cell_set_unique_axis(UnitCell *cell, char unique_axis);
 
 extern const char *cell_rep(UnitCell *cell);
 
-extern UnitCell *cell_transform(UnitCell *cell, UnitCellTransformation *t);
-extern UnitCell *cell_transform_inverse(UnitCell *cell,
-                                        UnitCellTransformation *t);
-
-extern UnitCellTransformation *tfn_identity(void);
-extern UnitCellTransformation *tfn_from_intmat(IntegerMatrix *m);
-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);
+extern UnitCell *cell_transform_gsl_direct(UnitCell *in, gsl_matrix *m);
+
+extern UnitCell *cell_transform_rational(UnitCell *cell, RationalMatrix *m);
+extern UnitCell *cell_transform_rational_inverse(UnitCell *cell, RationalMatrix *m);
+
+extern UnitCell *cell_transform_intmat(UnitCell *cell, IntegerMatrix *m);
+extern UnitCell *cell_transform_intmat_inverse(UnitCell *cell, IntegerMatrix *m);
 
 #ifdef __cplusplus
 }
diff --git a/libcrystfel/src/index.c b/libcrystfel/src/index.c
index f72334c4..ad9c0de3 100644
--- a/libcrystfel/src/index.c
+++ b/libcrystfel/src/index.c
@@ -698,9 +698,9 @@ static int try_indexer(struct image *image, IndexingMethod indm,
 			/* Don't do similarity check against bad crystals */
 			if ( crystal_get_user_flag(that_cr) ) continue;
 
-			if ( compare_cells(crystal_get_cell(cr),
-			                   crystal_get_cell(that_cr),
-			                   0.1, deg2rad(0.5), NULL) )
+			if ( compare_reindexed_cell_parameters_and_orientation(crystal_get_cell(cr),
+			                                                       crystal_get_cell(that_cr),
+			                                                       0.1, deg2rad(0.5), NULL) )
 			{
 				crystal_set_user_flag(cr, 1);
 			}
diff --git a/libcrystfel/src/integer_matrix.c b/libcrystfel/src/integer_matrix.c
index c31072b4..90d58b35 100644
--- a/libcrystfel/src/integer_matrix.c
+++ b/libcrystfel/src/integer_matrix.c
@@ -35,7 +35,9 @@
 #include <string.h>
 #include <assert.h>
 
+#include "rational.h"
 #include "integer_matrix.h"
+#include "utils.h"
 
 
 /**
@@ -95,7 +97,7 @@ IntegerMatrix *intmat_new(unsigned int rows, unsigned int cols)
  *
  * Returns: a newly allocated copy of @m, or NULL on error/
  **/
-IntegerMatrix *intmat_copy(IntegerMatrix *m)
+IntegerMatrix *intmat_copy(const IntegerMatrix *m)
 {
 	IntegerMatrix *p;
 	int i, j;
@@ -127,6 +129,19 @@ void intmat_free(IntegerMatrix *m)
 }
 
 
+void intmat_size(const IntegerMatrix *m, unsigned int *rows, unsigned int *cols)
+{
+	if ( m == NULL ) {
+		*rows = 0;
+		*cols = 0;
+		return;
+	}
+
+	*rows = m->rows;
+	*cols = m->cols;
+}
+
+
 /**
  * intmat_set:
  * @m: An %IntegerMatrix
@@ -163,32 +178,38 @@ signed int intmat_get(const IntegerMatrix *m, unsigned int i, unsigned int j)
 
 
 /**
- * intmat_intvec_mult:
- * @m: An %IntegerMatrix
- * @vec: An array of signed integers
+ * transform_indices:
+ * @P: An %IntegerMatrix
+ * @hkl: An array of signed integers
+ *
+ * Apply transformation matrix P to a set of reciprocal space Miller indices.
  *
- * Multiplies the matrix @m by the vector @vec.  The size of @vec must equal the
- * number of columns in @m, and the size of the result equals the number of rows
- * in @m.
+ * In other words:
+ * Multiplies the matrix @P by the row vector @hkl.  The size of @vec must equal
+ * the number of columns in @P, and the size of the result equals the number of
+ * rows in @P.
+ *
+ * The multiplication looks like this:
+ *    (a1, a2, a3) = (hkl1, hkl2, hkl3) P
+ * Therefore matching the notation in ITA chapter 5.1.
  *
  * Returns: a newly allocated array of signed integers containing the answer,
  * or NULL on error.
  **/
-signed int *intmat_intvec_mult(const IntegerMatrix *m, const signed int *vec)
+signed int *transform_indices(const IntegerMatrix *P, const signed int *hkl)
 {
 	signed int *ans;
-	unsigned int i;
+	unsigned int j;
 
-	ans = malloc(m->rows * sizeof(signed int));
+	ans = malloc(P->rows * sizeof(signed int));
 	if ( ans == NULL ) return NULL;
 
-	for ( i=0; i<m->rows; i++ ) {
-
-		unsigned int j;
+	for ( j=0; j<P->cols; j++ ) {
 
-		ans[i] = 0;
-		for ( j=0; j<m->cols; j++ ) {
-			ans[i] += intmat_get(m, i, j) * vec[j];
+		unsigned int i;
+		ans[j] = 0;
+		for ( i=0; i<P->rows; i++ ) {
+			ans[i] += intmat_get(P, i, j) * hkl[j];
 		}
 
 	}
@@ -342,6 +363,9 @@ static IntegerMatrix *intmat_cofactors(const IntegerMatrix *m)
  *
  * Calculates the inverse of @m.  Inefficiently.
  *
+ * Works only if the inverse of the matrix is also an integer matrix,
+ * i.e. if the determinant of @m is +/- 1.
+ *
  * Returns: the inverse of @m, or NULL on error.
  **/
 IntegerMatrix *intmat_inverse(const IntegerMatrix *m)
@@ -532,3 +556,39 @@ IntegerMatrix *intmat_identity(int size)
 
 	return m;
 }
+
+
+/**
+ * intmat_set_all_3x3
+ * @m: An %IntegerMatrix
+ *
+ * Returns: an identity %IntegerMatrix with side length @size, or NULL on error.
+ *
+ */
+void intmat_set_all_3x3(IntegerMatrix *m,
+                        signed int m11, signed int m12, signed int m13,
+                        signed int m21, signed int m22, signed int m23,
+                        signed int m31, signed int m32, signed int m33)
+{
+	intmat_set(m, 0, 0, m11);
+	intmat_set(m, 0, 1, m12);
+	intmat_set(m, 0, 2, m13);
+
+	intmat_set(m, 1, 0, m21);
+	intmat_set(m, 1, 1, m22);
+	intmat_set(m, 1, 2, m23);
+
+	intmat_set(m, 2, 0, m31);
+	intmat_set(m, 2, 1, m32);
+	intmat_set(m, 2, 2, m33);
+}
+
+
+IntegerMatrix *intmat_create_3x3(signed int m11, signed int m12, signed int m13,
+                                 signed int m21, signed int m22, signed int m23,
+                                 signed int m31, signed int m32, signed int m33)
+{
+	IntegerMatrix *t = intmat_new(3, 3);
+	intmat_set_all_3x3(t, m11, m12, m13, m21, m22, m23, m31, m32, m33);
+	return t;
+}
diff --git a/libcrystfel/src/integer_matrix.h b/libcrystfel/src/integer_matrix.h
index 6616b5e8..6fb3e399 100644
--- a/libcrystfel/src/integer_matrix.h
+++ b/libcrystfel/src/integer_matrix.h
@@ -40,25 +40,38 @@
  **/
 typedef struct _integermatrix IntegerMatrix;
 
+#include "rational.h"
+
 #ifdef __cplusplus
 extern "C" {
 #endif
 
 /* Alloc/dealloc */
 extern IntegerMatrix *intmat_new(unsigned int rows, unsigned int cols);
-extern IntegerMatrix *intmat_copy(IntegerMatrix *m);
+extern IntegerMatrix *intmat_copy(const IntegerMatrix *m);
 extern IntegerMatrix *intmat_identity(int size);
 extern void intmat_free(IntegerMatrix *m);
 
 /* Get/set */
+extern void intmat_size(const IntegerMatrix *m, unsigned int *rows,
+                        unsigned int *cols);
+
 extern void intmat_set(IntegerMatrix *m, unsigned int i, unsigned int j,
                        signed int v);
 extern signed int intmat_get(const IntegerMatrix *m,
                              unsigned int i, unsigned int j);
 
-/* Matrix-(int)vector multiplication */
-extern signed int *intmat_intvec_mult(const IntegerMatrix *m,
-                                      const signed int *vec);
+extern void intmat_set_all_3x3(IntegerMatrix *m,
+                               signed int m11, signed int m12, signed int m13,
+                               signed int m21, signed int m22, signed int m23,
+                               signed int m31, signed int m32, signed int m33);
+
+extern IntegerMatrix *intmat_create_3x3(signed int m11, signed int m12, signed int m13,
+                                        signed int m21, signed int m22, signed int m23,
+                                        signed int m31, signed int m32, signed int m33);
+
+/* Matrix-vector multiplication */
+extern signed int *transform_indices(const IntegerMatrix *P, const signed int *hkl);
 
 /* Matrix-matrix multiplication */
 extern IntegerMatrix *intmat_intmat_mult(const IntegerMatrix *a,
diff --git a/libcrystfel/src/rational.c b/libcrystfel/src/rational.c
new file mode 100644
index 00000000..65b209ff
--- /dev/null
+++ b/libcrystfel/src/rational.c
@@ -0,0 +1,675 @@
+/*
+ * rational.c
+ *
+ * A small rational number library
+ *
+ * Copyright © 2019 Deutsches Elektronen-Synchrotron DESY,
+ *                  a research centre of the Helmholtz Association.
+ *
+ * Authors:
+ *   2019 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 <string.h>
+#include <assert.h>
+#include <locale.h>
+
+#include "rational.h"
+#include "integer_matrix.h"
+#include "utils.h"
+
+
+/**
+ * SECTION:rational
+ * @short_description: Rational numbers
+ * @title: Rational numbers
+ * @section_id:
+ * @see_also:
+ * @include: "rational.h"
+ * @Image:
+ *
+ * A rational number library
+ */
+
+
+/* Eucliden algorithm for finding greatest common divisor */
+static signed int gcd(signed int a, signed int b)
+{
+	while ( b != 0 ) {
+		signed int t = b;
+		b = a % b;
+		a = t;
+	}
+	return a;
+}
+
+
+static void squish(Rational *rt)
+{
+	signed int g;
+
+	if ( rt->num == 0 ) {
+		rt->den = 1;
+		return;
+	}
+
+	g = gcd(rt->num, rt->den);
+	assert(g != 0);
+	rt->num /= g;
+	rt->den /= g;
+
+	if ( rt->den < 0 ) {
+		rt->num = -rt->num;
+		rt->den = -rt->den;
+	}
+}
+
+
+Rational rtnl_zero()
+{
+	Rational r;
+	r.num = 0;
+	r.den = 1;
+	return r;
+}
+
+
+Rational rtnl(signed long long int num, signed long long int den)
+{
+	Rational r;
+	r.num = num;
+	r.den = den;
+	squish(&r);
+	return r;
+}
+
+
+double rtnl_as_double(Rational r)
+{
+	return (double)r.num/r.den;
+}
+
+
+static void overflow(long long int c, long long int a, long long int b)
+{
+	setlocale(LC_ALL, "");
+	ERROR("Overflow detected in rational number library.\n");
+	ERROR("%'lli < %'lli * %'lli\n", c, a, b);
+	abort();
+}
+
+
+static void check_overflow(long long int c, long long int a, long long int b)
+{
+	if ( (a==0) || (b==0) ) {
+		if ( c != 0 ) overflow(c,a,b);
+	} else if ( (llabs(c) < llabs(a)) || (llabs(c) < llabs(b)) ) {
+		overflow(c,a,b);
+	}
+}
+
+
+Rational rtnl_mul(Rational a, Rational b)
+{
+	Rational r;
+	r.num = a.num * b.num;
+	r.den = a.den * b.den;
+	check_overflow(r.num, a.num, b.num);
+	check_overflow(r.den, a.den, b.den);
+	squish(&r);
+	return r;
+}
+
+
+Rational rtnl_div(Rational a, Rational b)
+{
+	signed int t = b.num;
+	b.num = b.den;
+	b.den = t;
+	return rtnl_mul(a, b);
+}
+
+
+Rational rtnl_add(Rational a, Rational b)
+{
+	Rational r, trt1, trt2;
+
+	trt1.num = a.num * b.den;
+	trt2.num = b.num * a.den;
+	check_overflow(trt1.num, a.num, b.den);
+	check_overflow(trt2.num, b.num, a.den);
+
+	trt1.den = a.den * b.den;
+	trt2.den = trt1.den;
+	check_overflow(trt1.den, a.den, b.den);
+
+	r.num = trt1.num + trt2.num;
+	r.den = trt1.den;
+	squish(&r);
+	return r;
+}
+
+
+Rational rtnl_sub(Rational a, Rational b)
+{
+	b.num = -b.num;
+	return rtnl_add(a, b);
+}
+
+
+/* -1, 0 +1 respectively for a<b, a==b, a>b */
+signed int rtnl_cmp(Rational a, Rational b)
+{
+	Rational trt1, trt2;
+
+	trt1.num = a.num * b.den;
+	trt2.num = b.num * a.den;
+
+	trt1.den = a.den * b.den;
+	trt2.den = a.den * b.den;
+
+	if ( trt1.num > trt2.num ) return +1;
+	if ( trt1.num < trt2.num ) return -1;
+	return 0;
+}
+
+
+Rational rtnl_abs(Rational a)
+{
+	Rational r = a;
+	squish(&r);
+	if ( r.num < 0 ) r.num = -r.num;
+	return r;
+}
+
+
+/**
+ * rtnl_format
+ * @rt: A %Rational
+ *
+ * Formats @rt as a string
+ *
+ * Returns: a string which should be freed by the caller
+ */
+char *rtnl_format(Rational rt)
+{
+	char *v = malloc(32);
+	if ( v == NULL ) return NULL;
+	if ( rt.den == 1 ) {
+		snprintf(v, 31, "%lli", rt.num);
+	} else {
+		snprintf(v, 31, "%lli/%lli", rt.num, rt.den);
+	}
+	return v;
+}
+
+
+Rational *rtnl_list(signed int num_min, signed int num_max,
+                    signed int den_min, signed int den_max,
+                    int *pn)
+{
+	signed int num, den;
+	Rational *list;
+	int n = 0;
+
+	list = malloc((1+num_max-num_min)*(1+den_max-den_min)*sizeof(Rational));
+	if ( list == NULL ) return NULL;
+
+	for ( num=num_min; num<=num_max; num++ ) {
+		for ( den=den_min; den<=den_max; den++ ) {
+
+			Rational r = rtnl(num, den);
+
+			/* Denominator zero? */
+			if ( den == 0 ) continue;
+
+			/* Same as last entry? */
+			if ( (n>0) && (rtnl_cmp(list[n-1], r)==0) ) continue;
+
+			/* Can be reduced? */
+			if ( gcd(num, den) != 1 ) continue;
+
+			list[n++] = r;
+		}
+	}
+	*pn = n;
+	return list;
+}
+
+
+/**
+ * SECTION:rational_matrix
+ * @short_description: Rational matrices
+ * @title: Rational matrices
+ * @section_id:
+ * @see_also:
+ * @include: "rational.h"
+ * @Image:
+ *
+ * A rational matrix library
+ */
+
+
+struct _rationalmatrix
+{
+	unsigned int rows;
+	unsigned int cols;
+	Rational *v;
+};
+
+
+/**
+ * rtnl_mtx_new:
+ * @rows: Number of rows that the new matrix is to have
+ * @cols: Number of columns that the new matrix is to have
+ *
+ * Allocates a new %RationalMatrix with all elements set to zero.
+ *
+ * Returns: a new %RationalMatrix, or NULL on error.
+ **/
+RationalMatrix *rtnl_mtx_new(unsigned int rows, unsigned int cols)
+{
+	RationalMatrix *m;
+	int i;
+
+	m = malloc(sizeof(RationalMatrix));
+	if ( m == NULL ) return NULL;
+
+	m->v = calloc(rows*cols, sizeof(Rational));
+	if ( m->v == NULL ) {
+		free(m);
+		return NULL;
+	}
+
+	m->rows = rows;
+	m->cols = cols;
+
+	for ( i=0; i<m->rows*m->cols; i++ ) {
+		m->v[i] = rtnl_zero();
+	}
+
+	return m;
+}
+
+
+RationalMatrix *rtnl_mtx_identity(int rows)
+{
+	int i;
+	RationalMatrix *m = rtnl_mtx_new(rows, rows);
+	for ( i=0; i<rows; i++ ) {
+		rtnl_mtx_set(m, i, i, rtnl(1,1));
+	}
+	return m;
+}
+
+
+RationalMatrix *rtnl_mtx_copy(const RationalMatrix *m)
+{
+	RationalMatrix *n;
+	int i;
+
+	n = rtnl_mtx_new(m->rows, m->cols);
+	if ( n == NULL ) return NULL;
+
+	for ( i=0; i<m->rows*m->cols; i++ ) {
+		n->v[i] = m->v[i];
+	}
+
+	return n;
+}
+
+
+Rational rtnl_mtx_get(const RationalMatrix *m, int i, int j)
+{
+	assert(m != NULL);
+	return m->v[j+m->cols*i];
+}
+
+
+void rtnl_mtx_set(const RationalMatrix *m, int i, int j, Rational v)
+{
+	assert(m != NULL);
+	m->v[j+m->cols*i] = v;
+}
+
+
+RationalMatrix *rtnl_mtx_from_intmat(const IntegerMatrix *m)
+{
+	RationalMatrix *n;
+	unsigned int rows, cols;
+	int i, j;
+
+	intmat_size(m, &rows, &cols);
+	n = rtnl_mtx_new(rows, cols);
+	if ( n == NULL ) return NULL;
+
+	for ( i=0; i<rows; i++ ) {
+		for ( j=0; j<cols; j++ ) {
+			rtnl_mtx_set(n, i, j, rtnl(intmat_get(m, i, j), 1));
+		}
+	}
+
+	return n;
+}
+
+
+IntegerMatrix *intmat_from_rtnl_mtx(const RationalMatrix *m)
+{
+	IntegerMatrix *n;
+	int i, j;
+
+	n = intmat_new(m->rows, m->cols);
+	if ( n == NULL ) return NULL;
+
+	for ( i=0; i<m->rows; i++ ) {
+		for ( j=0; j<m->cols; j++ ) {
+			Rational v = rtnl_mtx_get(m, i, j);
+			squish(&v);
+			if ( v.den != 1 ) {
+				ERROR("Rational matrix can't be converted to integers\n");
+				intmat_free(n);
+				return NULL;
+			}
+			intmat_set(n, i, j, v.num);
+		}
+	}
+
+	return n;
+}
+
+
+void rtnl_mtx_free(RationalMatrix *mtx)
+{
+	if ( mtx == NULL ) return;
+	free(mtx->v);
+	free(mtx);
+}
+
+
+/* rtnl_mtx_solve:
+ * @P: A %RationalMatrix
+ * @vec: An array of %Rational
+ * @ans: An array of %Rational in which to store the result
+ *
+ * Solves the matrix equation m*ans = vec, where @ans and @vec are
+ * vectors of %Rational.
+ *
+ * In this version, @m must be square.
+ *
+ * The number of columns in @m must equal the length of @ans, and the number
+ * of rows in @m must equal the length of @vec, but note that there is no way
+ * for this function to check that this is the case.
+ *
+ * Given that P is transformation of the unit cell axes (see ITA chapter 5.1),
+ * this function calculates the fractional coordinates of a point in the
+ * transformed axes, given its fractional coordinates in the original axes.
+ *
+ * Returns: non-zero on error
+ **/
+int transform_fractional_coords_rtnl(const RationalMatrix *P,
+                                     const Rational *ivec, Rational *ans)
+{
+	RationalMatrix *cm;
+	Rational *vec;
+	int h, k;
+	int i;
+
+	if ( P->rows != P->cols ) return 1;
+
+	/* Copy the matrix and vector because the calculation will
+	 * be done in-place */
+	cm = rtnl_mtx_copy(P);
+	if ( cm == NULL ) return 1;
+
+	vec = malloc(cm->rows*sizeof(Rational));
+	if ( vec == NULL ) return 1;
+	for ( h=0; h<cm->rows; h++ ) vec[h] = ivec[h];
+
+	/* Gaussian elimination with partial pivoting */
+	h = 0;
+	k = 0;
+	while ( h<=cm->rows && k<=cm->cols ) {
+
+		int prow = 0;
+		Rational pval = rtnl_zero();
+		Rational t;
+
+		/* Find the row with the largest value in column k */
+		for ( i=h; i<cm->rows; i++ ) {
+			if ( rtnl_cmp(rtnl_abs(rtnl_mtx_get(cm, i, k)), pval) > 0 ) {
+				pval = rtnl_abs(rtnl_mtx_get(cm, i, k));
+				prow = i;
+			}
+		}
+
+		if ( rtnl_cmp(pval, rtnl_zero()) == 0 ) {
+			k++;
+			continue;
+		}
+
+		/* Swap 'prow' with row h */
+		for ( i=0; i<cm->cols; i++ ) {
+			t = rtnl_mtx_get(cm, h, i);
+			rtnl_mtx_set(cm, h, i, rtnl_mtx_get(cm, prow, i));
+			rtnl_mtx_set(cm, prow, i, t);
+		}
+		t = vec[prow];
+		vec[prow] = vec[h];
+		vec[h] = t;
+
+		/* Divide and subtract rows below */
+		for ( i=h+1; i<cm->rows; i++ ) {
+
+			int j;
+			Rational dval;
+
+			dval = rtnl_div(rtnl_mtx_get(cm, i, k),
+			                rtnl_mtx_get(cm, h, k));
+
+			for ( j=0; j<cm->cols; j++ ) {
+				Rational t = rtnl_mtx_get(cm, i, j);
+				Rational p = rtnl_mul(dval, rtnl_mtx_get(cm, h, j));
+				t = rtnl_sub(t, p);
+				rtnl_mtx_set(cm, i, j, t);
+			}
+
+			/* Divide the right hand side as well */
+			Rational t = vec[i];
+			Rational p = rtnl_mul(dval, vec[h]);
+			vec[i] = rtnl_sub(t, p);
+		}
+
+		h++;
+		k++;
+
+
+	}
+
+	/* Back-substitution */
+	for ( i=cm->rows-1; i>=0; i-- ) {
+		int j;
+		Rational sum = rtnl_zero();
+		for ( j=i+1; j<cm->cols; j++ ) {
+			Rational av;
+			av = rtnl_mul(rtnl_mtx_get(cm, i, j), ans[j]);
+			sum = rtnl_add(sum, av);
+		}
+		sum = rtnl_sub(vec[i], sum);
+		ans[i] = rtnl_div(sum, rtnl_mtx_get(cm, i, i));
+	}
+
+	free(vec);
+	rtnl_mtx_free(cm);
+
+	return 0;
+}
+
+
+/**
+ * rtnl_mtx_print
+ * @m: A %RationalMatrix
+ *
+ * Prints @m to stderr.
+ *
+ */
+void rtnl_mtx_print(const RationalMatrix *m)
+{
+	unsigned int i, j;
+
+	for ( i=0; i<m->rows; i++ ) {
+
+		fprintf(stderr, "[ ");
+		for ( j=0; j<m->cols; j++ ) {
+			char *v = rtnl_format(rtnl_mtx_get(m, i, j));
+			fprintf(stderr, "%4s ", v);
+			free(v);
+		}
+		fprintf(stderr, "]\n");
+	}
+}
+
+
+void rtnl_mtx_mtxmult(const RationalMatrix *A, const RationalMatrix *B,
+                      RationalMatrix *ans)
+{
+	int i, j;
+
+	assert(ans->cols == B->cols);
+	assert(ans->rows == A->rows);
+	assert(A->cols == B->rows);
+
+	for ( i=0; i<ans->rows; i++ ) {
+		for ( j=0; j<ans->cols; j++ ) {
+			int k;
+			Rational sum = rtnl_zero();
+			for ( k=0; k<A->rows; k++ ) {
+				Rational add;
+				add = rtnl_mul(rtnl_mtx_get(A, i, k),
+				               rtnl_mtx_get(B, k, j));
+				sum = rtnl_add(sum, add);
+			}
+			rtnl_mtx_set(ans, i, j, sum);
+		}
+	}
+}
+
+
+/* Given a "P-matrix" (see ITA chapter 5.1), calculate the fractional
+ * coordinates of point "vec" in the original axes, given its fractional
+ * coordinates in the transformed axes.
+ * To transform in the opposite direction, we would multiply by Q = P^-1.
+ * Therefore, this direction is the "easy way" and needs just a matrix
+ * multiplication. */
+void transform_fractional_coords_rtnl_inverse(const RationalMatrix *P,
+                                              const Rational *vec,
+                                              Rational *ans)
+{
+	int i, j;
+
+	for ( i=0; i<P->rows; i++ ) {
+		ans[i] = rtnl_zero();
+		for ( j=0; j<P->cols; j++ ) {
+			Rational add;
+			add = rtnl_mul(rtnl_mtx_get(P, i, j), vec[j]);
+			ans[i] = rtnl_add(ans[i], add);
+		}
+	}
+}
+
+
+static RationalMatrix *delete_row_and_column(const RationalMatrix *m,
+                                             unsigned int di, unsigned int dj)
+{
+	RationalMatrix *n;
+	unsigned int i, j;
+
+	n = rtnl_mtx_new(m->rows-1, m->cols-1);
+	if ( n == NULL ) return NULL;
+
+	for ( i=0; i<n->rows; i++ ) {
+	for ( j=0; j<n->cols; j++ ) {
+
+		Rational val;
+		unsigned int gi, gj;
+
+		gi = (i>=di) ? i+1 : i;
+		gj = (j>=dj) ? j+1 : j;
+		val = rtnl_mtx_get(m, gi, gj);
+		rtnl_mtx_set(n, i, j, val);
+
+	}
+	}
+
+	return n;
+}
+
+
+static Rational cofactor(const RationalMatrix *m,
+                         unsigned int i, unsigned int j)
+{
+	RationalMatrix *n;
+	Rational t, C;
+
+	n = delete_row_and_column(m, i, j);
+	if ( n == NULL ) {
+		fprintf(stderr, "Failed to allocate matrix.\n");
+		return rtnl_zero();
+	}
+
+	/* -1 if odd, +1 if even */
+	t = (i+j) & 0x1 ? rtnl(-1, 1) : rtnl(1, 1);
+
+	C = rtnl_mul(t, rtnl_mtx_det(n));
+	rtnl_mtx_free(n);
+
+	return C;
+}
+
+
+
+Rational rtnl_mtx_det(const RationalMatrix *m)
+{
+	unsigned int i, j;
+	Rational det;
+
+	assert(m->rows == m->cols);  /* Otherwise determinant doesn't exist */
+
+	if ( m->rows == 2 ) {
+		Rational a, b;
+		a = rtnl_mul(rtnl_mtx_get(m, 0, 0), rtnl_mtx_get(m, 1, 1));
+		b = rtnl_mul(rtnl_mtx_get(m, 0, 1), rtnl_mtx_get(m, 1, 0));
+		return rtnl_sub(a, b);
+	}
+
+	i = 0;  /* Fixed */
+	det = rtnl_zero();
+	for ( j=0; j<m->cols; j++ ) {
+		Rational a;
+		a = rtnl_mul(rtnl_mtx_get(m, i, j), cofactor(m, i, j));
+		det = rtnl_add(det, a);
+	}
+
+	return det;
+}
diff --git a/libcrystfel/src/rational.h b/libcrystfel/src/rational.h
new file mode 100644
index 00000000..23c918cf
--- /dev/null
+++ b/libcrystfel/src/rational.h
@@ -0,0 +1,107 @@
+/*
+ * rational.h
+ *
+ * A small rational number library
+ *
+ * Copyright © 2019 Deutsches Elektronen-Synchrotron DESY,
+ *                  a research centre of the Helmholtz Association.
+ *
+ * Authors:
+ *   2019 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/>.
+ *
+ */
+
+#ifndef RATIONAL_H
+#define RATIONAL_H
+
+#ifdef HAVE_CONFIG_H
+#include <config.h>
+#endif
+
+/**
+ * Rational
+ *
+ * The Rational is an opaque-ish data structure representing a rational number.
+ *
+ * "Opaque-ish" means that the structure isn't technically opaque, allowing you
+ * to assign and allocate them easily.  But you shouldn't look at or set its
+ * contents, except by using the accessor functions.
+ **/
+typedef struct {
+	/* Private, don't modify */
+	signed long long int num;
+	signed long long int den;
+} Rational;
+
+
+/**
+ * RationalMatrix
+ *
+ * The RationalMatrix is an opaque data structure representing a matrix of
+ * rational numbers.
+ **/
+typedef struct _rationalmatrix RationalMatrix;
+
+#include "integer_matrix.h"
+
+#ifdef __cplusplus
+extern "C" {
+#endif
+
+extern Rational rtnl_zero(void);
+extern Rational rtnl(signed long long int num, signed long long int den);
+extern double rtnl_as_double(Rational r);
+
+extern Rational rtnl_mul(Rational a, Rational b);
+extern Rational rtnl_div(Rational a, Rational b);
+extern Rational rtnl_add(Rational a, Rational b);
+extern Rational rtnl_sub(Rational a, Rational b);
+
+extern signed int rtnl_cmp(Rational a, Rational b);
+extern Rational rtnl_abs(Rational a);
+
+extern char *rtnl_format(Rational rt);
+
+extern Rational *rtnl_list(signed int num_min, signed int num_max,
+                           signed int den_min, signed int den_max,
+                           int *pn);
+
+extern RationalMatrix *rtnl_mtx_new(unsigned int rows, unsigned int cols);
+extern RationalMatrix *rtnl_mtx_copy(const RationalMatrix *m);
+extern Rational rtnl_mtx_get(const RationalMatrix *m, int i, int j);
+extern void rtnl_mtx_set(const RationalMatrix *m, int i, int j, Rational v);
+extern RationalMatrix *rtnl_mtx_from_intmat(const IntegerMatrix *m);
+extern RationalMatrix *rtnl_mtx_identity(int rows);
+extern IntegerMatrix *intmat_from_rtnl_mtx(const RationalMatrix *m);
+extern void rtnl_mtx_free(RationalMatrix *mtx);
+extern void rtnl_mtx_mtxmult(const RationalMatrix *A, const RationalMatrix *B,
+                             RationalMatrix *ans);
+extern int transform_fractional_coords_rtnl(const RationalMatrix *P,
+                                            const Rational *ivec,
+                                            Rational *ans);
+extern void transform_fractional_coords_rtnl_inverse(const RationalMatrix *P,
+                                                     const Rational *vec,
+                                                     Rational *ans);
+extern void rtnl_mtx_print(const RationalMatrix *m);
+extern Rational rtnl_mtx_det(const RationalMatrix *m);
+
+#ifdef __cplusplus
+}
+#endif
+
+#endif	/* RATIONAL_H */
diff --git a/libcrystfel/src/symmetry.c b/libcrystfel/src/symmetry.c
index 06cc368c..69dada16 100644
--- a/libcrystfel/src/symmetry.c
+++ b/libcrystfel/src/symmetry.c
@@ -3,12 +3,12 @@
  *
  * Symmetry
  *
- * Copyright © 2012-2016 Deutsches Elektronen-Synchrotron DESY,
+ * Copyright © 2012-2019 Deutsches Elektronen-Synchrotron DESY,
  *                       a research centre of the Helmholtz Association.
  *
  * Authors:
- *   2010-2014,2016 Thomas White <taw@physics.org>
- *   2014           Kenneth Beyerlein <kenneth.beyerlein@desy.de>
+ *   2010-2019 Thomas White <taw@physics.org>
+ *   2014      Kenneth Beyerlein <kenneth.beyerlein@desy.de>
  *
  * This file is part of CrystFEL.
  *
@@ -41,6 +41,8 @@
 #include "symmetry.h"
 #include "utils.h"
 #include "integer_matrix.h"
+#include "symop-parse.h"
+#include "symop-lex.h"
 
 
 /**
@@ -195,9 +197,9 @@ static void add_symop_v(SymOpList *ops,
 	m = intmat_new(3, 3);
 	assert(m != NULL);
 
-	for ( i=0; i<3; i++ ) intmat_set(m, 0, i, h[i]);
-	for ( i=0; i<3; i++ ) intmat_set(m, 1, i, k[i]);
-	for ( i=0; i<3; i++ ) intmat_set(m, 2, i, l[i]);
+	for ( i=0; i<3; i++ ) intmat_set(m, i, 0, h[i]);
+	for ( i=0; i<3; i++ ) intmat_set(m, i, 1, k[i]);
+	for ( i=0; i<3; i++ ) intmat_set(m, i, 2, l[i]);
 
 	free(h);
 	free(k);
@@ -369,13 +371,13 @@ static void expand_ops(SymOpList *s)
 /* Transform all the operations in a SymOpList by a given matrix.
  * The matrix must have a determinant of +/- 1 (otherwise its inverse would
  * not also be an integer matrix). */
-static void transform_ops(SymOpList *s, IntegerMatrix *t)
+static void transform_ops(SymOpList *s, IntegerMatrix *P)
 {
 	int n, i;
-	IntegerMatrix *inv;
+	IntegerMatrix *Pi;
 	signed int det;
 
-	det = intmat_det(t);
+	det = intmat_det(P);
 	if ( det == -1 ) {
 		ERROR("WARNING: mirrored SymOpList.\n");
 	} else if ( det != 1 ) {
@@ -383,8 +385,8 @@ static void transform_ops(SymOpList *s, IntegerMatrix *t)
 		return;
 	}
 
-	inv = intmat_inverse(t);
-	if ( inv == NULL ) {
+	Pi = intmat_inverse(P);
+	if ( Pi == NULL ) {
 		ERROR("Failed to invert matrix.\n");
 		return;
 	}
@@ -394,13 +396,13 @@ static void transform_ops(SymOpList *s, IntegerMatrix *t)
 
 		IntegerMatrix *r, *f;
 
-		r = intmat_intmat_mult(s->ops[i], t);
+		r = intmat_intmat_mult(P, s->ops[i]);
 		if ( r == NULL ) {
 			ERROR("Matrix multiplication failed.\n");
 			return;
 		}
 
-		f = intmat_intmat_mult(inv, r);
+		f = intmat_intmat_mult(r, Pi);
 		if ( f == NULL ) {
 			ERROR("Matrix multiplication failed.\n");
 			return;
@@ -413,7 +415,7 @@ static void transform_ops(SymOpList *s, IntegerMatrix *t)
 
 	}
 
-	intmat_free(inv);
+	intmat_free(Pi);
 }
 
 
@@ -1012,14 +1014,14 @@ static SymOpList *getpg_arbitrary_ua(const char *sym, size_t s)
 
 		case 'a' :
 		intmat_set(t, 0, 2, 1);
-		intmat_set(t, 1, 1, 1);
-		intmat_set(t, 2, 0, -1);
+		intmat_set(t, 1, 0, 1);
+		intmat_set(t, 2, 1, 1);
 		break;
 
 		case 'b' :
-		intmat_set(t, 0, 0, 1);
+		intmat_set(t, 0, 1, 1);
 		intmat_set(t, 1, 2, 1);
-		intmat_set(t, 2, 1, -1);
+		intmat_set(t, 2, 0, 1);
 
 		break;
 
@@ -1124,7 +1126,7 @@ static void do_op(const IntegerMatrix *op,
 
 	v[0] = h;  v[1] = k;  v[2] = l;
 
-	ans = intmat_intvec_mult(op, v);
+	ans = transform_indices(op, v);
 	assert(ans != NULL);
 
 	*he = ans[0];  *ke = ans[1];  *le = ans[2];
@@ -1636,105 +1638,76 @@ SymOpList *get_ambiguities(const SymOpList *source, const SymOpList *target)
 }
 
 
-static IntegerMatrix *parse_symmetry_operation(const char *s)
+/* Parse a single symmetry operation, e.g. 'h,-2k,(h+l)/3' */
+RationalMatrix *parse_symmetry_operation(const char *s)
 {
-	IntegerMatrix *m;
-	char **els;
-	int n, i;
+	YY_BUFFER_STATE b;
+	RationalMatrix *m;
+	int r;
 
+	m = rtnl_mtx_new(3, 3);
+	b = symop_scan_string(s);
+	r = symopparse(m, NULL);
+	symop_delete_buffer(b);
 
-	n = assplode(s, ",", &els, ASSPLODE_NONE);
-	if ( n != 3 ) {
-		for ( i=0; i<n; i++ ) free(els[i]);
-		free(els);
+	if ( r ) {
+		ERROR("Failed to parse '%s'\n", s);
+		rtnl_mtx_free(m);
 		return NULL;
 	}
 
-	m = intmat_new(3, 3);
-	if ( m == NULL ) return NULL;
-
-	for ( i=0; i<n; i++ ) {
-
-		int c;
-		size_t cl;
-		signed int nh = 0;
-		signed int nk = 0;
-		signed int nl = 0;
-		signed int mult = 1;
-		int ndigit = 0;
-		signed int sign = +1;
-
-		/* We have one expression something like "-2h+k" */
-		cl = strlen(els[i]);
-		for ( c=0; c<cl; c++ ) {
-
-			if ( els[i][c] == '-' ) sign *= -1;
-			if ( els[i][c] == 'h' ) {
-				nh = mult*sign;
-				mult = 1;
-				ndigit = 0;
-				sign = +1;
-			}
-			if ( els[i][c] == 'k' ) {
-				nk = mult*sign;
-				mult = 1;
-				ndigit = 0;
-				sign = +1;
-			}
-			if ( els[i][c] == 'l' ) {
-				nl = mult*sign;
-				mult = 1;
-				ndigit = 0;
-				sign = +1;
-			}
-			if ( isdigit(els[i][c]) ) {
-				if ( ndigit > 0 ) {
-					mult *= 10;
-					mult += els[i][c] - '0';
-				} else {
-					mult *= els[i][c] - '0';
-				}
-				ndigit++;
-			}
-		}
-
-		intmat_set(m, i, 0, nh);
-		intmat_set(m, i, 1, nk);
-		intmat_set(m, i, 2, nl);
-
-		free(els[i]);
+	return m;
+}
 
-	}
-	free(els);
 
-	return m;
+/**
+ * parse_cell_transformation
+ * @s: Textual representation of cell transformation
+ *
+ * Parses @s, for example 'a,(b+a)/2,c', and returns the corresponding
+ * %RationalMatrix.
+ *
+ * Returns: A %RationalMatrix describing the transformation, or NULL on error.
+ *
+ */
+RationalMatrix *parse_cell_transformation(const char *s)
+{
+	return parse_symmetry_operation(s);
 }
 
 
+/**
+ * parse_symmetry_operations
+ * @s: Textual representation of a list of symmetry operations
+ *
+ * Parses @s, for example 'h,k,l;k,h,-l', and returns the corresponding
+ * %SymOpList
+ *
+ * Returns: A %SymOpList, or NULL on error.
+ *
+ */
 SymOpList *parse_symmetry_operations(const char *s)
 {
-	SymOpList *sol;
-	char **ops;
-	int n, i;
+	YY_BUFFER_STATE b;
+	RationalMatrix *m;
+	SymOpList *list;
+	int r;
 
-	sol = new_symoplist();
-	if ( sol == NULL ) return NULL;
+	m = rtnl_mtx_new(3, 3); /* Scratch space for parser */
+	list = new_symoplist(); /* The result we want */
 
-	n = assplode(s, ";:", &ops, ASSPLODE_NONE);
-	for ( i=0; i<n; i++ ) {
-		IntegerMatrix *m;
-		m = parse_symmetry_operation(ops[i]);
-		if ( m != NULL ) {
-			add_symop(sol, m);
-		} else {
-			ERROR("Invalid symmetry operation '%s'\n", ops[i]);
-			/* Try the next one */
-		}
-		free(ops[i]);
+	b = symop_scan_string(s);
+	r = symopparse(m, list);
+	symop_delete_buffer(b);
+	rtnl_mtx_free(m);
+
+	if ( r ) {
+		ERROR("Failed to parse '%s'\n", s);
+		free_symoplist(list);
+		return NULL;
 	}
-	free(ops);
 
-	return sol;
+	return list;
 }
 
 
diff --git a/libcrystfel/src/symmetry.h b/libcrystfel/src/symmetry.h
index a61ca96f..ab2aa934 100644
--- a/libcrystfel/src/symmetry.h
+++ b/libcrystfel/src/symmetry.h
@@ -3,12 +3,12 @@
  *
  * Symmetry
  *
- * Copyright © 2012-2016 Deutsches Elektronen-Synchrotron DESY,
+ * Copyright © 2012-2019 Deutsches Elektronen-Synchrotron DESY,
  *                       a research centre of the Helmholtz Association.
  *
  * Authors:
- *   2010-2014,2016 Thomas White <taw@physics.org>
- *   2014           Kenneth Beyerlein <kenneth.beyerlein@desy.de>
+ *   2010-2019 Thomas White <taw@physics.org>
+ *   2014      Kenneth Beyerlein <kenneth.beyerlein@desy.de>
  *
  * This file is part of CrystFEL.
  *
@@ -36,6 +36,7 @@
 
 
 #include "integer_matrix.h"
+#include "rational.h"
 
 /**
  * SymOpList
@@ -93,7 +94,9 @@ extern int is_centric(signed int h, signed int k, signed int l,
 extern void pointgroup_warning(const char *sym);
 
 extern void add_symop(SymOpList *ops, IntegerMatrix *m);
+extern RationalMatrix *parse_symmetry_operation(const char *s);
 extern SymOpList *parse_symmetry_operations(const char *s);
+extern RationalMatrix *parse_cell_transformation(const char *s);
 extern char *get_matrix_name(const IntegerMatrix *m, int row);
 
 #ifdef __cplusplus
diff --git a/libcrystfel/src/symop.l b/libcrystfel/src/symop.l
new file mode 100644
index 00000000..ed38a65e
--- /dev/null
+++ b/libcrystfel/src/symop.l
@@ -0,0 +1,52 @@
+/*
+ * symop.l
+ *
+ * Lexical scanner for symmetry operations
+ *
+ * Copyright © 2019 Deutsches Elektronen-Synchrotron DESY,
+ *                  a research centre of the Helmholtz Association.
+ *
+ * Authors:
+ *   2019 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/>.
+ *
+ */
+
+%{
+  #define YYDEBUG 1
+  #include "rational.h"
+  #include "symop-parse.h"
+%}
+
+%option prefix="symop"
+%option noyywrap nounput noinput
+
+%%
+
+[,]        { return COMMA; }
+[0-9]+     { symoplval.n = atoi(yytext); return NUMBER; }
+[/]        { return DIVIDE; }
+[+]        { return PLUS; }
+[-]        { return MINUS; }
+[ahx]      { return H; }
+[bky]      { return K; }
+[clz]      { return L; }
+[(]        { return OPENB; }
+[)]        { return CLOSEB; }
+[;]        { return SEMICOLON; }
+
+%%
diff --git a/libcrystfel/src/symop.y b/libcrystfel/src/symop.y
new file mode 100644
index 00000000..be31c21d
--- /dev/null
+++ b/libcrystfel/src/symop.y
@@ -0,0 +1,140 @@
+/*
+ * symop.y
+ *
+ * Parser for symmetry operations
+ *
+ * Copyright © 2019 Deutsches Elektronen-Synchrotron DESY,
+ *                  a research centre of the Helmholtz Association.
+ *
+ * Authors:
+ *   2019 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/>.
+ *
+ */
+
+%{
+  #include <stdio.h>
+
+  #include "rational.h"
+  #include "symmetry.h"
+
+  extern int symoplex();
+  extern int symopparse(RationalMatrix *m, SymOpList *list);
+  void symoperror(RationalMatrix *m, SymOpList *list, const char *s);
+%}
+
+%define api.prefix {symop}
+
+%parse-param {RationalMatrix *m} {SymOpList *list}
+
+%code requires {
+  #include "symmetry.h"
+}
+
+%union {
+  RationalMatrix *m;  /* Full rational matrix */
+  Rational rv[3];     /* Rational vector, e.g. '1/2h+3k' */
+  Rational r;         /* Rational number */
+  int n;              /* Just a number */
+}
+
+%token SEMICOLON
+%token COMMA
+%token NUMBER
+%token OPENB CLOSEB
+%token H K L
+
+%left PLUS MINUS
+%left DIVIDE
+%precedence MUL
+%precedence NEG
+
+%type <m> symop
+%type <rv> axexpr
+%type <rv> part
+%type <n> NUMBER
+%type <r> fraction
+
+%{
+static int try_add_symop(SymOpList *list, RationalMatrix *m, int complain)
+{
+	if ( list == NULL ) {
+		/* Only complain if this isn't the only operation provided */
+		if ( complain ) {
+			yyerror(m, list, "Must be a single symmetry operation");
+		}
+		return 1;
+	} else {
+		IntegerMatrix *im;
+		im = intmat_from_rtnl_mtx(m);
+		if ( im == NULL ) {
+			yyerror(m, list, "Symmetry operations must all be integer");
+			return 1;
+		} else {
+			add_symop(list, im);
+		}
+	}
+	return 0;
+}
+%}
+
+%%
+
+symoplist:
+  symop                       { try_add_symop(list, m, 0); }
+| symoplist SEMICOLON symop   { if ( try_add_symop(list, m, 1) ) YYERROR; }
+;
+
+symop:
+  axexpr COMMA axexpr COMMA axexpr { rtnl_mtx_set(m, 0, 0, $1[0]);
+                                     rtnl_mtx_set(m, 1, 0, $1[1]);
+                                     rtnl_mtx_set(m, 2, 0, $1[2]);
+                                     rtnl_mtx_set(m, 0, 1, $3[0]);
+                                     rtnl_mtx_set(m, 1, 1, $3[1]);
+                                     rtnl_mtx_set(m, 2, 1, $3[2]);
+                                     rtnl_mtx_set(m, 0, 2, $5[0]);
+                                     rtnl_mtx_set(m, 1, 2, $5[1]);
+                                     rtnl_mtx_set(m, 2, 2, $5[2]);
+                                   }
+;
+
+axexpr:
+  part                      { int i;  for ( i=0; i<3; i++ ) $$[i] = $1[i]; }
+| axexpr PLUS axexpr        { int i;  for ( i=0; i<3; i++ ) $$[i] = rtnl_add($1[i], $3[i]); }
+| axexpr MINUS axexpr       { int i;  for ( i=0; i<3; i++ ) $$[i] = rtnl_sub($1[i], $3[i]); }
+| MINUS axexpr %prec NEG    { int i;  for ( i=0; i<3; i++ ) $$[i] = rtnl_sub(rtnl_zero(), $2[i]); }
+| OPENB axexpr CLOSEB       { int i;  for ( i=0; i<3; i++ ) $$[i] = $2[i]; }
+| axexpr DIVIDE NUMBER      { int i;  for ( i=0; i<3; i++ ) $$[i] = rtnl_div($1[i], rtnl($3, 1)); }
+| NUMBER axexpr %prec MUL   { int i;  for ( i=0; i<3; i++ ) $$[i] = rtnl_mul($2[i], rtnl($1, 1)); }
+| fraction axexpr %prec MUL { int i;  for ( i=0; i<3; i++ ) $$[i] = rtnl_mul($2[i], $1); }
+;
+
+part:
+  H  { $$[0] = rtnl(1, 1); $$[1] = rtnl_zero(); $$[2] = rtnl_zero(); }
+| K  { $$[1] = rtnl(1, 1); $$[0] = rtnl_zero(); $$[2] = rtnl_zero(); }
+| L  { $$[2] = rtnl(1, 1); $$[0] = rtnl_zero(); $$[1] = rtnl_zero(); }
+;
+
+fraction:
+  NUMBER DIVIDE NUMBER     { $$ = rtnl($1, $3); }
+;
+
+%%
+
+void symoperror(RationalMatrix *m, SymOpList *list, const char *s) {
+	printf("Error: %s\n", s);
+}
diff --git a/libcrystfel/src/taketwo.c b/libcrystfel/src/taketwo.c
index 203e5a05..933fa76a 100644
--- a/libcrystfel/src/taketwo.c
+++ b/libcrystfel/src/taketwo.c
@@ -98,6 +98,7 @@
 #include <assert.h>
 #include <time.h>
 
+#include "cell.h"
 #include "cell-utils.h"
 #include "index.h"
 #include "taketwo.h"
@@ -355,6 +356,54 @@ static void show_rvec(struct rvec r2)
  * functions called under the core functions, still specialised (Level 3)
  * ------------------------------------------------------------------------*/
 
+/* cell_transform_gsl_direct() doesn't do quite what we want here.
+ * The matrix m should be post-multiplied by a matrix of real or reciprocal
+ * basis vectors (it doesn't matter which because it's just a rotation).
+ * M contains the basis vectors written in columns:  M' = mM */
+static UnitCell *cell_post_smiley_face(UnitCell *in, gsl_matrix *m)
+{
+	gsl_matrix *c;
+	double asx, asy, asz;
+	double bsx, bsy, bsz;
+	double csx, csy, csz;
+	gsl_matrix *res;
+	UnitCell *out;
+
+	cell_get_cartesian(in, &asx, &asy, &asz,
+	                       &bsx, &bsy, &bsz,
+	                       &csx, &csy, &csz);
+
+	c = gsl_matrix_alloc(3, 3);
+	gsl_matrix_set(c, 0, 0, asx);
+	gsl_matrix_set(c, 1, 0, asy);
+	gsl_matrix_set(c, 2, 0, asz);
+	gsl_matrix_set(c, 0, 1, bsx);
+	gsl_matrix_set(c, 1, 1, bsy);
+	gsl_matrix_set(c, 2, 1, bsz);
+	gsl_matrix_set(c, 0, 2, csx);
+	gsl_matrix_set(c, 1, 2, csy);
+	gsl_matrix_set(c, 2, 2, csz);
+
+	res = gsl_matrix_calloc(3, 3);
+	gsl_blas_dgemm(CblasNoTrans, CblasNoTrans, 1.0, m, c, 0.0, res);
+
+	out = cell_new_from_cell(in);
+	cell_set_cartesian(out, gsl_matrix_get(res, 0, 0),
+	                        gsl_matrix_get(res, 1, 0),
+	                        gsl_matrix_get(res, 2, 0),
+	                        gsl_matrix_get(res, 0, 1),
+	                        gsl_matrix_get(res, 1, 1),
+	                        gsl_matrix_get(res, 2, 1),
+	                        gsl_matrix_get(res, 0, 2),
+	                        gsl_matrix_get(res, 1, 2),
+	                        gsl_matrix_get(res, 2, 2));
+
+	gsl_matrix_free(res);
+	gsl_matrix_free(c);
+	return out;
+}
+
+
 static void rotation_around_axis(struct rvec c, double th,
 				 gsl_matrix *res)
 {
@@ -463,35 +512,6 @@ static void closest_rot_mat(struct rvec vec1, struct rvec vec2,
 	rotation_around_axis(axis, bestAngle, twizzle);
 }
 
-/*
-static double matrix_angle(gsl_matrix *m)
-{
-	double a = gsl_matrix_get(m, 0, 0);
-	double b = gsl_matrix_get(m, 1, 1);
-	double c = gsl_matrix_get(m, 2, 2);
-
-	double cos_t = (a + b + c - 1) / 2;
-	double theta = acos(cos_t);
-
-	return theta;
-}
-
-static struct rvec matrix_axis(gsl_matrix *a)
-{
-	double ang = matrix_angle(a);
-	double cos_t = cos(ang);
-	double p = gsl_matrix_get(a, 0, 0);
-	double q = gsl_matrix_get(a, 1, 1);
-	double r = gsl_matrix_get(a, 2, 2);
-	double x = sqrt((p - cos_t) / (1 - cos_t));
-	double y = sqrt((q - cos_t) / (1 - cos_t));
-	double z = sqrt((r - cos_t) / (1 - cos_t));
-
-	struct rvec v = new_rvec(x, y, z);
-	return v;
-}
-*/
-
 static double matrix_trace(gsl_matrix *a)
 {
 	int i;
@@ -1601,7 +1621,8 @@ static unsigned int start_seeds(gsl_matrix **rotation, struct TakeTwoCell *cell)
 }
 
 
-static void set_gsl_matrix(gsl_matrix *mat, double asx, double asy, double asz,
+static void set_gsl_matrix(gsl_matrix *mat,
+                           double asx, double asy, double asz,
                            double bsx, double bsy, double bsz,
                            double csx, double csy, double csz)
 {
@@ -1630,15 +1651,23 @@ static int generate_rotation_sym_ops(struct TakeTwoCell *ttCell)
 
 	gsl_matrix *recip = gsl_matrix_alloc(3, 3);
 	gsl_matrix *cart = gsl_matrix_alloc(3, 3);
-	cell_get_reciprocal(ttCell->cell, &asx, &asy, &asz, &bsx, &bsy,
-						&bsz, &csx, &csy, &csz);
 
-	set_gsl_matrix(recip, asx, asy, asz, bsx, bsy, bsz, csx, csy, csz);
+	cell_get_reciprocal(ttCell->cell, &asx, &asy, &asz,
+	                                  &bsx, &bsy, &bsz,
+	                                  &csx, &csy, &csz);
+
+	set_gsl_matrix(recip, asx, asy, asz,
+	                      asx, bsy, bsz,
+	                      csx, csy, csz);
+
+	cell_get_cartesian(ttCell->cell, &asx, &asy, &asz,
+	                                 &bsx, &bsy, &bsz,
+	                                 &csx, &csy, &csz);
 
-	cell_get_cartesian(ttCell->cell, &asx, &asy, &asz, &bsx, &bsy,
-						&bsz, &csx, &csy, &csz);
+	set_gsl_matrix(cart, asx, bsx, csx,
+	                     asy, bsy, csy,
+	                     asz, bsz, csz);
 
-	set_gsl_matrix(cart, asx, asy, asz, bsx, bsy, bsz, csx, csy, csz);
 
 	int i, j, k;
 	int numOps = num_equivs(rawList, NULL);
@@ -2058,7 +2087,7 @@ static UnitCell *run_taketwo(UnitCell *cell, const struct taketwo_options *opts,
 	tp->numPrevs++;
 
 	/* Prepare the solution for CrystFEL friendliness */
-	result = transform_cell_gsl(cell, solution);
+	result = cell_post_smiley_face(cell, solution);
 	cleanup_taketwo_cell(&ttCell);
 
 	return result;
diff --git a/libcrystfel/src/xgandalf.c b/libcrystfel/src/xgandalf.c
index 5f31df33..51819956 100644
--- a/libcrystfel/src/xgandalf.c
+++ b/libcrystfel/src/xgandalf.c
@@ -46,7 +46,7 @@ struct xgandalf_private_data {
 	IndexingMethod indm;
 	UnitCell *cellTemplate;
 	Lattice_t sampleRealLattice_A;   //same as cellTemplate
-	UnitCellTransformation *uncenteringTransformation;
+	IntegerMatrix *centeringTransformation;
 	LatticeTransform_t latticeReductionTransform;
 };
 
@@ -114,7 +114,7 @@ int run_xgandalf(struct image *image, void *ipriv)
 		if(xgandalf_private_data->cellTemplate != NULL){
 			restoreCell(uc, &xgandalf_private_data->latticeReductionTransform);
 
-			UnitCell *new_cell_trans = cell_transform_inverse(uc, xgandalf_private_data->uncenteringTransformation);
+			UnitCell *new_cell_trans = cell_transform_intmat(uc, xgandalf_private_data->centeringTransformation);
 			cell_free(uc);
 			uc = new_cell_trans;
 
@@ -147,7 +147,7 @@ void *xgandalf_prepare(IndexingMethod *indm, UnitCell *cell,
 	allocReciprocalPeaks(&(xgandalf_private_data->reciprocalPeaks_1_per_A));
 	xgandalf_private_data->indm = *indm;
 	xgandalf_private_data->cellTemplate = NULL;
-	xgandalf_private_data->uncenteringTransformation = NULL;
+	xgandalf_private_data->centeringTransformation = NULL;
 
 	float tolerance = xgandalf_opts->tolerance;
 	samplingPitch_t samplingPitch = xgandalf_opts->sampling_pitch;
@@ -157,7 +157,7 @@ void *xgandalf_prepare(IndexingMethod *indm, UnitCell *cell,
 
 		xgandalf_private_data->cellTemplate = cell;
 
-		UnitCell* primitiveCell = uncenter_cell(cell, &xgandalf_private_data->uncenteringTransformation);
+		UnitCell* primitiveCell = uncenter_cell(cell, &xgandalf_private_data->centeringTransformation, NULL);
 
 		UnitCell *uc = cell_new_from_cell(primitiveCell);
 		reduceCell(primitiveCell, &xgandalf_private_data->latticeReductionTransform);
@@ -250,8 +250,8 @@ void xgandalf_cleanup(void *pp)
 
 	freeReciprocalPeaks(xgandalf_private_data->reciprocalPeaks_1_per_A);
 	IndexerPlain_delete(xgandalf_private_data->indexer);
-	if(xgandalf_private_data->uncenteringTransformation != NULL){
-		tfn_free(xgandalf_private_data->uncenteringTransformation);
+	if(xgandalf_private_data->centeringTransformation != NULL){
+		intmat_free(xgandalf_private_data->centeringTransformation);
 	}
 	free(xgandalf_private_data);
 }