Commit 7b4161e4 authored by Thomas White's avatar Thomas White
Browse files

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

parent da7591b2
......@@ -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 \
......
......@@ -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 \
......
......@@ -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);
}
}
}
......
......@@ -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();
......
......@@ -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 */
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment