Symmetry stuff

11 files changed, 382 insertions, 208 deletions

diff --git a/Makefile.am b/Makefile.am
index a2fdb3c0..f9615e37 100644
--- a/Makefile.am
+++ b/Makefile.am
@@ -9,11 +9,12 @@ bin_PROGRAMS = src/pattern_sim src/process_hkl src/get_hkl src/indexamajig \
                contrib/alter_stream
 
 noinst_PROGRAMS = tests/list_check tests/integration_check \
-                  tests/pr_gradient_check
+                  tests/pr_gradient_check tests/symmetry_check
 
 TESTS = tests/list_check tests/first_merge_check tests/second_merge_check \
         tests/third_merge_check tests/fourth_merge_check \
-        tests/integration_check tests/pr_gradient_check
+        tests/integration_check tests/pr_gradient_check \
+        tests/symmetry_check
 
 EXTRA_DIST += tests/first_merge_check tests/second_merge_check \
               tests/third_merge_check tests/fourth_merge_check \
@@ -146,6 +147,10 @@ tests_integration_check_SOURCES = tests/integration_check.c src/peaks.c \
                                   src/thread-pool.c src/utils.c src/image.c \
                                   src/hdf5-file.c
 
+tests_symmetry_check_SOURCES = tests/symmetry_check.c src/symmetry.c \
+                               src/utils.c src/thread-pool.c
+
+
 tests_pr_gradient_check_SOURCES = tests/pr_gradient_check.c src/detector.c \
                                   src/cell.c src/geometry.c src/reflist.c \
                                   src/thread-pool.c src/utils.c src/peaks.c \
diff --git a/Makefile.in b/Makefile.in
index df5dd25d..b94918e7 100644
--- a/Makefile.in
+++ b/Makefile.in
@@ -44,11 +44,13 @@ bin_PROGRAMS = src/pattern_sim$(EXEEXT) src/process_hkl$(EXEEXT) \
 	contrib/alter_stream$(EXEEXT) $(am__EXEEXT_1) $(am__EXEEXT_2)
 noinst_PROGRAMS = tests/list_check$(EXEEXT) \
 	tests/integration_check$(EXEEXT) \
-	tests/pr_gradient_check$(EXEEXT) $(am__EXEEXT_3)
+	tests/pr_gradient_check$(EXEEXT) tests/symmetry_check$(EXEEXT) \
+	$(am__EXEEXT_3)
 TESTS = tests/list_check$(EXEEXT) tests/first_merge_check \
 	tests/second_merge_check tests/third_merge_check \
 	tests/fourth_merge_check tests/integration_check$(EXEEXT) \
-	tests/pr_gradient_check$(EXEEXT) $(am__EXEEXT_3)
+	tests/pr_gradient_check$(EXEEXT) tests/symmetry_check$(EXEEXT) \
+	$(am__EXEEXT_3)
 @BUILD_HDFSEE_TRUE@am__append_1 = src/hdfsee
 @BUILD_CUBEIT_TRUE@am__append_2 = src/cubeit
 @HAVE_OPENCL_TRUE@am__append_3 = src/diffraction-gpu.c src/cl-utils.c
@@ -306,6 +308,12 @@ tests_pr_gradient_check_OBJECTS =  \
 	$(am_tests_pr_gradient_check_OBJECTS)
 tests_pr_gradient_check_LDADD = $(LDADD)
 tests_pr_gradient_check_DEPENDENCIES = $(top_builddir)/lib/libgnu.a
+am_tests_symmetry_check_OBJECTS = tests/symmetry_check.$(OBJEXT) \
+	src/symmetry.$(OBJEXT) src/utils.$(OBJEXT) \
+	src/thread-pool.$(OBJEXT)
+tests_symmetry_check_OBJECTS = $(am_tests_symmetry_check_OBJECTS)
+tests_symmetry_check_LDADD = $(LDADD)
+tests_symmetry_check_DEPENDENCIES = $(top_builddir)/lib/libgnu.a
 DEFAULT_INCLUDES = -I.@am__isrc@
 depcomp = $(SHELL) $(top_srcdir)/build-aux/depcomp
 am__depfiles_maybe = depfiles
@@ -336,7 +344,8 @@ SOURCES = $(contrib_alter_stream_SOURCES) \
 	$(src_render_hkl_SOURCES) $(src_sum_stack_SOURCES) \
 	$(tests_gpu_sim_check_SOURCES) \
 	$(tests_integration_check_SOURCES) $(tests_list_check_SOURCES) \
-	$(tests_pr_gradient_check_SOURCES)
+	$(tests_pr_gradient_check_SOURCES) \
+	$(tests_symmetry_check_SOURCES)
 DIST_SOURCES = $(contrib_alter_stream_SOURCES) \
 	$(src_calibrate_detector_SOURCES) $(src_check_hkl_SOURCES) \
 	$(src_compare_hkl_SOURCES) $(am__src_cubeit_SOURCES_DIST) \
@@ -348,7 +357,8 @@ DIST_SOURCES = $(contrib_alter_stream_SOURCES) \
 	$(src_sum_stack_SOURCES) \
 	$(am__tests_gpu_sim_check_SOURCES_DIST) \
 	$(tests_integration_check_SOURCES) $(tests_list_check_SOURCES) \
-	$(tests_pr_gradient_check_SOURCES)
+	$(tests_pr_gradient_check_SOURCES) \
+	$(tests_symmetry_check_SOURCES)
 RECURSIVE_TARGETS = all-recursive check-recursive dvi-recursive \
 	html-recursive info-recursive install-data-recursive \
 	install-dvi-recursive install-exec-recursive \
@@ -751,6 +761,9 @@ tests_integration_check_SOURCES = tests/integration_check.c src/peaks.c \
                                   src/thread-pool.c src/utils.c src/image.c \
                                   src/hdf5-file.c
 
+tests_symmetry_check_SOURCES = tests/symmetry_check.c src/symmetry.c \
+                               src/utils.c src/thread-pool.c
+
 tests_pr_gradient_check_SOURCES = tests/pr_gradient_check.c src/detector.c \
                                   src/cell.c src/geometry.c src/reflist.c \
                                   src/thread-pool.c src/utils.c src/peaks.c \
@@ -1040,6 +1053,11 @@ tests/pr_gradient_check.$(OBJEXT): tests/$(am__dirstamp) \
 tests/pr_gradient_check$(EXEEXT): $(tests_pr_gradient_check_OBJECTS) $(tests_pr_gradient_check_DEPENDENCIES) tests/$(am__dirstamp)
 	@rm -f tests/pr_gradient_check$(EXEEXT)
 	$(AM_V_CCLD)$(LINK) $(tests_pr_gradient_check_OBJECTS) $(tests_pr_gradient_check_LDADD) $(LIBS)
+tests/symmetry_check.$(OBJEXT): tests/$(am__dirstamp) \
+	tests/$(DEPDIR)/$(am__dirstamp)
+tests/symmetry_check$(EXEEXT): $(tests_symmetry_check_OBJECTS) $(tests_symmetry_check_DEPENDENCIES) tests/$(am__dirstamp)
+	@rm -f tests/symmetry_check$(EXEEXT)
+	$(AM_V_CCLD)$(LINK) $(tests_symmetry_check_OBJECTS) $(tests_symmetry_check_LDADD) $(LIBS)
 
 mostlyclean-compile:
 	-rm -f *.$(OBJEXT)
@@ -1089,6 +1107,7 @@ mostlyclean-compile:
 	-rm -f tests/integration_check.$(OBJEXT)
 	-rm -f tests/list_check.$(OBJEXT)
 	-rm -f tests/pr_gradient_check.$(OBJEXT)
+	-rm -f tests/symmetry_check.$(OBJEXT)
 
 distclean-compile:
 	-rm -f *.tab.c
@@ -1139,6 +1158,7 @@ distclean-compile:
 @AMDEP_TRUE@@am__include@ @am__quote@tests/$(DEPDIR)/integration_check.Po@am__quote@
 @AMDEP_TRUE@@am__include@ @am__quote@tests/$(DEPDIR)/list_check.Po@am__quote@
 @AMDEP_TRUE@@am__include@ @am__quote@tests/$(DEPDIR)/pr_gradient_check.Po@am__quote@
+@AMDEP_TRUE@@am__include@ @am__quote@tests/$(DEPDIR)/symmetry_check.Po@am__quote@
 
 .c.o:
 @am__fastdepCC_TRUE@	$(AM_V_CC)depbase=`echo $@ | sed 's|[^/]*$$|$(DEPDIR)/&|;s|\.o$$||'`;\
diff --git a/src/diffraction.c b/src/diffraction.c
index aa04fd52..4133cd0b 100644
--- a/src/diffraction.c
+++ b/src/diffraction.c
@@ -100,14 +100,14 @@ static double sym_lookup_intensity(const double *intensities,
 	int i;
 	double ret = 0.0;
 
-	for ( i=0; i<num_equivs(sym); i++ ) {
+	for ( i=0; i<num_equivs(sym, NULL); i++ ) {
 
 		signed int he;
 		signed int ke;
 		signed int le;
 		double f, val;
 
-		get_equiv(sym, i, h, k, l, &he, &ke, &le);
+		get_equiv(sym, NULL, i, h, k, l, &he, &ke, &le);
 
 		f = (double)lookup_flag(flags, he, ke, le);
 		val = lookup_intensity(intensities, he, ke, le);
@@ -127,14 +127,14 @@ static double sym_lookup_phase(const double *phases,
 	int i;
 	double ret = 0.0;
 
-	for ( i=0; i<num_equivs(sym); i++ ) {
+	for ( i=0; i<num_equivs(sym, NULL); i++ ) {
 
 		signed int he;
 		signed int ke;
 		signed int le;
 		double f, val;
 
-		get_equiv(sym, i, h, k, l, &he, &ke, &le);
+		get_equiv(sym, NULL, i, h, k, l, &he, &ke, &le);
 
 		f = (double)lookup_flag(flags, he, ke, le);
 		val = lookup_phase(phases, he, ke, le);
diff --git a/src/get_hkl.c b/src/get_hkl.c
index 72e1c2a3..5edfa088 100644
--- a/src/get_hkl.c
+++ b/src/get_hkl.c
@@ -146,15 +146,17 @@ static RefList *twin_reflections(RefList *in,
 	Reflection *refl;
 	RefListIterator *iter;
 	RefList *out;
+	SymOpMask *m;
 
 	/* FIXME: Check properly by coset decomposition */
-	if ( num_equivs(holo) < num_equivs(mero) ) {
+	if ( num_equivs(holo, NULL) < num_equivs(mero, NULL) ) {
 		ERROR("%s is not a subgroup of %s!\n", symmetry_name(mero),
 		                                       symmetry_name(holo));
 		return NULL;
 	}
 
 	out = reflist_new();
+	m = new_symopmask(holo);
 
 	for ( refl = first_refl(in, &iter);
 	      refl != NULL;
@@ -178,13 +180,15 @@ static RefList *twin_reflections(RefList *in,
 		total = 0.0;
 		sigma = 0.0;
 		skip = 0;
-		n = num_equivs(holo);
+		special_position(holo, m, h, k, l);
+		n = num_equivs(holo, m);
+
 		for ( j=0; j<n; j++ ) {
 
 			signed int he, ke, le;
 			signed int hu, ku, lu;
 
-			get_equiv(holo, j, h, k, l, &he, &ke, &le);
+			get_equiv(holo, m, j, h, k, l, &he, &ke, &le);
 
 			/* Do we have this reflection?
 			 * We might not have the particular (merohedral)
@@ -229,15 +233,17 @@ static RefList *expand_reflections(RefList *in, const SymOpList *target,
 	Reflection *refl;
 	RefListIterator *iter;
 	RefList *out;
+	SymOpMask *m;
 
 	/* FIXME: Check properly */
-	if ( num_equivs(target) > num_equivs(initial) ) {
+	if ( num_equivs(target, NULL) > num_equivs(initial, NULL) ) {
 		ERROR("%s is not a subgroup of %s!\n", symmetry_name(initial),
 		                                       symmetry_name(target));
 		return NULL;
 	}
 
 	out = reflist_new();
+	m = new_symopmask(initial);
 
 	for ( refl = first_refl(in, &iter);
 	      refl != NULL;
@@ -250,7 +256,8 @@ static RefList *expand_reflections(RefList *in, const SymOpList *target,
 		get_indices(refl, &h, &k, &l);
 		intensity = get_intensity(refl);
 
-		n = num_equivs(initial);
+		special_position(initial, m, h, k, l);
+		n = num_equivs(initial, m);
 
 		/* For each equivalent in the higher symmetry group */
 		for ( j=0; j<n; j++ ) {
@@ -259,7 +266,7 @@ static RefList *expand_reflections(RefList *in, const SymOpList *target,
 			Reflection *new;
 
 			/* Get the equivalent */
-			get_equiv(initial, j, h, k, l, &he, &ke, &le);
+			get_equiv(initial, m, j, h, k, l, &he, &ke, &le);
 
 			/* Put it into the asymmetric unit for the target */
 			get_asymm(target, he, ke, le, &he, &ke, &le);
@@ -272,6 +279,8 @@ static RefList *expand_reflections(RefList *in, const SymOpList *target,
 
 	}
 
+	free_symopmask(m);
+
 	return out;
 }
 
@@ -460,6 +469,9 @@ int main(int argc, char *argv[])
 
 		Reflection *refl;
 		RefListIterator *iter;
+		SymOpMask *m;
+
+		m = new_symopmask(mero);
 
 		for ( refl = first_refl(input, &iter);
 		      refl != NULL;
@@ -471,10 +483,13 @@ int main(int argc, char *argv[])
 			get_indices(refl, &h, &k, &l);
 			inty = get_intensity(refl);
 
-			inty *= (double)num_equivs(mero);
+			special_position(mero, m, h, k, l);
+			inty *= (double)num_equivs(mero, m);
 			set_int(refl, inty);
 
 		}
+
+		free_symopmask(m);
 	}
 
 	if ( template ) {
diff --git a/src/povray.c b/src/povray.c
index 9fb267f4..d20c18e5 100644
--- a/src/povray.c
+++ b/src/povray.c
@@ -42,6 +42,7 @@ int povray_render_animation(UnitCell *cell, RefList *list, unsigned int nproc,
 	int i;
 	Reflection *refl;
 	RefListIterator *iter;
+	SymOpMask *m;
 
 	if ( (nproc > MAX_PROC) || (nproc < 1) ) {
 		ERROR("Number of processes must be a number between 1 and %i\n",
@@ -167,6 +168,8 @@ int povray_render_animation(UnitCell *cell, RefList *list, unsigned int nproc,
 
 	fprintf(fh, "}\n");
 
+	m = new_symopmask(sym);
+
 	max = 0.0;
 	for ( refl = first_refl(list, &iter);
 	      refl != NULL;
@@ -174,10 +177,9 @@ int povray_render_animation(UnitCell *cell, RefList *list, unsigned int nproc,
 
 		float val;
 		signed int h, k, l;
-		SymOpList *sp;
 
 		get_indices(refl, &h, &k, &l);
-		sp = special_position(sym, h, k, l);
+		special_position(sym, m, h, k, l);
 
 		switch ( wght ) {
 		case WGHT_I :
@@ -189,7 +191,7 @@ int povray_render_animation(UnitCell *cell, RefList *list, unsigned int nproc,
 			break;
 		case WGHT_COUNTS :
 			val = get_redundancy(refl);
-			val /= (double)num_equivs(sp);
+			val /= (double)num_equivs(sym, m);
 			break;
 		case WGHT_RAWCOUNTS :
 			val = get_redundancy(refl);
@@ -201,8 +203,6 @@ int povray_render_animation(UnitCell *cell, RefList *list, unsigned int nproc,
 
 		if ( val > max ) max = val;
 
-		free_symoplist(sp);
-
 	}
 	max /= boost;
 
@@ -211,7 +211,6 @@ int povray_render_animation(UnitCell *cell, RefList *list, unsigned int nproc,
 		max = scale_top;
 	}
 
-
 	for ( refl = first_refl(list, &iter);
 	      refl != NULL;
 	      refl = next_refl(refl, iter) ) {
@@ -220,12 +219,11 @@ int povray_render_animation(UnitCell *cell, RefList *list, unsigned int nproc,
 		int s;
 		float val, p, r, g, b, trans;
 		int j;
-		SymOpList *sp;
 		int neq;
 
 		get_indices(refl, &h, &k, &l);
-		sp = special_position(sym, h, k, l);
-		neq = num_equivs(sp);
+		special_position(sym, m, h, k, l);
+		neq = num_equivs(sym, m);
 
 		switch ( wght ) {
 		case WGHT_I :
@@ -294,7 +292,7 @@ int povray_render_animation(UnitCell *cell, RefList *list, unsigned int nproc,
 			signed int he, ke, le;
 			float x, y, z;
 
-			get_equiv(sp, j, h, k, l, &he, &ke, &le);
+			get_equiv(sym, m, j, h, k, l, &he, &ke, &le);
 
 			x = asx*he + bsx*ke + csx*le;
 			y = asy*he + bsy*ke + csy*le;
@@ -311,10 +309,10 @@ int povray_render_animation(UnitCell *cell, RefList *list, unsigned int nproc,
 
 		}
 
-		free_symoplist(sp);
-
 	}
 
+	free_symopmask(m);
+
 	fprintf(fh, "\n");
 	fclose(fh);
 
diff --git a/src/powder_plot.c b/src/powder_plot.c
index 6065da07..d073f2ef 100644
--- a/src/powder_plot.c
+++ b/src/powder_plot.c
@@ -318,6 +318,9 @@ static unsigned int process_hkl(struct image *image, const SymOpList *sym,
 	int h, k, l, redundancy;
 	double q, intensity;
 	unsigned int nref;
+	SymOpMask *m;
+
+	m = new_symopmask(sym);
 
 	nref = num_reflections(image->reflections);
 
@@ -330,9 +333,8 @@ static unsigned int process_hkl(struct image *image, const SymOpList *sym,
 		if ( use_redundancy ) {
 			redundancy = get_redundancy(refl);
 		} else {
-			SymOpList *sp = special_position(sym, h, k, l);
-			redundancy = num_equivs(sp);
-			free_symoplist(sp);
+			special_position(sym, m, h, k, l);
+			redundancy = num_equivs(sym, m);
 		}
 
 		/* Multiply by 2 to get 1/d (in m^-1) */
@@ -348,6 +350,8 @@ static unsigned int process_hkl(struct image *image, const SymOpList *sym,
 
 	}
 
+	free_symopmask(m);
+
 	return n_peaks;
 }
 
diff --git a/src/process_hkl.c b/src/process_hkl.c
index 06c5a54b..d34ea60a 100644
--- a/src/process_hkl.c
+++ b/src/process_hkl.c
@@ -521,6 +521,7 @@ int main(int argc, char *argv[])
 
 	if ( sym_str == NULL ) sym_str = strdup("1");
 	sym = get_pointgroup(sym_str);
+	STATUS("%s -> %p\n", sym_str, sym);
 	free(sym_str);
 
 	/* Open the data stream */
@@ -553,7 +554,7 @@ int main(int argc, char *argv[])
 			ERROR("Invalid indices for '--histogram'\n");
 			return 1;
 		}
-		space_for_hist = n_total_patterns * num_equivs(sym);
+		space_for_hist = n_total_patterns * num_equivs(sym, NULL);
 		hist_vals = malloc(space_for_hist * sizeof(double));
 		free(histo);
 		STATUS("Histogramming %i %i %i -> ", hist_h, hist_k, hist_l);
diff --git a/src/reflist-utils.c b/src/reflist-utils.c
index c0784dbb..5a9ec44f 100644
--- a/src/reflist-utils.c
+++ b/src/reflist-utils.c
@@ -109,8 +109,14 @@ int check_list_symmetry(RefList *list, const SymOpList *sym)
 	unsigned char *flags;
 	Reflection *refl;
 	RefListIterator *iter;
+	SymOpMask *mask;
 
 	flags = flags_from_list(list);
+	mask = new_symopmask(sym);
+	if ( mask == NULL ) {
+		ERROR("Couldn't create mask for list symmetry check.\n");
+		return 1;
+	}
 
 	for ( refl = first_refl(list, &iter);
 	      refl != NULL;
@@ -119,13 +125,19 @@ int check_list_symmetry(RefList *list, const SymOpList *sym)
 		int j;
 		int found = 0;
 		signed int h, k, l;
+		int n;
 
 		get_indices(refl, &h, &k, &l);
 
-		for ( j=0; j<num_equivs(sym); j++ ) {
+		special_position(sym, mask, h, k, l);
+		n = num_equivs(sym, mask);
+		STATUS("%i equivs: %3i %3i %3i\n", n, h, k, l);
+
+		for ( j=0; j<n; j++ ) {
 
 			signed int he, ke, le;
-			get_equiv(sym, j, h, k, l, &he, &ke, &le);
+			get_equiv(sym, mask, j, h, k, l, &he, &ke, &le);
+			STATUS("%3i: %3i %3i %3i\n", j, he, ke, le);
 
 			if ( abs(he) > INDMAX ) continue;
 			if ( abs(le) > INDMAX ) continue;
@@ -137,12 +149,15 @@ int check_list_symmetry(RefList *list, const SymOpList *sym)
 
 		if ( found > 1 ) {
 			free(flags);
+			free_symopmask(mask);
+			STATUS("Found %i %i %i twice\n", h, k, l);
 			return 1;  /* Symmetry is wrong! */
 		}
 
 	}
 
 	free(flags);
+	free_symopmask(mask);
 
 	return 0;
 }
@@ -155,11 +170,11 @@ int find_equiv_in_list(RefList *list, signed int h, signed int k,
 	int i;
 	int found = 0;
 
-	for ( i=0; i<num_equivs( sym); i++ ) {
+	for ( i=0; i<num_equivs(sym, NULL); i++ ) {
 
 		signed int he, ke, le;
 		Reflection *f;
-		get_equiv(sym, i, h, k, l, &he, &ke, &le);
+		get_equiv(sym, NULL, i, h, k, l, &he, &ke, &le);
 		f = find_refl(list, he, ke, le);
 
 		/* There must only be one equivalent.  If there are more, it
diff --git a/src/render_hkl.c b/src/render_hkl.c
index 823b6f1a..18d6bb81 100644
--- a/src/render_hkl.c
+++ b/src/render_hkl.c
@@ -96,12 +96,15 @@ static void draw_circles(signed int xh, signed int xk, signed int xl,
 {
 	Reflection *refl;
 	RefListIterator *iter;
+	SymOpMask *m;
 
 	if ( dctx == NULL ) {
 		*max_u = 0.0;  *max_v = 0.0;  *max_val = 0.0;
 		*max_res = 0.0;  *max_ux = 0;  *max_uy = 0;
 	}
 
+	m = new_symopmask(sym);
+
 	/* Iterate over all reflections */
 	for ( refl = first_refl(list, &iter);
 	      refl != NULL;
@@ -111,18 +114,17 @@ static void draw_circles(signed int xh, signed int xk, signed int xl,
 		signed int ha, ka, la;
 		int xi, yi;
 		int i, n;
-		SymOpList *sp;
 
 		get_indices(refl, &ha, &ka, &la);
 
-		sp = special_position(sym, ha, ka, la);
-		n = num_equivs(sp);
+		special_position(sym, m, ha, ka, la);
+		n = num_equivs(sym, m);
 
 		for ( i=0; i<n; i++ ) {
 
 			signed int h, k, l;
 
-			get_equiv(sp, i, ha, ka, la, &h, &k, &l);
+			get_equiv(sym, m, i, ha, ka, la, &h, &k, &l);
 
 			/* Is the reflection in the zone? */
 			if ( h*zh + k*zk + l*zl != 0 ) continue;
@@ -194,9 +196,9 @@ static void draw_circles(signed int xh, signed int xk, signed int xl,
 
 		}
 
-		free_symoplist(sp);
-
 	}
+
+	free_symopmask(m);
 }
 
 
diff --git a/src/symmetry.c b/src/symmetry.c
index 19daf0ec..3b995967 100644
--- a/src/symmetry.c
+++ b/src/symmetry.c
@@ -37,19 +37,6 @@
  */
 
 
-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;
@@ -84,6 +71,13 @@ struct _symoplist
 };
 
 
+struct _symopmask
+{
+	const SymOpList *list;
+	int *mask;
+};
+
+
 
 static void alloc_ops(SymOpList *ops)
 {
@@ -92,6 +86,35 @@ static void alloc_ops(SymOpList *ops)
 }
 
 
+/**
+ * new_symopmask:
+ *
+ * Returns: a new %SymOpMask, which you can use when filtering out special
+ * reflections.
+ **/
+SymOpMask *new_symopmask(const SymOpList *list)
+{
+	SymOpMask *m;
+	int i;
+
+	m = malloc(sizeof(struct _symopmask));
+	if ( m == NULL ) return NULL;
+
+	m->list = list;
+	m->mask = malloc(sizeof(int)*list->n_ops);
+	if ( m->mask == NULL ) {
+		free(m);
+		return NULL;
+	}
+
+	for ( i=0; i<list->n_ops; i++ ) {
+		m->mask[i] = 1;
+	}
+
+	return m;
+}
+
+
 /* Creates a new SymOpList */
 static SymOpList *new_symoplist()
 {
@@ -101,6 +124,7 @@ static SymOpList *new_symoplist()
 	new->max_ops = 16;
 	new->n_ops = 0;
 	new->ops = NULL;
+	new->divisors = NULL;
 	new->name = NULL;
 	new->num_equivs = 1;
 	alloc_ops(new);
@@ -128,88 +152,33 @@ void free_symoplist(SymOpList *ops)
 	free(ops);
 }
 
-
-static int is_identity(signed int *h, signed int *k, signed int *l)
+/**
+ * free_symopmask:
+ *
+ * Frees a %SymOpMask and all associated resources.
+ **/
+void free_symopmask(SymOpMask *m)
 {
-	if ( (h[0]!=1) || (h[1]!=0) || (h[2]!=0) ) return 0;
-	if ( (k[0]!=0) || (k[1]!=1) || (k[2]!=0) ) return 0;
-	if ( (l[0]!=0) || (l[1]!=0) || (l[2]!=1) ) return 0;
-	return 1;
+	if ( m == NULL ) return;
+	free(m->mask);
+	free(m);
 }
 
 
-/* Calculate the order of the operation "M", which is the lowest
- * integer n such that M^n = I. */
-static int order_of_op(signed int *hin, signed int *kin, signed int *lin)
-{
-	int n;
-	signed int h[3];
-	signed int k[3];
-	signed int l[3];
-
-	memcpy(h, hin, 3*sizeof(signed int));
-	memcpy(k, kin, 3*sizeof(signed int));
-	memcpy(l, lin, 3*sizeof(signed int));
-
-	for ( n=1; n<6; n++ ) {
-
-		signed int hnew[3];
-		signed int knew[3];
-		signed int lnew[3];
-
-		/* Yay matrices */
-		hnew[0] = h[0]*h[0] + h[1]*k[0] + h[2]*l[0];
-		hnew[1] = h[0]*h[1] + h[1]*k[1] + h[2]*l[1];
-		hnew[2] = h[0]*h[2] + h[1]*k[2] + h[2]*l[2];
-
-		knew[0] = k[0]*h[0] + k[1]*k[0] + k[2]*l[0];
-		knew[1] = k[0]*h[1] + k[1]*k[1] + k[2]*l[1];
-		knew[2] = k[0]*h[2] + k[1]*k[2] + k[2]*l[2];
-
-		lnew[0] = l[0]*h[0] + l[1]*k[0] + l[2]*l[0];
-		lnew[1] = l[0]*h[1] + l[1]*k[1] + l[2]*l[1];
-		lnew[2] = l[0]*h[2] + l[1]*k[2] + l[2]*l[2];
-
-		if ( is_identity(hnew, knew, lnew) ) break;
-
-		memcpy(h, hnew, 3*sizeof(signed int));
-		memcpy(k, knew, 3*sizeof(signed int));
-		memcpy(l, lnew, 3*sizeof(signed int));
-
-	}
-
-	return n+1;
-}
-
-
-/* This returns the number of operations in "ops".  To get the number of
- * symmetric equivalents this generates, use num_equivs() instead. */
+/* This returns the number of operations in "ops".  This might be different
+ * to num_equivs() if the point group is being constructed. */
 static int num_ops(const SymOpList *ops)
 {
 	return ops->n_ops;
 }
 
 
-static void update_num_equivs(SymOpList *ops)
-{
-	int i, n, tot;
-
-	n = num_ops(ops);
-	tot = 1;
-
-	for ( i=0; i<n; i++ ) {
-		tot *= ops->ops[i].order;
-	}
-
-	ops->num_equivs = tot;
-}
-
-
 /* Add a operation to a SymOpList */
 static void add_symop(SymOpList *ops,
-                      signed int *h, signed int *k, signed int *l)
+                      signed int *h, signed int *k, signed int *l,
+                      int order)
 {
-	int n, i;
+	int n;
 
 	if ( ops->n_ops == ops->max_ops ) {
 		/* Pretty sure this never happens, but still... */
@@ -221,27 +190,7 @@ static void add_symop(SymOpList *ops,
 	ops->ops[n].h = h;
 	ops->ops[n].k = k;
 	ops->ops[n].l = l;
-	ops->ops[n].order = order_of_op(h, k, l);
-	ops->n_ops++;
-
-	ops->divisors[0] = 1;
-	for ( i=1; i<ops->n_ops; i++ ) {
-		ops->divisors[i] = ops->divisors[i-1]*ops->ops[i].order;
-	}
-
-	update_num_equivs(ops);
-}
-
-
-static void add_copied_symop(SymOpList *ops, struct sym_op *copyme)
-{
-	if ( ops->n_ops == ops->max_ops ) {
-		/* Pretty sure this never happens, but still... */
-		ops->max_ops += 16;
-		alloc_ops(ops);
-	}
-
-	memcpy(&ops->ops[ops->n_ops], copyme, sizeof(*copyme));
+	ops->ops[n].order = order;
 	ops->n_ops++;
 }
 
@@ -252,9 +201,9 @@ static void add_copied_symop(SymOpList *ops, struct sym_op *copyme)
  * Returns: the number of equivalent reflections for a general reflection
  * in point group "ops".
  **/
-int num_equivs(const SymOpList *ops)
+int num_equivs(const SymOpList *ops, const SymOpMask *m)
 {
-	return ops->num_equivs;
+	return num_ops(ops);
 }
 
 
@@ -267,22 +216,167 @@ static signed int *v(signed int h, signed int k, signed int i, signed int l)
 }
 
 
+static void combine_ops(signed int *h1, signed int *k1, signed int *l1,
+                        signed int *h2, signed int *k2, signed int *l2,
+                        signed int *hnew, signed int *knew, signed int *lnew)
+{
+	/* Yay matrices */
+	hnew[0] = h1[0]*h2[0] + h1[1]*k2[0] + h1[2]*l2[0];
+	hnew[1] = h1[0]*h2[1] + h1[1]*k2[1] + h1[2]*l2[1];
+	hnew[2] = h1[0]*h2[2] + h1[1]*k2[2] + h1[2]*l2[2];
+
+	knew[0] = k1[0]*h2[0] + k1[1]*k2[0] + k1[2]*l2[0];
+	knew[1] = k1[0]*h2[1] + k1[1]*k2[1] + k1[2]*l2[1];
+	knew[2] = k1[0]*h2[2] + k1[1]*k2[2] + k1[2]*l2[2];
+
+	lnew[0] = l1[0]*h2[0] + l1[1]*k2[0] + l1[2]*l2[0];
+	lnew[1] = l1[0]*h2[1] + l1[1]*k2[1] + l1[2]*l2[1];
+	lnew[2] = l1[0]*h2[2] + l1[1]*k2[2] + l1[2]*l2[2];
+}
+
+
+static void expand_all_ops(signed int *hs, signed int *ks, signed int *ls,
+                           int skip, SymOpList *s, SymOpList *expanded)
+{
+	int i, n;
+
+	n = num_ops(s);
+
+	for ( i=0; i<n; i++ ) {
+
+		int oi;
+		struct sym_op *opi = &s->ops[i];
+
+		for ( oi=0; oi<opi->order; oi++ ) {
+
+			int oi_it;
+			signed int *h_ordered, *k_ordered, *l_ordered;
+
+			h_ordered = malloc(3*sizeof(signed int));
+			k_ordered = malloc(3*sizeof(signed int));
+			l_ordered = malloc(3*sizeof(signed int));
+			assert(h_ordered != NULL);
+			assert(k_ordered != NULL);
+			assert(l_ordered != NULL);
+
+			memcpy(h_ordered, hs, 3*sizeof(signed int));
+			memcpy(k_ordered, ks, 3*sizeof(signed int));
+			memcpy(l_ordered, ls, 3*sizeof(signed int));
+
+			/* Apply "opi" this number of times */
+			for ( oi_it=0; oi_it<oi; oi_it++ ) {
+
+				signed int hfs[3], kfs[3], lfs[3];
+
+				combine_ops(h_ordered, k_ordered, l_ordered,
+				            opi->h, opi->k, opi->l,
+				            hfs, kfs, lfs);
+
+				if ( i != skip ) {
+
+					memcpy(h_ordered, hfs,
+					       3*sizeof(signed int));
+					memcpy(k_ordered, kfs,
+					       3*sizeof(signed int));
+					memcpy(l_ordered, lfs,
+					       3*sizeof(signed int));
+
+				}
+			}
+
+			STATUS("i=%i, oi=%i\n", i, oi);
+			STATUS("Creating %3i %3i %3i\n", h_ordered[0],
+			                                 h_ordered[1],
+			                                 h_ordered[2]);
+			STATUS("         %3i %3i %3i\n", k_ordered[0],
+			                                 k_ordered[1],
+			                                 k_ordered[2]);
+			STATUS("         %3i %3i %3i\n", l_ordered[0],
+			                                 l_ordered[1],
+			                                 l_ordered[2]);
+			STATUS("\n");
+			add_symop(expanded, h_ordered, k_ordered, l_ordered, 1);
+
+		}
+
+	}
+}
+
+
+/* Fill in the other operations for a point group starting from its
+ * generators */
+static SymOpList *expand_ops(SymOpList *s)
+{
+	int n, i;
+	SymOpList *expanded;
+
+	n = num_ops(s);
+	expanded = new_symoplist();
+	if ( expanded == NULL ) return NULL;
+	expanded->name = strdup(symmetry_name(s));
+	STATUS("%i ops to expand.\n", n);
+
+	for ( i=0; i<n; i++ ) {
+
+		int oi;
+		struct sym_op *opi = &s->ops[i];
+		STATUS("Op %i, order=%i\n", i, opi->order);
+
+		for ( oi=0; oi<opi->order; oi++ ) {
+
+			int oi_it;
+			signed int h_ordered[3], k_ordered[3], l_ordered[3];
+
+			h_ordered[0] = 1;  h_ordered[1] = 0;  h_ordered[2] = 0;
+			k_ordered[0] = 0;  k_ordered[1] = 1;  k_ordered[2] = 0;
+			l_ordered[0] = 0;  l_ordered[1] = 0;  l_ordered[2] = 1;
+
+			/* Apply "opi" this number of times */
+			for ( oi_it=0; oi_it<oi; oi_it++ ) {
+
+				signed int hfs[3], kfs[3], lfs[3];
+
+				combine_ops(h_ordered, k_ordered, l_ordered,
+				            opi->h, opi->k, opi->l,
+				            hfs, kfs, lfs);
+
+				memcpy(h_ordered, hfs, 3*sizeof(signed int));
+				memcpy(k_ordered, kfs, 3*sizeof(signed int));
+				memcpy(l_ordered, lfs, 3*sizeof(signed int));
+
+			}
+
+			STATUS("Upper: i=%i, oi=%i\n", i, oi);
+			expand_all_ops(h_ordered, k_ordered, l_ordered, i,
+			               s, expanded);
+
+		}
+
+	}
+
+	free_symoplist(s);
+
+	return expanded;
+}
+
+
 /********************************* 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));  /* -I */
+	add_symop(new, v(1,0,0,0), v(0,1,0,0), v(0,0,0,1), 2);  /* -I */
 	new->name = strdup("-1");
-	return new;
+	return expand_ops(new);
 }
 
 
 static SymOpList *make_1()
 {
 	SymOpList *new = new_symoplist();
+	add_symop(new, v(1,0,0,0), v(0,1,0,0), v(0,0,0,1), 1);  /* I */
 	new->name = strdup("1");
-	return new;
+	return expand_ops(new);
 }
 
 
@@ -291,28 +385,28 @@ static SymOpList *make_1()
 static SymOpList *make_2m()
 {
 	SymOpList *new = new_symoplist();
-	add_symop(new, v(-1,0,0,0), v(0,1,0,0), v(0,0,0,-1)); /* 2 */
-	add_symop(new, v(1,0,0,0), v(0,-1,0,0), v(0,0,0,1));  /* m */
+	add_symop(new, v(-1,0,0,0), v(0,1,0,0), v(0,0,0,-1), 2); /* 2 */
+	add_symop(new, v(1,0,0,0), v(0,-1,0,0), v(0,0,0,1), 2);  /* m */
 	new->name = strdup("2/m");
-	return new;
+	return expand_ops(new);
 }
 
 
 static SymOpList *make_2()
 {
 	SymOpList *new = new_symoplist();
-	add_symop(new, v(-1,0,0,0), v(0,1,0,0), v(0,0,0,-1)); /* 2 */
+	add_symop(new, v(-1,0,0,0), v(0,1,0,0), v(0,0,0,-1), 2); /* 2 */
 	new->name = strdup("2");
-	return new;
+	return expand_ops(new);
 }
 
 
 static SymOpList *make_m()
 {
 	SymOpList *new = new_symoplist();
-	add_symop(new, v(1,0,0,0), v(0,-1,0,0), v(0,0,0,1));  /* m */
+	add_symop(new, v(1,0,0,0), v(0,-1,0,0), v(0,0,0,1), 2);  /* m */
 	new->name = strdup("m");
-	return new;
+	return expand_ops(new);
 }
 
 
@@ -321,31 +415,31 @@ static SymOpList *make_m()
 static SymOpList *make_mmm()
 {
 	SymOpList *new = new_symoplist();
-	add_symop(new, v(-1,0,0,0), v(0,-1,0,0), v(0,0,0,1)); /* 2 */
-	add_symop(new, v(-1,0,0,0), v(0,1,0,0), v(0,0,0,-1)); /* 2 */
-	add_symop(new, v(1,0,0,0), v(0,-1,0,0), v(0,0,0,1));  /* m */
+	add_symop(new, v(-1,0,0,0), v(0,-1,0,0), v(0,0,0,1), 2); /* 2 */
+	add_symop(new, v(-1,0,0,0), v(0,1,0,0), v(0,0,0,-1), 2); /* 2 */
+	add_symop(new, v(1,0,0,0), v(0,-1,0,0), v(0,0,0,1), 2);  /* m */
 	new->name = strdup("mmm");
-	return new;
+	return expand_ops(new);
 }
 
 
 static SymOpList *make_222()
 {
 	SymOpList *new = new_symoplist();
-	add_symop(new, v(-1,0,0,0), v(0,-1,0,0), v(0,0,0,1)); /* 2 */
-	add_symop(new, v(-1,0,0,0), v(0,1,0,0), v(0,0,0,-1)); /* 2 */
+	add_symop(new, v(-1,0,0,0), v(0,-1,0,0), v(0,0,0,1), 2); /* 2 */
+	add_symop(new, v(-1,0,0,0), v(0,1,0,0), v(0,0,0,-1), 2); /* 2 */
 	new->name = strdup("222");
-	return new;
+	return expand_ops(new);
 }
 
 
 static SymOpList *make_mm2()
 {
 	SymOpList *new = new_symoplist();
-	add_symop(new, v(-1,0,0,0), v(0,-1,0,0), v(0,0,0,1)); /* 2 */
-	add_symop(new, v(1,0,0,0), v(0,-1,0,0), v(0,0,0,1));  /* m */
+	add_symop(new, v(-1,0,0,0), v(0,-1,0,0), v(0,0,0,1), 2); /* 2 */
+	add_symop(new, v(1,0,0,0), v(0,-1,0,0), v(0,0,0,1), 2);  /* m */
 	new->name = strdup("mm2");
-	return new;
+	return expand_ops(new);
 }
 
 
@@ -360,7 +454,7 @@ static SymOpList *make_4m()
 static SymOpList *make_4()
 {
 	SymOpList *new = new_symoplist();
-	add_symop(new, v(0,1,0,0), v(1,0,0,0), v(0,0,0,1)); /* 4 */
+	add_symop(new, v(0,1,0,0), v(1,0,0,0), v(0,0,0,1), 4); /* 4 */
 	return NULL;
 }
 
@@ -488,8 +582,8 @@ static void do_op(const struct sym_op *op,
                   signed int *he, signed int *ke, signed int *le)
 {
 	*he = h*op->h[0] + k*op->h[1] + l*op->h[2];
-	*ke = h*op->k[0] + k*op->h[1] + l*op->k[2];
-	*le = h*op->l[0] + k*op->h[1] + l*op->l[2];
+	*ke = h*op->k[0] + k*op->k[1] + l*op->k[2];
+	*le = h*op->l[0] + k*op->l[1] + l*op->l[2];
 }
 
 
@@ -513,37 +607,31 @@ static void do_op(const struct sym_op *op,
  * given reflection is a special high-symmetry one), call special_position()
  * first to get a "specialised" SymOpList and use that instead.
  **/
-void get_equiv(const SymOpList *ops, int idx,
+void get_equiv(const SymOpList *ops, const SymOpMask *m, int idx,
                signed int h, signed int k, signed int l,
                signed int *he, signed int *ke, signed int *le)
 {
-	int sig[32];
-	int i, n, r;
+	int i, n;
 
 	n = num_ops(ops);
-	r = idx;
-	for ( i=n-1; i>=0; i-- ) {
-		sig[i] = r / ops->divisors[i];
-		r = r % ops->divisors[i];
-		assert(sig[i] < ops->ops[i].order);
+	for ( i=idx; i<n; i++ ) {
+		if ( (m == NULL) || m->mask[i] ) {
+			do_op(&ops->ops[i], h, k, l, he, ke, le);
+			return;
+		}
 	}
 
-	for ( i=0; i<n; i++ ) {
-
-		int s;
+	ERROR("Index %i out of range for point group '%s'\n", idx,
+	      symmetry_name(ops));
 
-		/* Do this operation "sig[i]" times */
-		for ( s=0; s<sig[i]; s++ ) {
-			do_op(&ops->ops[i], h, k, l, &h, &k, &l);
-		}
-
-	}
+	*he = 0;  *ke = 0;  *le = 0;
 }
 
 
 /**
  * special_position:
  * @ops: A %SymOpList, usually corresponding to a point group
+ * @m: A %SymOpMask created with new_symopmask()
  * @h: index of a reflection
  * @k: index of a reflection
  * @l: index of a reflection
@@ -554,26 +642,47 @@ void get_equiv(const SymOpList *ops, int idx,
  *
  * Returns: the "specialised" %SymOpList.
  **/
-SymOpList *special_position(const SymOpList *ops,
-                            signed int h, signed int k, signed int l)
+void special_position(const SymOpList *ops, SymOpMask *m,
+                      signed int h, signed int k, signed int l)
 {
 	int i, n;
-	SymOpList *specialised;
+	signed int *htest;
+	signed int *ktest;
+	signed int *ltest;
 
-	n = num_ops(ops);
-	specialised = new_symoplist();
+	assert(m->list = ops);
+
+	n = num_equivs(ops, NULL);
+	htest = malloc(n*sizeof(signed int));
+	ktest = malloc(n*sizeof(signed int));
+	ltest = malloc(n*sizeof(signed int));
 
 	for ( i=0; i<n; i++ ) {
 
-		signed int ht, kt, lt;
+		signed int he, ke, le;
+		int j;
+
+		get_equiv(ops, NULL, i, h, k, l, &he, &ke, &le);
+
+		m->mask[i] = 1;
+		for ( j=0; j<i; j++ ) {
+			if ( (he==htest[j]) && (ke==ktest[j])
+			  && (le==ltest[j]) )
+			{
+				m->mask[i] = 0;
+				break;  /* Only need to find one */
+			}
+		}
 
-		do_op(&ops->ops[i], h, k, l, &ht, &kt, &lt);
-		if ( (h==ht) || (k==kt) || (l==lt) ) continue;
-		add_copied_symop(specialised, &ops->ops[i]);
+		htest[i] = he;
+		ktest[i] = ke;
+		ltest[i] = le;
 
 	}
 
-	return specialised;
+	free(htest);
+	free(ktest);
+	free(ltest);
 }
 
 
@@ -585,12 +694,12 @@ void get_asymm(const SymOpList *ops,
 	int p;
 	signed int best_h, best_k, best_l;
 
-	nequiv = num_equivs(ops);
+	nequiv = num_equivs(ops, NULL);
 
 	best_h = h;  best_k = k;  best_l = l;
 	for ( p=0; p<nequiv; p++ ) {
 
-		get_equiv(ops, p, h, k, l, hp, kp, lp);
+		get_equiv(ops, NULL, p, h, k, l, hp, kp, lp);
 
 		if ( h > best_h ) {
 			best_h = h;  best_k = k;  best_l = l;
diff --git a/src/symmetry.h b/src/symmetry.h
index 4dcaae06..45ffe318 100644
--- a/src/symmetry.h
+++ b/src/symmetry.h
@@ -24,22 +24,27 @@
  **/
 typedef struct _symoplist SymOpList;
 
-extern void free_symoplist(SymOpList *ops);
+/**
+ * SymOpMask
+ *
+ * Opaque type.
+ **/
+typedef struct _symopmask SymOpMask;
 
+extern void free_symoplist(SymOpList *ops);
 extern SymOpList *get_pointgroup(const char *sym);
-
 extern const char *symmetry_name(const SymOpList *ops);
 
-extern SymOpList *special_position(const SymOpList *ops,
-                                   signed int h, signed int k, signed int l);
+extern SymOpMask *new_symopmask(const SymOpList *ops);
+extern void free_symopmask(SymOpMask *m);
 
+extern void special_position(const SymOpList *ops, SymOpMask *m,
+                             signed int h, signed int k, signed int l);
 extern void get_asymm(const SymOpList *ops,
                       signed int h, signed int k, signed int l,
                       signed int *hp, signed int *kp, signed int *lp);
-
-extern int num_equivs(const SymOpList *ops);
-
-extern void get_equiv(const SymOpList *ops, int idx,
+extern int num_equivs(const SymOpList *ops, const SymOpMask *m);
+extern void get_equiv(const SymOpList *ops, const SymOpMask *m, int idx,
                       signed int h, signed int k, signed int l,
                       signed int *he, signed int *ke, signed int *le);