Handle weird unique axis variants of all point groups

5 files changed, 703 insertions, 41 deletions

diff --git a/libcrystfel/Makefile.am b/libcrystfel/Makefile.am
index 28d483ac..2a02a20b 100644
--- a/libcrystfel/Makefile.am
+++ b/libcrystfel/Makefile.am
@@ -8,7 +8,7 @@ libcrystfel_la_SOURCES = src/reflist.c src/utils.c src/cell.c src/detector.c \
                          src/symmetry.c src/stream.c src/peaks.c \
                          src/reflist-utils.c src/filters.c \
                          src/render.c src/index.c src/dirax.c src/mosflm.c \
-                         src/cell-utils.c
+                         src/cell-utils.c src/integer_matrix.c
 
 if HAVE_FFTW
 libcrystfel_la_SOURCES += src/reax.c
@@ -23,7 +23,8 @@ libcrystfel_la_include_HEADERS = src/beam-parameters.h src/hdf5-file.h \
                                  src/geometry.h src/peaks.h src/stream.h \
                                  src/render.h src/index.h src/image.h \
                                  src/filters.h src/dirax.h src/mosflm.h \
-                                 src/index-priv.h src/reax.h src/cell-utils.h
+                                 src/index-priv.h src/reax.h src/cell-utils.h \
+                                 src/integer_matrix.h
 
 INCLUDES = "-I$(top_srcdir)/data"
 AM_CPPFLAGS = -DDATADIR=\""$(datadir)"\" -I$(top_builddir)/lib -Wall
diff --git a/libcrystfel/src/integer_matrix.c b/libcrystfel/src/integer_matrix.c
new file mode 100644
index 00000000..5aefb958
--- /dev/null
+++ b/libcrystfel/src/integer_matrix.c
@@ -0,0 +1,386 @@
+/*
+ * integer_matrix.c
+ *
+ * A small integer matrix library
+ *
+ * Copyright © 2012 Deutsches Elektronen-Synchrotron DESY,
+ *                  a research centre of the Helmholtz Association.
+ *
+ * Authors:
+ *   2012 Thomas White <taw@physics.org>
+ *
+ * This file is part of CrystFEL.
+ *
+ * CrystFEL is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU General Public License as published by
+ * the Free Software Foundation, either version 3 of the License, or
+ * (at your option) any later version.
+ *
+ * CrystFEL is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU General Public License
+ * along with CrystFEL.  If not, see <http://www.gnu.org/licenses/>.
+ *
+ */
+
+#ifdef HAVE_CONFIG_H
+#include <config.h>
+#endif
+
+#include <stdlib.h>
+#include <stdio.h>
+#include <string.h>
+#include <assert.h>
+
+#include "integer_matrix.h"
+
+
+/**
+ * SECTION:integer_matrix
+ * @short_description: Integer matrices
+ * @title: Integer matrices
+ * @section_id:
+ * @see_also:
+ * @include: "integer_matrix.h"
+ * @Image:
+ *
+ * Routines to handle integer matrix
+ */
+
+
+struct _integermatrix
+{
+	unsigned int rows;
+	unsigned int cols;
+
+	signed int *v;
+};
+
+
+/**
+ * intmat_new:
+ *
+ * Allocates a new %IntegerMatrix with all elements set to zero.
+ *
+ * Returns: a new %IntegerMatrix, or NULL on error.
+ **/
+IntegerMatrix *intmat_new(unsigned int rows, unsigned int cols)
+{
+	IntegerMatrix *m;
+
+	m = malloc(sizeof(IntegerMatrix));
+	if ( m == NULL ) return NULL;
+
+	m->v = calloc(rows*cols, sizeof(signed int));
+	if ( m->v == NULL ) {
+		free(m);
+		return NULL;
+	}
+
+	m->rows = rows;
+	m->cols = cols;
+
+	return m;
+}
+
+
+/**
+ * intmat_free:
+ * @m: An %IntegerMatrix
+ *
+ * Frees @m.
+ **/
+void intmat_free(IntegerMatrix *m)
+{
+	free(m->v);
+	free(m);
+}
+
+
+/**
+ * intmat_set:
+ * @m: An %IntegerMatrix
+ * @i: row number to set
+ * @j: column number to set
+ * @v: value to set to
+ *
+ * Sets the @i,@j element of @m to @v.
+ **/
+void intmat_set(IntegerMatrix *m, unsigned int i, unsigned int j, signed int v)
+{
+	assert(i < m->rows);
+	assert(j < m->cols);
+	m->v[j + m->cols*i] = v;
+}
+
+
+/**
+ * intmat_get:
+ * @m: An %IntegerMatrix
+ * @i: column number to set
+ * @j: row number to set
+ *
+ * Gets the @i,@j element of @m.
+ *
+ * Returns: the @i,@j element of @m.
+ **/
+signed int intmat_get(IntegerMatrix *m, unsigned int i, unsigned int j)
+{
+	assert(i < m->rows);
+	assert(j < m->cols);
+	return m->v[j + m->cols*i];
+}
+
+
+/**
+ * intmat_intvec_mult:
+ * @m: An %IntegerMatrix
+ * @vec: An array of signed integers
+ *
+ * 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.
+ *
+ * Returns: a newly allocated array of signed integers containing the answer,
+ * or NULL on error.
+ **/
+signed int *intmat_intvec_mult(IntegerMatrix *m, signed int *vec)
+{
+	signed int *ans;
+	unsigned int i;
+
+	ans = malloc(m->rows * sizeof(signed int));
+	if ( ans == NULL ) return NULL;
+
+	for ( i=0; i<m->rows; i++ ) {
+
+		unsigned int j;
+
+		ans[i] = 0;
+		for ( j=0; j<m->cols; j++ ) {
+			ans[i] += intmat_get(m, i, j) * vec[j];
+		}
+
+	}
+
+	return ans;
+}
+
+
+/**
+ * intmat_intmat_mult:
+ * @a: An %IntegerMatrix
+ * @b: An %IntegerMatrix
+ *
+ * Multiplies the matrix @a by the matrix @b.
+ *
+ * Returns: a newly allocated %IntegerMatrix containing the answer, or NULL on
+ * error.
+ **/
+IntegerMatrix *intmat_intmat_mult(IntegerMatrix *a, IntegerMatrix *b)
+{
+	unsigned int i, j;
+	IntegerMatrix *ans;
+
+	if ( a->cols != b->rows ) return NULL;
+
+	ans = intmat_new(a->rows, a->cols);
+	if ( ans == NULL ) return NULL;
+
+	for ( i=0; i<ans->rows; i++ ) {
+	for ( j=0; j<ans->cols; j++ ) {
+
+		unsigned int k;
+		signed int r = 0;
+
+		for ( k=0; k<a->cols; k++ ) {  /* a->cols == b->rows */
+			r += intmat_get(a, i, k) * intmat_get(b, k, j);
+		}
+		intmat_set(ans, i, j, r);
+
+	}
+	}
+
+	return ans;
+}
+
+
+static IntegerMatrix *delete_row_and_column(IntegerMatrix *m,
+                                            unsigned int di, unsigned int dj)
+{
+	IntegerMatrix *n;
+	unsigned int i, j;
+
+	n = intmat_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++ ) {
+
+		signed int val;
+		unsigned int gi, gj;
+
+		gi = (i>=di) ? i+1 : i;
+		gj = (j>=dj) ? j+1 : j;
+		val = intmat_get(m, gi, gj);
+		intmat_set(n, i, j, val);
+
+	}
+	}
+
+	return n;
+}
+
+
+/**
+ * intmat_det:
+ * @m: An %IntegerMatrix
+ *
+ * Calculates the determinant of @m.  Inefficiently.
+ *
+ * Returns: the determinant of @m.
+ **/
+signed int intmat_det(IntegerMatrix *m)
+{
+	unsigned int i, j;
+	signed int det = 0;
+
+	assert(m->rows == m->cols);  /* Otherwise determinant doesn't exist */
+
+	if ( m->rows == 2 ) {
+		return intmat_get(m, 0, 0)*intmat_get(m, 1, 1)
+		     - intmat_get(m, 0, 1)*intmat_get(m, 1, 0);
+	}
+
+	i = 0;  /* Fixed */
+	for ( j=0; j<m->cols; j++ ) {
+
+		signed int C;
+		IntegerMatrix *n;
+		signed int t;
+
+		n = delete_row_and_column(m, i, j);
+		if ( n == NULL ) {
+			fprintf(stderr, "Failed to allocate matrix.\n");
+			return 0;
+		}
+
+		/* -1 if odd, +1 if even */
+		t = (i+j) & 0x1 ? -1 : +1;
+
+		C = t * intmat_det(n);
+		intmat_free(n);
+
+		det += intmat_get(m, i, j) * C;
+
+	}
+
+	return det;
+}
+
+
+static IntegerMatrix *intmat_cofactors(IntegerMatrix *m)
+{
+	IntegerMatrix *n;
+	signed int i, j;
+
+	n = intmat_new(m->rows, m->cols);
+	if ( n == NULL ) return NULL;
+
+	for ( i=0; i<n->rows; i++ ) {
+	for ( j=0; j<n->cols; j++ ) {
+
+		signed int C;
+		IntegerMatrix *M;
+		signed int t;
+
+		M = delete_row_and_column(m, i, j);
+		if ( M == NULL ) {
+			fprintf(stderr, "Failed to allocate matrix.\n");
+			return 0;
+		}
+
+		/* -1 if odd, +1 if even */
+		t = (i+j) & 0x1 ? -1 : +1;
+
+		C = t * intmat_det(M);
+		intmat_free(M);
+
+		intmat_set(n, i, j, C);
+
+	}
+	}
+
+	return n;
+}
+
+
+/**
+ * intmat_inverse:
+ * @m: An %IntegerMatrix
+ *
+ * Calculates the inverse of @m.  Inefficiently.
+ *
+ * Returns: the inverse of @m, or NULL on error.
+ **/
+IntegerMatrix *intmat_inverse(IntegerMatrix *m)
+{
+	IntegerMatrix *adjugateT;
+	IntegerMatrix *inverse;
+	unsigned int i, j;
+	signed int det;
+
+	det = intmat_det(m);
+	if ( (det != +1) && (det != -1) ) {
+		fprintf(stderr, "Inverse matrix not an integer matrix.\n");
+		return NULL;
+	}
+
+	adjugateT = intmat_cofactors(m);
+	if ( adjugateT == NULL ) return NULL;
+
+	inverse = intmat_new(m->cols, m->rows);  /* The other way round */
+	if ( inverse == NULL ) return NULL;
+
+	for ( i=0; i<inverse->rows; i++ ) {
+	for ( j=0; j<inverse->cols; j++ ) {
+
+		signed int v;
+
+		v = intmat_get(adjugateT, j, i);
+
+		/* 1/-1 = -1 and 1/+1 = +1, and these are the only two cases */
+		intmat_set(inverse, i, j, v*det);
+
+	}
+	}
+
+	intmat_free(adjugateT);
+
+	return inverse;
+}
+
+
+/**
+ * intmat_print
+ * @m: An %IntegerMatrix
+ *
+ * Prints @m to stderr.
+ *
+ */
+void intmat_print(IntegerMatrix *m)
+{
+	unsigned int i, j;
+
+	for ( i=0; i<m->rows; i++ ) {
+
+		fprintf(stderr, "[ ");
+		for ( j=0; j<m->cols; j++ ) {
+			fprintf(stderr, "%4i ", intmat_get(m, i, j));
+		}
+		fprintf(stderr, "]\n");
+	}
+}
diff --git a/libcrystfel/src/integer_matrix.h b/libcrystfel/src/integer_matrix.h
new file mode 100644
index 00000000..41780629
--- /dev/null
+++ b/libcrystfel/src/integer_matrix.h
@@ -0,0 +1,68 @@
+/*
+ * integer_matrix.h
+ *
+ * A small integer matrix library
+ *
+ * Copyright © 2012 Deutsches Elektronen-Synchrotron DESY,
+ *                  a research centre of the Helmholtz Association.
+ *
+ * Authors:
+ *   2012 Thomas White <taw@physics.org>
+ *
+ * This file is part of CrystFEL.
+ *
+ * CrystFEL is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU General Public License as published by
+ * the Free Software Foundation, either version 3 of the License, or
+ * (at your option) any later version.
+ *
+ * CrystFEL is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU General Public License
+ * along with CrystFEL.  If not, see <http://www.gnu.org/licenses/>.
+ *
+ */
+
+#ifndef INTEGER_MATRIX_H
+#define INTEGER_MATRIX_H
+
+#ifdef HAVE_CONFIG_H
+#include <config.h>
+#endif
+
+/**
+ * IntegerMatrix
+ *
+ * The IntegerMatrix is an opaque data structure representing an integer matrix.
+ **/
+typedef struct _integermatrix IntegerMatrix;
+
+
+/* Alloc/dealloc */
+extern IntegerMatrix *intmat_new(unsigned int rows, unsigned int cols);
+extern void intmat_free(IntegerMatrix *m);
+
+/* Get/set */
+extern void intmat_set(IntegerMatrix *m, unsigned int i, unsigned int j,
+                       signed int v);
+extern signed int intmat_get(IntegerMatrix *m, unsigned int i, unsigned int j);
+
+/* Matrix-(int)vector multiplication */
+extern signed int *intmat_intvec_mult(IntegerMatrix *m, signed int *vec);
+
+/* Matrix-matrix multiplication */
+extern IntegerMatrix *intmat_intmat_mult(IntegerMatrix *a, IntegerMatrix *b);
+
+/* Inverse */
+extern IntegerMatrix *intmat_inverse(IntegerMatrix *m);
+
+/* Determinant */
+extern signed int intmat_det(IntegerMatrix *m);
+
+/* Diagnostics */
+extern void intmat_print(IntegerMatrix *m);
+
+#endif	/* INTEGER_MATRIX_H */
diff --git a/libcrystfel/src/symmetry.c b/libcrystfel/src/symmetry.c
index 52a9aae6..14681434 100644
--- a/libcrystfel/src/symmetry.c
+++ b/libcrystfel/src/symmetry.c
@@ -38,6 +38,7 @@
 
 #include "symmetry.h"
 #include "utils.h"
+#include "integer_matrix.h"
 
 
 /**
@@ -57,7 +58,7 @@ struct sym_op
 {
 	signed int *h;
 	signed int *k;
-	signed int *l;  /* Contributions to h, k and l from h, k, i and l */
+	signed int *l;  /* Contributions to h, k and l from h, k and l */
 	int order;
 };
 
@@ -359,6 +360,100 @@ static SymOpList *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, signed int *na,
+                                        signed int *nb,
+                                        signed int *nc)
+{
+	int n, i;
+	IntegerMatrix *t, *inv;
+	signed int det;
+
+	t = intmat_new(3, 3);
+	if ( t == NULL ) {
+		ERROR("Failed to allocate matrix.\n");
+		return;
+	}
+
+	intmat_set(t, 0, 0, na[0]);
+	intmat_set(t, 1, 0, na[1]);
+	intmat_set(t, 2, 0, na[2]);
+	intmat_set(t, 0, 1, nb[0]);
+	intmat_set(t, 1, 1, nb[1]);
+	intmat_set(t, 2, 1, nb[2]);
+	intmat_set(t, 0, 2, nc[0]);
+	intmat_set(t, 1, 2, nc[1]);
+	intmat_set(t, 2, 2, nc[2]);
+
+	det = intmat_det(t);
+	if ( det == -1 ) {
+		ERROR("Warning: mirrored SymOpList.\n");
+	} else if ( det != 1 ) {
+		ERROR("Invalid transformation for SymOpList.\n");
+		return;
+	}
+
+	inv = intmat_inverse(t);
+	if ( inv == NULL ) {
+		ERROR("Failed to invert matrix.\n");
+		return;
+	}
+
+	n = num_ops(s);
+	for ( i=0; i<n; i++ ) {
+
+		IntegerMatrix *m, *r, *f;
+
+		m = intmat_new(3, 3);
+		if ( m == NULL ) {
+			ERROR("Failed to allocate matrix.\n");
+			return;
+		}
+
+		intmat_set(m, 0, 0, s->ops[i].h[0]);
+		intmat_set(m, 1, 0, s->ops[i].h[1]);
+		intmat_set(m, 2, 0, s->ops[i].h[2]);
+		intmat_set(m, 0, 1, s->ops[i].k[0]);
+		intmat_set(m, 1, 1, s->ops[i].k[1]);
+		intmat_set(m, 2, 1, s->ops[i].k[2]);
+		intmat_set(m, 0, 2, s->ops[i].l[0]);
+		intmat_set(m, 1, 2, s->ops[i].l[1]);
+		intmat_set(m, 2, 2, s->ops[i].l[2]);
+
+		r = intmat_intmat_mult(m, t);
+		if ( r == NULL ) {
+			ERROR("Matrix multiplication failed.\n");
+			return;
+		}
+		intmat_free(m);
+
+		f = intmat_intmat_mult(inv, r);
+		if ( f == NULL ) {
+			ERROR("Matrix multiplication failed.\n");
+			return;
+		}
+		intmat_free(r);
+
+		s->ops[i].h[0] = intmat_get(f, 0, 0);
+		s->ops[i].h[1] = intmat_get(f, 1, 0);
+		s->ops[i].h[2] = intmat_get(f, 2, 0);
+		s->ops[i].k[0] = intmat_get(f, 0, 1);
+		s->ops[i].k[1] = intmat_get(f, 1, 1);
+		s->ops[i].k[2] = intmat_get(f, 2, 1);
+		s->ops[i].l[0] = intmat_get(f, 0, 2);
+		s->ops[i].l[1] = intmat_get(f, 1, 2);
+		s->ops[i].l[2] = intmat_get(f, 2, 2);
+		intmat_free(f);
+
+	}
+
+	intmat_free(t);
+	intmat_free(inv);
+}
+
+
 /********************************* Triclinic **********************************/
 
 static SymOpList *make_1bar()
@@ -390,15 +485,6 @@ static SymOpList *make_2m()
 }
 
 
-static SymOpList *make_2_uab()
-{
-	SymOpList *new = new_symoplist();
-	add_symop(new, v(-1,0,0,0), v(0,1,0,0), v(0,0,0,-1), 2); /* 2 // k */
-	new->name = strdup("2_uab");
-	return expand_ops(new);
-}
-
-
 static SymOpList *make_2()
 {
 	SymOpList *new = new_symoplist();
@@ -531,17 +617,6 @@ static SymOpList *make_4mmm()
 }
 
 
-static SymOpList *make_4mmm_uaa()
-{
-	SymOpList *new = new_symoplist();
-	add_symop(new, v(1,0,0,0), v(0,0,0,1), v(0,-1,0,0), 4); /* 4 // h */
-	add_symop(new, v(1,0,0,0), v(0,-1,0,0), v(0,0,0,1), 2); /* m -| k */
-	add_symop(new, v(-1,0,0,0), v(0,1,0,0), v(0,0,0,1), 2); /* m -| h */
-	new->name = strdup("4/mmm_uaa");
-	return expand_ops(new);
-}
-
-
 /************************** Trigonal (Rhombohedral) ***************************/
 
 static SymOpList *make_3_R()
@@ -817,21 +892,7 @@ static SymOpList *make_m3barm()
 }
 
 
-/**
- * get_pointgroup:
- * @sym: A string representation of a point group
- *
- * This function parses @sym and returns the corresponding %SymOpList.
- * In the string representation of the point group, use a preceding minus sign
- * for any character which would have a "bar".  Trigonal groups must be suffixed
- * with either "_H" or "_R" for a hexagonal or rhombohedral lattice
- * respectively.
- *
- * Examples: -1 1 2/m 2 m mmm 222 mm2 4/m 4 -4 4/mmm 422 -42m -4m2 4mm
- * 3_R -3_R 32_R 3m_R -3m_R 3_H -3_H 321_H 312_H 3m1_H 31m_H -3m1_H -31m_H
- * 6/m 6 -6 6/mmm 622 -62m -6m2 6mm 23 m-3 432 -43m m-3m.
- **/
-SymOpList *get_pointgroup(const char *sym)
+static SymOpList *getpg_old(const char *sym)
 {
 	/* Triclinic */
 	if ( strcmp(sym, "-1") == 0 ) return make_1bar();
@@ -840,7 +901,6 @@ SymOpList *get_pointgroup(const char *sym)
 	/* Monoclinic */
 	if ( strcmp(sym, "2/m") == 0 ) return make_2m();
 	if ( strcmp(sym, "2") == 0 ) return make_2();
-	if ( strcmp(sym, "2_uab") == 0 ) return make_2_uab();
 	if ( strcmp(sym, "m") == 0 ) return make_m();
 
 	/* Orthorhombic */
@@ -853,7 +913,6 @@ SymOpList *get_pointgroup(const char *sym)
 	if ( strcmp(sym, "4") == 0 ) return make_4();
 	if ( strcmp(sym, "-4") == 0 ) return make_4bar();
 	if ( strcmp(sym, "4/mmm") == 0 ) return make_4mmm();
-	if ( strcmp(sym, "4/mmm_uaa") == 0 ) return make_4mmm_uaa();
 	if ( strcmp(sym, "422") == 0 ) return make_422();
 	if ( strcmp(sym, "-42m") == 0 ) return make_4bar2m();
 	if ( strcmp(sym, "-4m2") == 0 ) return make_4barm2();
@@ -898,6 +957,136 @@ SymOpList *get_pointgroup(const char *sym)
 }
 
 
+static int char_count(const char *a, char b)
+{
+	size_t i;
+	int n;
+
+	i = 0;  n = 0;
+	do {
+		if ( a[i] == b ) n++;
+		if ( a[i] == '\0' ) return n;
+		i++;
+	} while ( 1 );
+}
+
+
+static SymOpList *getpg_old_ua(const char *sym, size_t s)
+{
+	char ua;
+	char *pg_type;
+	SymOpList *pg;
+
+	if ( strncmp(sym+s, "ua", 2) == 0 ) {
+		ua = sym[s+2];
+	} else {
+		ERROR("Unrecognised point group '%s'\n", sym);
+		return NULL;
+	}
+
+	pg_type = strndup(sym, s-1);
+	if ( pg_type == NULL ) {
+		ERROR("Couldn't allocate string.\n");
+		return NULL;
+	}
+
+	pg = getpg_old(pg_type);
+	if ( pg == NULL ) {
+		ERROR("Unrecognised point group type '%s'\n",
+		      pg_type);
+		return NULL;
+	}
+	free(pg_type);
+
+	switch ( ua ) {
+
+		case 'a' :
+		transform_ops(pg, v(0,0,0,1),
+		                  v(0,1,0,0),
+		                  v(-1,0,0,0));
+		break;
+
+		case 'b' :
+		transform_ops(pg, v(1,0,0,0),
+		                  v(0,0,0,1),
+		                  v(0,-1,0,0));
+		break;
+
+		case 'c' :
+		/* No transformation needed */
+		break;
+
+		default :
+		ERROR("Bad unique axis '%c'\n", ua);
+		free_symoplist(pg);
+		return NULL;
+
+	}
+
+	return pg;
+}
+
+
+/**
+ * get_pointgroup:
+ * @sym: A string representation of a point group
+ *
+ * This function parses @sym and returns the corresponding %SymOpList.
+ * In the string representation of the point group, use a preceding minus sign
+ * for any character which would have a "bar".  Trigonal groups must be suffixed
+ * with either "_H" or "_R" for a hexagonal or rhombohedral lattice
+ * respectively.
+ *
+ * Examples: -1 1 2/m 2 m mmm 222 mm2 4/m 4 -4 4/mmm 422 -42m -4m2 4mm
+ * 3_R -3_R 32_R 3m_R -3m_R 3_H -3_H 321_H 312_H 3m1_H 31m_H -3m1_H -31m_H
+ * 6/m 6 -6 6/mmm 622 -62m -6m2 6mm 23 m-3 432 -43m m-3m.
+ **/
+SymOpList *get_pointgroup(const char *sym)
+{
+	int n_space, n_underscore;
+
+	n_space = char_count(sym, ' ');
+	n_underscore = char_count(sym, '_');
+
+	if ( n_space == 0 ) {
+
+		/* No spaces nor underscores -> old system */
+		if ( n_underscore == 0 ) return getpg_old(sym);
+
+		/* No spaces and 1 underscore -> old system + lattice or UA */
+		if ( n_underscore == 1 ) {
+
+			const char *s;
+
+			s = strchr(sym, '_');
+			assert(s != NULL);
+			s++;
+
+			/* Old system with H/R lattice? */
+			if ( (s[0] == 'H') || (s[0] == 'R') ) {
+				return getpg_old(sym);
+			}
+
+			/* Old system with unique axis */
+			return getpg_old_ua(sym, s-sym);
+
+		}
+
+	}
+
+	if ( n_space > 4 ) {
+		ERROR("Unrecognised point group '%s'\n", sym);
+		return NULL;
+	}
+
+	/* 1-4 spaces -> new system, possibly with UA and/or lattice */
+
+	/* FIXME: Implementation */
+
+	return NULL;
+}
+
+
 static void do_op(const struct sym_op *op,
                   signed int h, signed int k, signed int l,
                   signed int *he, signed int *ke, signed int *le)
diff --git a/tests/symmetry_check.c b/tests/symmetry_check.c
index edbacbf1..648814d0 100644
--- a/tests/symmetry_check.c
+++ b/tests/symmetry_check.c
@@ -152,6 +152,7 @@ static void check_pg_props(const char *pg, int answer, int centro, int *fail)
 
 	if ( sym == NULL ) {
 		*fail = 1;
+		STATUS("%15s : NULL!\n", pg);
 		return;
 	}
 
@@ -273,7 +274,6 @@ int main(int argc, char *argv[])
 	check_pg_props(  "-42m",  8,  0, &fail);
 	check_pg_props(  "-4m2",  8,  0, &fail);
 	check_pg_props( "4/mmm", 16,  1, &fail);
-	check_pg_props( "4/mmm_uaa", 16,  1, &fail);
 	STATUS("\n");
 
 	check_pg_props(  "3_R",   3,  0, &fail);
@@ -392,5 +392,23 @@ int main(int argc, char *argv[])
 	check_subgroup("-3_R",  "-1", 1, 1,  3, &fail);
 	check_subgroup("6",  "2", 1, 1,  3, &fail);
 
+	/* Check some weird settings */
+	STATUS("\nWeird settings:\n");
+	check_pg_props( "2_uaa", 2, 0, &fail);
+	check_pg_props( "2_uab", 2, 0, &fail);
+	check_pg_props( "2_uac", 2, 0, &fail);
+	check_pg_props( "4_uaa", 4, 0, &fail);
+	check_pg_props( "4_uab", 4, 0, &fail);
+	check_pg_props( "4_uac", 4, 0, &fail);
+	check_pg_props( "4/m_uaa", 8, 1, &fail);
+	check_pg_props( "4/m_uab", 8, 1, &fail);
+	check_pg_props( "4/m_uac", 8, 1, &fail);
+
+	/* Check "new style" parsing */
+	STATUS("\nNew style:\n");
+	check_pg_props( "2 1 1", 2, 0, &fail);
+	check_pg_props( "1 2 1", 2, 0, &fail);
+	check_pg_props( "1 1 2", 2, 0, &fail);
+
 	return fail;
 }