Fix symmetry so that it works, and use it in render_hkl

5 files changed, 148 insertions, 78 deletions

diff --git a/src/Makefile.am b/src/Makefile.am
index 259542fb..8d4cdec2 100644
--- a/src/Makefile.am
+++ b/src/Makefile.am
@@ -46,7 +46,8 @@ powder_plot_SOURCES = powder_plot.c cell.c utils.c image.c hdf5-file.c \
                       detector.c
 powder_plot_LDADD = @LIBS@
 
-render_hkl_SOURCES = render_hkl.c cell.c reflections.c utils.c povray.c
+render_hkl_SOURCES = render_hkl.c cell.c reflections.c utils.c povray.c \
+                     symmetry.c
 render_hkl_LDADD = @LIBS@
 
 calibrate_detector_SOURCES = calibrate_detector.c utils.c hdf5-file.c image.c \
diff --git a/src/Makefile.in b/src/Makefile.in
index 950b5789..5de55c96 100644
--- a/src/Makefile.in
+++ b/src/Makefile.in
@@ -108,7 +108,8 @@ am_process_hkl_OBJECTS = process_hkl.$(OBJEXT) sfac.$(OBJEXT) \
 process_hkl_OBJECTS = $(am_process_hkl_OBJECTS)
 process_hkl_DEPENDENCIES =
 am_render_hkl_OBJECTS = render_hkl.$(OBJEXT) cell.$(OBJEXT) \
-	reflections.$(OBJEXT) utils.$(OBJEXT) povray.$(OBJEXT)
+	reflections.$(OBJEXT) utils.$(OBJEXT) povray.$(OBJEXT) \
+	symmetry.$(OBJEXT)
 render_hkl_OBJECTS = $(am_render_hkl_OBJECTS)
 render_hkl_DEPENDENCIES =
 DEFAULT_INCLUDES = -I.@am__isrc@ -I$(top_builddir)
@@ -261,7 +262,9 @@ powder_plot_SOURCES = powder_plot.c cell.c utils.c image.c hdf5-file.c \
                       detector.c
 
 powder_plot_LDADD = @LIBS@
-render_hkl_SOURCES = render_hkl.c cell.c reflections.c utils.c povray.c
+render_hkl_SOURCES = render_hkl.c cell.c reflections.c utils.c povray.c \
+                     symmetry.c
+
 render_hkl_LDADD = @LIBS@
 calibrate_detector_SOURCES = calibrate_detector.c utils.c hdf5-file.c image.c \
                              filters.c peaks.c detector.c cell.c diffraction.c \
diff --git a/src/render_hkl.c b/src/render_hkl.c
index 06ac745d..106d1b58 100644
--- a/src/render_hkl.c
+++ b/src/render_hkl.c
@@ -26,6 +26,7 @@
 #include "utils.h"
 #include "reflections.h"
 #include "povray.h"
+#include "symmetry.h"
 
 
 static void show_help(const char *s)
@@ -55,6 +56,7 @@ static void render_za(UnitCell *cell, double *ref, unsigned int *c)
 	double csx, csy, csz;
 	signed int h, k;
 	float wh, ht;
+	const char *sym = "6/mmm";
 
 	wh = 1024;
 	ht = 1024;
@@ -94,22 +96,34 @@ static void render_za(UnitCell *cell, double *ref, unsigned int *c)
 
 		double u, v, intensity, res;
 		int ct;
+		int nequiv, p;
 
 		ct = lookup_count(c, h, k, 0);
 		if ( ct < 1 ) continue;
 
 		intensity = lookup_intensity(ref, h, k, 0) / (float)ct;
 
-		res = resolution(cell, h, k, 0);
-		if ( res > max_res ) max_res = res;
-
 		if ( intensity != 0 ) {
-			u =  (double)h*as*sin(theta);
-			v =  (double)h*as*cos(theta) + k*bs;
-			if ( fabs(u) > fabs(max_u) ) max_u = fabs(u);
-			if ( fabs(v) > fabs(max_v) ) max_v = fabs(v);
-			if ( fabs(intensity) > fabs(max_intensity) )
-						max_intensity = fabs(intensity);
+
+			nequiv = num_equivs(h, k, 0, sym);
+			for ( p=0; p<nequiv; p++ ) {
+
+				signed int he, ke, le;
+				get_equiv(h, k, 0, &he, &ke, &le, sym, p);
+
+				u =  (double)he*as*sin(theta);
+				v =  (double)he*as*cos(theta) + ke*bs;
+				if ( fabs(u) > fabs(max_u) ) max_u = fabs(u);
+				if ( fabs(v) > fabs(max_v) ) max_v = fabs(v);
+				if ( fabs(intensity) > fabs(max_intensity) ) {
+					max_intensity = fabs(intensity);
+				}
+				res = resolution(cell, he, ke, 0);
+				if ( res > max_res ) max_res = res;
+
+
+
+			}
 		}
 
 	}
@@ -140,6 +154,7 @@ static void render_za(UnitCell *cell, double *ref, unsigned int *c)
 
 		double u, v, intensity, val;
 		int ct;
+		int nequiv, p;
 
 		ct = lookup_count(c, h, k, 0);
 		if ( ct < 1 ) continue;
@@ -147,14 +162,23 @@ static void render_za(UnitCell *cell, double *ref, unsigned int *c)
 		intensity = lookup_intensity(ref, h, k, 0) / (float)ct;
 		val = 3.0*intensity/max_intensity;
 
-		u = (double)h*as*sin(theta);
-		v = (double)h*as*cos(theta) + k*bs;
+		nequiv = num_equivs(h, k, 0, sym);
+		for ( p=0; p<nequiv; p++ ) {
+
+			signed int he, ke, le;
+			get_equiv(h, k, 0, &he, &ke, &le, sym, p);
+
+			u = (double)he*as*sin(theta);
+			v = (double)he*as*cos(theta) + ke*bs;
 
-		cairo_arc(dctx, ((double)wh/2)+u*scale*2,
-				((double)ht/2)+v*scale*2, max_r, 0, 2*M_PI);
+			cairo_arc(dctx, ((double)wh/2)+u*scale*2,
+					((double)ht/2)+v*scale*2, max_r,
+					0, 2*M_PI);
 
-		cairo_set_source_rgb(dctx, val, val, val);
-		cairo_fill(dctx);
+			cairo_set_source_rgb(dctx, val, val, val);
+			cairo_fill(dctx);
+
+		}
 
 	}
 	}
diff --git a/src/symmetry.c b/src/symmetry.c
index bedc146c..06148de7 100644
--- a/src/symmetry.c
+++ b/src/symmetry.c
@@ -21,25 +21,29 @@
 #include "utils.h"
 
 
+#ifdef DEBUG
+#define SYM_DEBUG STATUS
+#else /* DEBUG */
+#define SYM_DEBUG(...)
+#endif /* DEBUG */
+
+
 /* Conditions for a reflection to be in the asymmetric unit cell */
 #define COND_1(h, k, l) (1)
-#define COND_6MMM(h, k, i, l) ( (h>=0) && (k>=0) && (l>=0) \
-                             && ((h>k)||((h==0)&&(k==0))) )
+#define COND_6(h, k, l) ( (h>=0) && (k>=0) )
+#define COND_6M(h, k, l) ( (h>=0) && (k>=0) && (l>=0) )
+#define COND_6MMM(h, k, l) ( (h>=0) && (k>=0) && (l>=0) \
+                          && (h>=k) )
 /* TODO: Add more groups here */
 
+
 /* Macros for checking the above conditions and returning if satisfied */
-#define CHECK_COND_FOURIDX(h, k, i, l, cond)       \
-	if ( COND_##cond((h), (k), (i), (l)) ) {   \
+#define CHECK_COND(h, k, l, cond)                  \
+	if ( COND_##cond((h), (k), (l)) ) {        \
 		*hp = (h);  *kp = (k);  *lp = (l); \
 		return;                            \
 	}
 
-#define CHECK_COND_THREEIDX(h, k, l, cond)          \
-	if ( COND_##cond((h), (k), (l)) ) {         \
-		*hp = (h);  *kp = (k);  *lp = (l);  \
-		return;                             \
-	}
-
 
 /* Abort macro if no match found */
 #define SYM_ABORT                                             \
@@ -49,50 +53,87 @@
 
 /* FIXME: Should take into account special indices
  * e.g. l==0 has fewer equivalent reflections */
-static int num_equivs4(signed int h, signed int k, signed int i, signed int l,
-                       const char *sym)
+int num_equivs(signed int h, signed int k, signed int l, const char *sym)
 {
 	if ( strcmp(sym, "6/mmm") == 0 ) return 24;
+	if ( strcmp(sym, "6") == 0 ) return 6;
+	if ( strcmp(sym, "6/m") == 0 ) return 12;
 	/* TODO: Add more groups here */
 
 	return 1;
 }
 
 
-static void get_equiv4(signed int h, signed int k, signed int i, signed int l,
-                       signed int *he, signed int *ke, signed int *ie,
-                       signed int *le, const char *sym, int idx)
+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)
 {
+	if ( strcmp(sym, "1") == 0 ) {
+		*he = h;   *ke = k;   *le = l;  return;
+	}
+
+	if ( strcmp(sym, "6") == 0 ) {
+		signed int i = -h-k;
+		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 ) {
+		signed int i = -h-k;
+		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;
+		}
+	}
+
 	if ( strcmp(sym, "6/mmm") == 0 ) {
+		signed int i = -h-k;
 		switch ( idx ) {
-		case  0 : *he = h;  *ke = k;  *ie = i;  *le = l;  return;
-		case  1 : *he = h;  *ke = i;  *ie = k;  *le = l;  return;
-		case  2 : *he = k;  *ke = h;  *ie = i;  *le = l;  return;
-		case  3 : *he = k;  *ke = i;  *ie = h;  *le = l;  return;
-		case  4 : *he = i;  *ke = h;  *ie = k;  *le = l;  return;
-		case  5 : *he = i;  *ke = k;  *ie = h;  *le = l;  return;
-		case  6 : *he = h;  *ke = k;  *ie = i;  *le = -l;  return;
-		case  7 : *he = h;  *ke = i;  *ie = k;  *le = -l;  return;
-		case  8 : *he = k;  *ke = h;  *ie = i;  *le = -l;  return;
-		case  9 : *he = k;  *ke = i;  *ie = h;  *le = -l;  return;
-		case 10 : *he = i;  *ke = h;  *ie = k;  *le = -l;  return;
-		case 11 : *he = i;  *ke = k;  *ie = h;  *le = -l;  return;
-		case 12 : *he = -h;  *ke = -k;  *ie = -i;  *le = l;  return;
-		case 13 : *he = -h;  *ke = -i;  *ie = -k;  *le = l;  return;
-		case 14 : *he = -k;  *ke = -h;  *ie = -i;  *le = l;  return;
-		case 15 : *he = -k;  *ke = -i;  *ie = -h;  *le = l;  return;
-		case 16 : *he = -i;  *ke = -h;  *ie = -k;  *le = l;  return;
-		case 17 : *he = -i;  *ke = -k;  *ie = -h;  *le = l;  return;
-		case 18 : *he = -h;  *ke = -k;  *ie = -i;  *le = -l;  return;
-		case 19 : *he = -h;  *ke = -i;  *ie = -k;  *le = -l;  return;
-		case 20 : *he = -k;  *ke = -h;  *ie = -i;  *le = -l;  return;
-		case 21 : *he = -k;  *ke = -i;  *ie = -h;  *le = -l;  return;
-		case 22 : *he = -i;  *ke = -h;  *ie = -k;  *le = -l;  return;
-		case 23 : *he = -i;  *ke = -k;  *ie = -h;  *le = -l;  return;
+		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;
 		}
 	}
 
-	*he = h;  *ke = k;  *ie = i;  *le = l;
+	/* Fallback for unrecognised groups */
+	*he = h;  *ke = k;  *le = l;
 }
 
 
@@ -100,28 +141,21 @@ void get_asymm(signed int h, signed int k, signed int l,
                signed int *hp, signed int *kp, signed int *lp,
                const char *sym)
 {
-	if ( strcmp(sym, "1") == 0 ) {
-		CHECK_COND_THREEIDX(h, k, l, 1);
-		SYM_ABORT;
+	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);
+		if ( strcmp(sym, "1") == 0 ) CHECK_COND(he, ke, le, 1);
+		if ( strcmp(sym, "6") == 0 ) CHECK_COND(he, ke, le, 6);
+		if ( strcmp(sym, "6/m") == 0 ) CHECK_COND(he, ke, le, 6M);
+		if ( strcmp(sym, "6/mmm") == 0 ) CHECK_COND(he, ke, le, 6MMM);
 	}
 
-	if ( strcmp(sym, "6/mmm") == 0 ) {
-
-		const signed int i = h+k;
-		int nequiv = num_equivs4(h, k, i, l, sym);
-		int p;
-
-		for ( p=0; p<nequiv; p++ ) {
-			signed int he, ke, ie, le;
-			get_equiv4(h, k, i, l, &he, &ke, &ie, &le, sym, p);
-			CHECK_COND_FOURIDX(he, ke, ie, le, 6MMM);
-		}
-
-		SYM_ABORT;  /* Should never reach here */
-
-	}
-
-	/* TODO: Add more groups here */
+	SYM_ABORT;  /* Should never reach here */
 
 	ERROR("Unknown point group '%s'\n", sym);
 	abort();
diff --git a/src/symmetry.h b/src/symmetry.h
index f1b6b127..ebe123c9 100644
--- a/src/symmetry.h
+++ b/src/symmetry.h
@@ -22,4 +22,12 @@ extern void get_asymm(signed int h, signed int k, signed int l,
                       signed int *hp, signed int *kp, signed int *lp,
                       const char *sym);
 
+extern int num_equivs(signed int h, signed int k, signed int l,
+                      const char *sym);
+
+extern 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);
+
+
 #endif	/* SYMMETRY_H */