Work in progress...

2 files changed, 267 insertions, 615 deletions

diff --git a/src/symmetry.c b/src/symmetry.c
index 3a9d505e..77ac3ec3 100644
--- a/src/symmetry.c
+++ b/src/symmetry.c
@@ -23,12 +23,6 @@
 #include "utils.h"
 
 
-#ifdef DEBUG
-#define SYM_DEBUG STATUS
-#else /* DEBUG */
-#define SYM_DEBUG(...)
-#endif /* DEBUG */
-
 /**
  * SECTION:symmetry
  * @short_description: Point symmetry handling
@@ -42,11 +36,26 @@
  */
 
 
-struct sym_op {
-	signed int h;
-	signed int k;
-	signed int l;
-	int op;
+enum lattice_type
+{
+	L_TRICLINIC,
+	L_MONOCLINIC,
+	L_ORTHORHOMBIC,
+	L_TETRAGONAL,
+	L_RHOMBOHEDRAL,
+	L_TRIGONAL,
+	L_HEXAGONAL,
+	L_CUBIC,
+};
+
+
+struct sym_op
+{
+	signed int *h;
+	signed int *k;
+	signed int *l;  /* Contributions to h, k and l from h, k, i and l */
+	enum lattice_type latt;
+	char op[6];
 };
 
 
@@ -59,416 +68,341 @@ struct sym_op {
  * @include: "symmetry.h"
  * @Image:
  *
- * Wibble
+ * The SymOpList is an opaque data structure containing a list of point symmetry
+ * operations.  It could represent an point group or a list of indexing
+ * ambiguities (twin laws), or similar.
  */
 
-struct _symoplist {
-	struct sym_op *items;
+struct _symoplist
+{
+	struct sym_op *ops;
 	int n_ops;
 	int max_ops;
 };
 
 
 
-static void alloc_items(SymOpList *items)
+static void alloc_ops(SymOpList *ops)
 {
-	items->items = realloc(items->items,
-	                       items->max_ops*sizeof(struct sym_op));
+	ops->ops = realloc(ops->ops, ops->max_ops*sizeof(struct sym_op));
 }
 
 
-/**
- * new_items:
- *
- * Creates a new %SymOpList.
- *
- * Returns: The new list, or NULL.
- **/
-SymOpList *new_items()
+/* Creates a new SymOpList */
+static SymOpList *new_symoplist()
 {
 	SymOpList *new;
 	new = malloc(sizeof(SymOpList));
 	if ( new == NULL ) return NULL;
-	new->max_ops = 1024;
+	new->max_ops = 16;
 	new->n_ops = 0;
-	new->items = NULL;
-	alloc_items(new);
+	new->ops = NULL;
+	alloc_ops(new);
 	return new;
 }
 
 
-void delete_items(SymOpList *items)
+/**
+ * free_symoplist:
+ *
+ * Frees a %SymOpList and all associated resources.
+ **/
+void free_symoplist(SymOpList *ops)
 {
-	if ( items == NULL ) return;
-	if ( items->items != NULL ) free(items->items);
-	free(items);
+	int i;
+
+	if ( ops == NULL ) return;
+	for ( i=0; i<ops->n_ops; i++ ) {
+		free(ops->ops[i].h);
+		free(ops->ops[i].k);
+		free(ops->ops[i].l);
+	}
+	if ( ops->ops != NULL ) free(ops->ops);
+	free(ops);
 }
 
 
-void add_item_with_op(SymOpList *items, signed int h, signed int k,
-                      signed int l, int op)
+/* Add a operation to a SymOpList */
+static void add_symop(SymOpList *ops,
+                      signed int *h, signed int *k, signed int *l)
 {
-	if ( items->n_ops == items->max_ops ) {
-		items->max_ops += 1024;
-		alloc_items(items);
+	if ( ops->n_ops == ops->max_ops ) {
+		/* Pretty sure this never happens, but still... */
+		ops->max_ops += 16;
+		alloc_ops(ops);
 	}
 
-	items->items[items->n_ops].h = h;
-	items->items[items->n_ops].k = k;
-	items->items[items->n_ops].l = l;
-	items->items[items->n_ops].op = op;
-	items->n_ops++;
+	ops->ops[ops->n_ops].h = h;
+	ops->ops[ops->n_ops].k = k;
+	ops->ops[ops->n_ops].l = l;
+	ops->n_ops++;
 }
 
 
-void add_item(SymOpList *items, signed int h, signed int k, signed int l)
+int num_ops(const SymOpList *ops)
 {
-	add_item_with_op(items, h, k, l, 0);
+	return ops->n_ops;
 }
 
 
-int find_item(SymOpList *items,
-                     signed int h, signed int k, signed int l)
+static signed int *v(signed int h, signed int k, signed int i, signed int l)
 {
-	int i;
+	signed int *vec = malloc(4*sizeof(signed int));
+	vec[0] = h;  vec[1] = k;  vec[2] = i;  vec[3] = l;
+	return vec;
+}
 
-	for ( i=0; i<items->n_ops; i++ ) {
-		if ( items->items[i].h != h ) continue;
-		if ( items->items[i].k != k ) continue;
-		if ( items->items[i].l != l ) continue;
-		return 1;
-	}
-	return 0;
+
+/********************************* Triclinic **********************************/
+
+static SymOpList *make_1bar()
+{
+	SymOpList *new = new_symoplist();
+	add_symop(new, v(1,0,0,0), v(0,1,0,0), v(0,0,0,1));
+	add_symop(new, v(-1,0,0,0), v(0,-1,0,0), v(0,0,0,-1));
+	return new;
 }
 
 
-static int find_op(SymOpList *items, int op)
+static SymOpList *make_1()
 {
-	int i;
+	SymOpList *new = new_symoplist();
+	add_symop(new, v(1,0,0,0), v(0,1,0,0), v(0,0,0,1));
+	return new;
+}
 
-	for ( i=0; i<items->n_ops; i++ ) {
-		if ( items->items[i].op == op ) return 1;
-	}
-	return 0;
+
+/********************************* Monoclinic *********************************/
+
+static SymOpList *make_2m()
+{
+	return NULL;
 }
 
 
-struct sym_op *get_item(SymOpList *items, int i)
+static SymOpList *make_2()
 {
-	if ( i >= items->n_ops ) return NULL;
-	return &items->items[i];
+	return NULL;
 }
 
 
-int num_items(const SymOpList *items)
+static SymOpList *make_m()
 {
-	return items->n_ops;
+	return NULL;
 }
 
 
-void union_op_items(SymOpList *items, SymOpList *newi)
+/******************************** Orthorhombic ********************************/
+
+static SymOpList *make_mmm()
 {
-	int n, i;
+	return NULL;
+}
 
-	n = num_items(newi);
-	for ( i=0; i<n; i++ ) {
 
-		struct sym_op *r = get_item(newi, i);
-		if ( find_op(items, r->op) ) continue;
+static SymOpList *make_222()
+{
+	return NULL;
+}
 
-		add_item_with_op(items, r->h, r->k, r->l, r->op);
 
-	}
+static SymOpList *make_mm2()
+{
+	return NULL;
 }
 
 
-void union_ops(SymOpList *items, SymOpList *newi)
+/********************************* Tetragonal *********************************/
+
+static SymOpList *make_4m()
 {
-	int n, i;
+	return NULL;
+}
 
-	n = num_items(newi);
-	for ( i=0; i<n; i++ ) {
 
-		struct sym_op *r = get_item(newi, i);
-		if ( find_item(items, r->h, r->k, r->l) ) continue;
+static SymOpList *make_4()
+{
+	return NULL;
+}
 
-		add_item_with_op(items, r->h, r->k, r->l, r->op);
 
-	}
+static SymOpList *make_4bar()
+{
+	return NULL;
 }
 
 
-SymOpList *intersection_ops(SymOpList *i1, SymOpList *i2)
+static SymOpList *make_4mmm()
 {
-	int n, i;
-	SymOpList *res = new_items();
+	return NULL;
+}
 
-	n = num_items(i1);
-	for ( i=0; i<n; i++ ) {
 
-		struct sym_op *r = get_item(i1, i);
-		if ( find_item(i2, r->h, r->k, r->l) ) {
-			add_item_with_op(res, r->h, r->k, r->l, r->op);
-		}
+static SymOpList *make_422()
+{
+	return NULL;
+}
 
-	}
 
-	return res;
+static SymOpList *make_4bar2m()
+{
+	return NULL;
 }
 
 
-/* Check if a reflection is in the asymmetric unit cell */
-static int check_cond(signed int h, signed int k, signed int l, const char *sym)
+static SymOpList *make_4mm()
 {
-	/* Triclinic */
-	if ( strcmp(sym, "1") == 0 )
-		return ( 1 );
-	if ( strcmp(sym, "-1") == 0 )
-		return ( (l>0)
-		      || ( (l==0) && (k>0) )
-		      || ( (l==0) && (k==0) && (h>=0) ) );
+	return NULL;
+}
 
-	/* Orthorhombic */
-	if ( strcmp(sym, "mmm") == 0 )
-		return ( (h>=0) && (k>=0) && (l>=0) );
-	if ( strcmp(sym, "222") == 0 )
-		return ( (h>=0) && (k>0) )
-		    || ( (h>0) && (k>=0) )
-		    || ( (h>=0) && (k==0) && (l>=0) )
-		    || ( (h==0) && (k>=0) && (l>=0) );
 
-	/* Tetragonal */
-	if ( strcmp(sym, "4/mmm") == 0 )
-		return ( (((h>0) && (k>=0)) || ((h==0) && (k==0))) && (l>=0)
-		         && (h>=k) );  /* Like 6/mmm */
-	if ( strcmp(sym, "422") == 0 )
-		return ( (((h>0) && (k>=0)) || ((h==0) && (k==0)))
-		         && (h>=k) );
-	if ( strcmp(sym, "4/m") == 0 )
-		return ( (((h>0) && (k>=0)) || ((h==0) && (k==0))) && (l>=0) );
-	if ( strcmp(sym, "4") == 0 )
-		return ( ((h>0) && (k>=0)) || ((h==0) && (k==0)) );
+/******************************** Rhombohedral ********************************/
 
-	/* Hexagonal */
-	if ( strcmp(sym, "6") == 0 )
-		return ( ((h>0) && (k>=0)) || ((h==0) && (k==0)) );
-	if ( strcmp(sym, "6/m") == 0 )
-		return ( (((h>0) && (k>=0)) || ((h==0) && (k==0))) && (l>=0) );
-	if ( strcmp(sym, "6/mmm") == 0 )
-		return ( (((h>0) && (k>=0)) || ((h==0) && (k==0))) && (l>=0)
-		         && (h>=k) );
 
-	/* TODO: Add more groups here */
+/********************************** Hexgonal **********************************/
 
-	return 1;
+static SymOpList *make_6m()
+{
+	return NULL;
 }
 
 
-/* Macros for checking the above conditions and returning if satisfied */
-#define CHECK_COND(h, k, l, sym)                   \
-	if ( check_cond((h), (k), (l), (sym)) ) {  \
-		*hp = (h);  *kp = (k);  *lp = (l); \
-		return;                            \
-	}
+static SymOpList *make_6()
+{
+	return NULL;
+}
 
 
-int num_general_equivs(const char *sym)
+static SymOpList *make_6bar()
 {
-	/* Triclinic */
-	if ( strcmp(sym, "1") == 0 ) return 1;
-	if ( strcmp(sym, "-1") == 0 ) return 2;
+	return NULL;
+}
 
-	/* Orthorhombic */
-	if ( strcmp(sym, "222") == 0 ) return 4;
-	if ( strcmp(sym, "mmm") == 0 ) return 8;
 
-	/* Tetragonal */
-	if ( strcmp(sym, "4") == 0 ) return 4;
-	if ( strcmp(sym, "4/m") == 0 ) return 8;
-	if ( strcmp(sym, "422") == 0 ) return 8;
-	if ( strcmp(sym, "4/mmm") == 0 ) return 16;
+static SymOpList *make_6mmm()
+{
+	return NULL;
+}
 
-	/* Hexagonal */
-	if ( strcmp(sym, "6") == 0 ) return 6;
-	if ( strcmp(sym, "6/m") == 0 ) return 12;
-	if ( strcmp(sym, "6/mmm") == 0 ) return 24;
 
-	/* TODO: Add more groups here */
+static SymOpList *make_622()
+{
+	return NULL;
+}
+
 
-	return 1;
+static SymOpList *make_6bar2m()
+{
+	return NULL;
 }
 
 
-void get_general_equiv(signed int h, signed int k, signed int l,
-                       signed int *he, signed int *ke, signed int *le,
-                       const char *sym, int idx)
+static SymOpList *make_6mm()
 {
-	signed int i = -h-k;
+	return NULL;
+}
 
-	/* The returned indices when idx=0 *must* be the same as the input.
-	 * After that, the order does not matter. */
 
-	if ( strcmp(sym, "1") == 0 ) {
-		*he = h;   *ke = k;   *le = l;  return;
-	}
+/************************************ Cubic ***********************************/
 
-	if ( strcmp(sym, "-1") == 0 ) {
-		switch ( idx ) {
-		case 0 : *he = h;   *ke = k;   *le = l;   return;
-		case 1 : *he = -h;  *ke = -k;  *le = -l;  return;
-		}
-	}
 
-	if ( strcmp(sym, "222") == 0 ) {
-		switch ( idx ) {
-		case  0 : *he = h;   *ke = k;   *le = l;  return;
-		case  1 : *he = -h;  *ke = -k;  *le = l;  return;
-		case  2 : *he = -h;  *ke = k;  *le = -l;  return;
-		case  3 : *he = h;  *ke = -k;  *le = -l;  return;
-		}
-	}
 
-	if ( strcmp(sym, "mmm") == 0 ) {
-		switch ( idx ) {
-		case  0 : *he = h;  *ke = k;  *le = l;  return;
-		case  1 : *he = -h; *ke = -k; *le = l;  return;
-		case  2 : *he = -h; *ke = k;  *le = -l; return;
-		case  3 : *he = h;  *ke = -k; *le = -l; return;
-		case  4 : *he = -h; *ke = -k; *le = -l; return;
-		case  5 : *he = h;  *ke = k;  *le = -l; return;
-		case  6 : *he = h;  *ke = -k; *le = l;  return;
-		case  7 : *he = -h; *ke = k;  *le = l;  return;
-		}
-	}
+SymOpList *get_pointgroup(const char *sym)
+{
+	/* Triclinic */
+	if ( strcmp(sym, "-1") == 0 ) return make_1bar();
+	if ( strcmp(sym, "1") == 0 ) return make_1();
 
-	if ( strcmp(sym, "4") == 0 ) {
-		switch ( idx ) {
-		case  0 : *he = h;   *ke = k;   *le = l;  return;
-		case  1 : *he = -h;  *ke = -k;  *le = l;  return;
-		case  2 : *he = -k;  *ke = h;   *le = l;  return;
-		case  3 : *he = k;   *ke = -h;  *le = l;  return;
-		}
-	}
+	/* Monoclinic */
+	if ( strcmp(sym, "2/m") == 0 ) return make_2m();
+	if ( strcmp(sym, "2") == 0 ) return make_2();
+	if ( strcmp(sym, "m") == 0 ) return make_m();
 
-	if ( strcmp(sym, "4/m") == 0 ) {
-		switch ( idx ) {
-		case  0 : *he = h;   *ke = k;   *le = l;  return;
-		case  1 : *he = -h;  *ke = -k;  *le = l;  return;
-		case  2 : *he = -k;  *ke = h;   *le = l;  return;
-		case  3 : *he = k;   *ke = -h;  *le = l;  return;
-		case  4 : *he = -h;  *ke = -k;  *le = -l; return;
-		case  5 : *he = h;   *ke = k;   *le = -l; return;
-		case  6 : *he = k;   *ke = -h;  *le = -l; return;
-		case  7 : *he = -k;  *ke = h;   *le = -l; return;
-		}
-	}
+	/* Orthorhombic */
+	if ( strcmp(sym, "mmm") == 0 ) return make_mmm();
+	if ( strcmp(sym, "222") == 0 ) return make_222();
+	if ( strcmp(sym, "mm2") == 0 ) return make_mm2();
 
-	if ( strcmp(sym, "422") == 0 ) {
-		switch ( idx ) {
-		case  0 : *he = h;   *ke = k;   *le = l;  return;
-		case  1 : *he = -h;  *ke = -k;  *le = l;  return;
-		case  2 : *he = -k;  *ke = h;   *le = l;  return;
-		case  3 : *he = k;   *ke = -h;  *le = l;  return;
-		case  4 : *he = -h;  *ke = k;   *le = -l; return;
-		case  5 : *he = h;   *ke = -k;  *le = -l; return;
-		case  6 : *he = k;   *ke = h;   *le = -l; return;
-		case  7 : *he = -k;  *ke = -h;  *le = -l; return;
-		}
-	}
+	/* Tetragonal */
+	if ( strcmp(sym, "4/m") == 0 ) return make_4m();
+	if ( strcmp(sym, "4") == 0 ) return make_4();
+	if ( strcmp(sym, "4bar") == 0 ) return make_4bar();
+	if ( strcmp(sym, "4/mmm") == 0 ) return make_4mmm();
+	if ( strcmp(sym, "422") == 0 ) return make_422();
+	if ( strcmp(sym, "4bar2m") == 0 ) return make_4bar2m();
+	if ( strcmp(sym, "4mm") == 0 ) return make_4mm();
 
-	if ( strcmp(sym, "4/mmm") == 0 ) {
-		switch ( idx ) {
-		case  0 : *he = h;   *ke = k;   *le = l;  return;
-		case  1 : *he = -h;  *ke = -k;  *le = l;  return;
-		case  2 : *he = -k;  *ke = h;   *le = l;  return;
-		case  3 : *he = k;   *ke = -h;  *le = l;  return;
-		case  4 : *he = -h;  *ke = k;   *le = -l; return;
-		case  5 : *he = h;   *ke = -k;  *le = -l; return;
-		case  6 : *he = k;   *ke = h;   *le = -l; return;
-		case  7 : *he = -k;  *ke = -h;  *le = -l; return;
-		case  8 : *he = -h;  *ke = -k;  *le = -l; return;
-		case  9 : *he = h;   *ke = k;   *le = -l; return;
-		case 10 : *he = k;   *ke = -h;  *le = -l; return;
-		case 11 : *he = -k;  *ke = h;   *le = -l; return;
-		case 12 : *he = h;   *ke = -k;  *le = l;  return;
-		case 13 : *he = -h;  *ke = k;   *le = l;  return;
-		case 14 : *he = -k;  *ke = -h;  *le = l;  return;
-		case 15 : *he = k;   *ke = h;   *le = l;  return;
-		}
-	}
+	/* Hexagonal */
+	if ( strcmp(sym, "6/m") == 0 ) return make_6m();
+	if ( strcmp(sym, "6") == 0 ) return make_6();
+	if ( strcmp(sym, "6bar") == 0 ) return make_6bar();
+	if ( strcmp(sym, "6/mmm") == 0 ) return make_6mmm();
+	if ( strcmp(sym, "622") == 0 ) return make_622();
+	if ( strcmp(sym, "6bar2m") == 0 ) return make_6bar2m();
+	if ( strcmp(sym, "6mm") == 0 ) return make_6mm();
+}
 
-	if ( strcmp(sym, "6") == 0 ) {
-		switch ( idx ) {
-		case  0 : *he = h;   *ke = k;   *le = l;  return;
-		case  1 : *he = i;   *ke = h;   *le = l;  return;
-		case  2 : *he = k;   *ke = i;   *le = l;  return;
-		case  3 : *he = -h;  *ke = -k;  *le = l;  return;
-		case  4 : *he = -i;  *ke = -h;  *le = l;  return;
-		case  5 : *he = -k;  *ke = -i;  *le = l;  return;
-		}
-	}
 
-	if ( strcmp(sym, "6/m") == 0 ) {
-		switch ( idx ) {
-		case  0 : *he = h;   *ke = k;   *le = l;   return;
-		case  1 : *he = i;   *ke = h;   *le = l;   return;
-		case  2 : *he = k;   *ke = i;   *le = l;   return;
-		case  3 : *he = -h;  *ke = -k;  *le = l;   return;
-		case  4 : *he = -i;  *ke = -h;  *le = l;   return;
-		case  5 : *he = -k;  *ke = -i;  *le = l;   return;
-		case  6 : *he = h;   *ke = k;   *le = -l;  return;
-		case  7 : *he = i;   *ke = h;   *le = -l;  return;
-		case  8 : *he = k;   *ke = i;   *le = -l;  return;
-		case  9 : *he = -h;  *ke = -k;  *le = -l;  return;
-		case 10 : *he = -i;  *ke = -h;  *le = -l;  return;
-		case 11 : *he = -k;  *ke = -i;  *le = -l;  return;
-		}
-	}
+/**
+ * num_ops:
+ * @ops: A %SymOpList
+ *
+ * Returns: the number of operations in @ops.
+ **/
+int num_ops(SymOpList *ops)
+{
+	return ops->n_ops;
+}
 
-	if ( strcmp(sym, "6/mmm") == 0 ) {
-		switch ( idx ) {
-		case  0 : *he = h;   *ke = k;   *le = l;   return;
-		case  1 : *he = i;   *ke = h;   *le = l;   return;
-		case  2 : *he = k;   *ke = i;   *le = l;   return;
-		case  3 : *he = -h;  *ke = -k;  *le = l;   return;
-		case  4 : *he = -i;  *ke = -h;  *le = l;   return;
-		case  5 : *he = -k;  *ke = -i;  *le = l;   return;
-		case  6 : *he = k;   *ke = h;   *le = -l;  return;
-		case  7 : *he = h;   *ke = i;   *le = -l;  return;
-		case  8 : *he = i;   *ke = k;   *le = -l;  return;
-		case  9 : *he = -k;  *ke = -h;  *le = -l;  return;
-		case 10 : *he = -h;  *ke = -i;  *le = -l;  return;
-		case 11 : *he = -i;  *ke = -k;  *le = -l;  return;
-		case 12 : *he = -h;  *ke = -k;  *le = -l;  return;
-		case 13 : *he = -i;  *ke = -h;  *le = -l;  return;
-		case 14 : *he = -k;  *ke = -i;  *le = -l;  return;
-		case 15 : *he = h;   *ke = k;   *le = -l;  return;
-		case 16 : *he = i;   *ke = h;   *le = -l;  return;
-		case 17 : *he = k;   *ke = i;   *le = -l;  return;
-		case 18 : *he = -k;  *ke = -h;  *le = l;   return;
-		case 19 : *he = -h;  *ke = -i;  *le = l;   return;
-		case 20 : *he = -i;  *ke = -k;  *le = l;   return;
-		case 21 : *he = k;   *ke = h;   *le = l;   return;
-		case 22 : *he = h;   *ke = i;   *le = l;   return;
-		case 23 : *he = i;   *ke = k;   *le = l;   return;
-		}
-	}
 
-	/* TODO: Add more groups here */
+/**
+ * get_equiv:
+ * @ops: A %SymOpList
+ * @idx: Index of the operation to use
+ * @h: index of reflection
+ * @k: index of reflection
+ * @l: index of reflection
+ * @he: location to store h index of equivalent reflection
+ * @ke: location to store k index of equivalent reflection
+ * @le: location to store l index of equivalent reflection
+ *
+ * This function applies the @idx-th symmetry operation from @ops to the
+ * reflection @h, @k, @l, and stores the result at @he, @ke and @le.
+ *
+ * If you don't mind that the same equivalent might appear twice, simply call
+ * this function the number of times returned by num_ops(), using the actual
+ * point group.  If repeating the same equivalent twice (for example, if the
+ * given reflection is a special high-symmetry one), call special_position()
+ * first to get a "specialised" SymOpList and use that instead.
+ **/
+void get_equiv(SymOpList *ops, int idx,
+               signed int h, signed int k, signed int l,
+               signed int *he, signed int *ke, signed int *le)
+{
+	signed int i = -h-k;
+	struct sym_op op = ops->ops[idx];
 
-	ERROR("Unrecognised symmetry '%s'\n", sym);
-	abort();
+	*he = h*op.h[0] + k*op.h[1] + i*op.h[2] + l*op.h[3];
+	*ke = h*op.k[0] + k*op.h[1] + i*op.k[2] + l*op.k[3];
+	*le = h*op.l[0] + k*op.h[1] + i*op.l[2] + l*op.l[3];
 }
 
 
-/* Given a reflection and a point group, this returns (by reference) the indices
- * of the "idx"th equivalent reflection, taking special positions into account.
- * It returns "idx" if successful.  Otherwise, it returns the number of
- * equivalents for the particular reflection (taking special positions into
- * account).  Therefore, set idx=-1 to get the number of equivalents. */
-static int special_position(signed int hs, signed int ks, signed int ls,
-                            signed int *hp, signed int *kp, signed int *lp,
-                            const char *sym, signed int idx)
+/**
+ * special_position:
+ * @ops: A %SymOpList, usually corresponding to a point group
+ * @h: index of a reflection
+ * @k: index of a reflection
+ * @l: index of a reflection
+ *
+ * This function determines which operations in @ops map the reflection @h, @k,
+ * @l onto itself, and returns a new %SymOpList containing only the operations
+ * from @ops which do not do so.
+ *
+ * Returns: the "specialised" %SymOpList.
+ **/
+SymOpList *special_position(SymOpList *ops,
+                            signed int h, signed int k, signed int l)
 {
 	int n_general;
 	int i;
@@ -481,15 +415,15 @@ static int special_position(signed int hs, signed int ks, signed int ls,
 		return 0;
 	}
 
-	equivs = new_items();
-	n_general = num_general_equivs(sym);
+	equivs = new_symoplist();
+	n_general = num_general_equivs(ops);
 
 	for ( i=0; i<n_general; i++ ) {
 
 		signed int h, k, l;
 
-		/* Get equivalent according to the holohedral group */
-		get_general_equiv(hs, ks, ls, &h, &k, &l, sym, i);
+		/* Get equivalent according to the point group */
+		get_general_equiv(ops, i, hs, ks, ls, &h, &k, &l);
 
 		/* Already got this one? */
 		if ( find_item(equivs, h, k, l) ) continue;
@@ -498,7 +432,7 @@ static int special_position(signed int hs, signed int ks, signed int ls,
 			*hp = h;
 			*kp = k;
 			*lp = l;
-			delete_items(equivs);
+			delete_ops(equivs);
 			return n_equivs;
 		}
 		add_item(equivs, h, k, l);
@@ -506,329 +440,47 @@ static int special_position(signed int hs, signed int ks, signed int ls,
 
 	}
 
-	delete_items(equivs);
+	delete_ops(equivs);
 	return n_equivs;
 }
 
 
-void get_equiv(signed int h, signed int k, signed int l,
-               signed int *he, signed int *ke, signed int *le,
-               const char *sym, int idx)
-{
-	special_position(h, k, l, he, ke, le, sym, idx);
-}
-
-
-int num_equivs(signed int h, signed int k, signed int l, const char *sym)
-{
-	return special_position(h, k, l, NULL, NULL, NULL, sym, -1);
-}
-
-
-void get_asymm(signed int h, signed int k, signed int l,
-               signed int *hp, signed int *kp, signed int *lp,
-               const char *sym)
+void get_asymm(SymOpList *ops, int idx,
+               signed int h, signed int k, signed int l,
+               signed int *hp, signed int *kp, signed int *lp)
 {
 	int nequiv = num_equivs(h, k, l, sym);
 	int p;
 
-	SYM_DEBUG("------ %i %i %i\n", h, k, l);
 	for ( p=0; p<nequiv; p++ ) {
 		signed int he, ke, le;
 		get_equiv(h, k, l, &he, &ke, &le, sym, p);
-		SYM_DEBUG("%i : %i %i %i\n", p, he, ke, le);
 		CHECK_COND(he, ke, le, sym);
 	}
-
-	/* Should never reach here */
-	ERROR("No match found in %s for %i %i %i\n", sym, h, k, l);
-	abort();
 }
 
 
-/* This is kind of like a "numerical" left coset decomposition.
- * Given a reflection index and a point group, it returns the list of twinning
- * possibilities.
+/**
+ * get_twins:
  *
- * To count the number of possibilities, use num_items() on the result.
+ * Calculate twinning laws.
+ *
+ * To count the number of possibilities, use num_ops() on the result.
  */
-static SymOpList *coset_decomp(signed int hs, signed int ks, signed int ls,
-                               const char *holo, const char *mero)
+SymOpList *get_twins(SymOpList *source, SymOpList *target)
 {
-	int n_mero, n_holo;
+	int n_src, n_tgt;
 	int i;
 	signed int h, k, l;
-	SymOpList *twins = new_items();
-
-	/* Start by putting the given reflection into the asymmetric cell
-	 * for its (probably merohedral) point group. */
-	get_asymm(hs, ks, ls, &h, &k, &l, mero);
-
-	/* How many equivalents in the holohedral point group are not
-	 * equivalents according to the (possibly) merohedral group? */
-	n_holo = num_general_equivs(holo);
-	n_mero = num_general_equivs(mero);
-
-	for ( i=0; i<n_holo; i++ ) {
+	SymOpList *twins = new_symoplist();
 
-		signed int h_holo, k_holo, l_holo;
-		signed int hs_holo, ks_holo, ls_holo;
+	n_src = num_ops(source);
+	n_tgt = num_ops(target);
 
-		/* Get equivalent according to the holohedral group */
-		get_general_equiv(h, k, l, &hs_holo, &ks_holo, &ls_holo,
-		                  holo, i);
+	for ( i=0; i<n_src; i++ ) {
 
-		/* Put it into the asymmetric cell for the merohedral group */
-		get_asymm(hs_holo, ks_holo, ls_holo,
-		          &h_holo, &k_holo, &l_holo, mero);
-
-		/* Already got this one?
-		 * Note: The list "twins" starts empty, so the first iteration
-		 * (i=0) will add the original reflection to the list along with
-		 * the identity operation. */
-		if ( find_item(twins, h_holo, k_holo, l_holo) ) continue;
-
-		add_item_with_op(twins, h_holo, k_holo, l_holo, i);
 
 	}
 
 	return twins;
 }
-
-
-/* Work out the twinning possibilities for this pattern.
- * To use the result, call get_general_equiv() on each reflection using
- * the holohedral point group (use get_holohedral() for this), and for "idx"
- * give each "op" field from the list returned by this function. */
-SymOpList *get_twins(const char *holo, const char *mero)
-{
-	int i;
-	SymOpList *ops = new_items();
-	int expected, actual;
-	SymOpList *items;
-
-	/* Run the coset decomposition for every reflection in the "pattern",
-	 * and see which gives the highest number of possibilities.  This
-	 * approach takes into account that a pattern consisting entirely of
-	 * special reflections might have fewer twin possibilities. */
-	for ( i=0; i<num_items(items); i++ ) {
-
-		signed int h, k, l;
-		struct sym_op *item;
-		SymOpList *new_ops;
-
-		item = get_item(items, i);
-
-		h = item->h;
-		k = item->k;
-		l = item->l;
-
-		new_ops = coset_decomp(h, k, l, holo, mero);
-		union_op_items(ops, new_ops);
-		delete_items(new_ops);
-
-	}
-
-	/* Idiot check */
-	actual = num_items(ops);
-	expected = num_general_equivs(holo) / num_general_equivs(mero);
-	if ( actual != expected ) {
-		ERROR("Found %i twin possibilities, but expected %i.\n",
-		       actual, expected);
-		ERROR("I couldn't find the number of twin laws that I expected."
-		      " This is an internal error, and shouldn't happen. "
-		      "Sorry.\n");
-		abort();
-	}
-
-	return ops;
-}
-
-
-int find_unique_equiv(SymOpList *items, signed int h, signed int k,
-                      signed int l, const char *mero, signed int *hu,
-                      signed int *ku, signed int *lu)
-{
-	int i;
-	int found = 0;
-
-	for ( i=0; i<num_equivs(h, k, l, mero); i++ ) {
-
-		signed int he, ke, le;
-		int f;
-		get_equiv(h, k, l, &he, &ke, &le, mero, i);
-		f = find_item(items, he, ke, le);
-
-		/* There must only be one equivalent.  If there are more, it
-		 * indicates that the user lied about the input symmetry.
-		 * This situation should have been checked for earlier by
-		 * calling check_symmetry() with 'items' and 'mero'. */
-
-		if ( f && !found ) {
-			*hu = he;  *ku = ke;  *lu = le;
-			return 1;
-		}
-
-	}
-
-	return 0;
-}
-
-
-/* Returns true if the point group is 23, m-3, 432, -43m or m-3m
- * (i.e. T, Th, O, Td or Oh) */
-int is_polyhedral(const char *sym)
-{
-	/* Triclinic */
-	if ( strcmp(sym, "1") == 0 ) return 0;
-	if ( strcmp(sym, "-1") == 0 ) return 0;
-
-	/* Orthorhombic */
-	if ( strcmp(sym, "222") == 0 ) return 0;
-	if ( strcmp(sym, "mmm") == 0 ) return 0;
-
-	/* Tetragonal */
-	if ( strcmp(sym, "4") == 0 ) return 0;
-	if ( strcmp(sym, "4/m") == 0 ) return 0;
-	if ( strcmp(sym, "422") == 0 ) return 0;
-	if ( strcmp(sym, "4/mmm") == 0 ) return 0;
-
-	/* Hexagonal */
-	if ( strcmp(sym, "6") == 0 ) return 0;
-	if ( strcmp(sym, "6/m") == 0 ) return 0;
-	if ( strcmp(sym, "6/mmm") == 0 ) return 0;
-
-	/* TODO: Add more groups here */
-
-	ERROR("Don't know if '%s' is polyhedral or not.\n", sym);
-	abort();
-}
-
-
-/* Returns the order of the characteristic axis of proper or improper rotation
- * for the point group.  This should be the rotation about the c axis. */
-int rotational_order(const char *sym)
-{
-	/* Triclinic */
-	if ( strcmp(sym, "1") == 0 ) return 1;
-	if ( strcmp(sym, "-1") == 0 ) return 2;
-
-	/* Orthorhombic */
-	if ( strcmp(sym, "222") == 0 ) return 2;
-	if ( strcmp(sym, "mmm") == 0 ) return 2;
-
-	/* Tetragonal */
-	if ( strcmp(sym, "4") == 0 ) return 4;
-	if ( strcmp(sym, "4/m") == 0 ) return 4;
-	if ( strcmp(sym, "422") == 0 ) return 4;
-	if ( strcmp(sym, "4/mmm") == 0 ) return 4;
-
-	/* Hexagonal */
-	if ( strcmp(sym, "6") == 0 ) return 6;
-	if ( strcmp(sym, "6/m") == 0 ) return 6;
-	if ( strcmp(sym, "6/mmm") == 0 ) return 6;
-
-	/* TODO: Add more groups here */
-
-	ERROR("Couldn't find rotational order for '%s'.\n", sym);
-	abort();
-}
-
-
-int has_perpendicular_mirror(const char *sym)
-{
-	/* Triclinic */
-	if ( strcmp(sym, "1") == 0 ) return 0;
-	if ( strcmp(sym, "-1") == 0 ) return 0;
-
-	/* Orthorhombic */
-	if ( strcmp(sym, "222") == 0 ) return 0;
-	if ( strcmp(sym, "mmm") == 0 ) return 1;
-
-	/* Tetragonal */
-	if ( strcmp(sym, "4") == 0 ) return 0;
-	if ( strcmp(sym, "4/m") == 0 ) return 1;
-	if ( strcmp(sym, "422") == 0 ) return 0;
-	if ( strcmp(sym, "4/mmm") == 0 ) return 1;
-
-	/* Hexagonal */
-	if ( strcmp(sym, "6") == 0 ) return 0;
-	if ( strcmp(sym, "6/m") == 0 ) return 1;
-	if ( strcmp(sym, "6/mmm") == 0 ) return 1;
-
-	/* TODO: Add more groups here */
-
-	ERROR("Couldn't find mirror definition for '%s'.\n", sym);
-	abort();
-}
-
-
-int has_bisecting_mirror_or_diad(const char *sym)
-{
-	/* Triclinic */
-	if ( strcmp(sym, "1") == 0 ) return 0;
-	if ( strcmp(sym, "-1") == 0 ) return 0;
-
-	/* Orthorhombic */
-	if ( strcmp(sym, "222") == 0 ) return 1;
-	if ( strcmp(sym, "mmm") == 0 ) return 1;
-
-	/* Tetragonal */
-	if ( strcmp(sym, "4") == 0 ) return 0;
-	if ( strcmp(sym, "4/m") == 0 ) return 0;
-	if ( strcmp(sym, "422") == 0 ) return 0;
-	if ( strcmp(sym, "4/mmm") == 0 ) return 1;
-
-	/* Hexagonal */
-	if ( strcmp(sym, "6") == 0 ) return 0;
-	if ( strcmp(sym, "6/m") == 0 ) return 1;
-	if ( strcmp(sym, "6/mmm") == 0 ) return 1;
-
-	/* TODO: Add more groups here */
-
-	ERROR("Couldn't find mirror definition for '%s'.\n", sym);
-	abort();
-}
-
-
-int check_symmetry(SymOpList *items, const char *sym)
-{
-	int i;
-	unsigned char *flags;
-
-	flags = new_list_flag();
-	for ( i=0; i<num_items(items); i++ ) {
-		struct sym_op *it = get_item(items, i);
-		set_flag(flags, it->h, it->k, it->l, 1);
-	}
-
-	for ( i=0; i<num_items(items); i++ ) {
-
-		int j;
-		struct sym_op *it = get_item(items, i);
-		int found = 0;
-
-		for ( j=0; j<num_equivs(it->h, it->k, it->l, sym); j++ ) {
-
-			signed int he, ke, le;
-			get_equiv(it->h, it->k, it->l, &he, &ke, &le, sym, j);
-
-			if ( abs(he) > INDMAX ) continue;
-			if ( abs(le) > INDMAX ) continue;
-			if ( abs(ke) > INDMAX ) continue;
-
-			found += lookup_flag(flags, he, ke, le);
-
-		}
-
-		if ( found > 1 ) {
-			free(flags);
-			return 1;  /* Symmetry is wrong! */
-		}
-
-	}
-
-	free(flags);
-
-	return 0;
-}
diff --git a/src/symmetry.h b/src/symmetry.h
index fbb80d26..23c00160 100644
--- a/src/symmetry.h
+++ b/src/symmetry.h
@@ -24,6 +24,12 @@
  **/
 typedef struct _symoplist SymOpList;
 
+extern void free_symoplist(SymOpList *ops);
+
+extern SymOpList *get_pointgroup(const char *sym);
+
+extern SymOpList *special_position(SymOpList *ops,
+                                   signed int h, signed int k, signed int l);
 
 extern void get_asymm(signed int h, signed int k, signed int l,
                       signed int *hp, signed int *kp, signed int *lp,
@@ -42,12 +48,6 @@ extern void get_general_equiv(signed int h, signed int k, signed int l,
                               signed int *he, signed int *ke, signed int *le,
                               const char *sym, int idx);
 
-extern SymOpList *get_twins(const char *holo, const char *mero);
-
-/* Properties of point groups */
-extern int is_polyhedral(const char *sym);
-extern int rotational_order(const char *sym);
-extern int has_perpendicular_mirror(const char *sym);
-extern int has_bisecting_mirror_or_diad(const char *sym);
+extern SymOpList *get_twins(SymOpList *source, SymOpList *target);
 
 #endif	/* SYMMETRY_H */