aboutsummaryrefslogtreecommitdiff
diff options
context:
space:
mode:
authorThomas White <taw@bitwiz.org.uk>2010-07-25 22:36:15 +0200
committerThomas White <taw@physics.org>2012-02-22 15:26:54 +0100
commit5d22ee1971a09a12910e2d03ac7dabfc159c1924 (patch)
tree7f45ab29b7367986a7ce18737c1f656e4ad8fb5a
parent034b20d9356d108bd8d7e8a9f874120b9600cad2 (diff)
render_hkl: Allow other zone axes
-rw-r--r--src/render_hkl.c187
1 files changed, 126 insertions, 61 deletions
diff --git a/src/render_hkl.c b/src/render_hkl.c
index 5a3014ec..06075131 100644
--- a/src/render_hkl.c
+++ b/src/render_hkl.c
@@ -76,7 +76,8 @@ static void show_help(const char *s)
#ifdef HAVE_CAIRO
-static void render_za(UnitCell *cell, double *ref, unsigned int *c,
+static void render_za(UnitCell *cell, ReflItemList *items,
+ double *ref, unsigned int *c,
double boost, const char *sym, int wght, int colscale)
{
cairo_surface_t *surface;
@@ -85,16 +86,36 @@ static void render_za(UnitCell *cell, double *ref, unsigned int *c,
double scale_u, scale_v, scale;
double sep_u, sep_v, max_r;
double u, v;
- signed int max_h, max_k;
+ signed int max_ux, max_uy;
double as, bs, theta;
double asx, asy, asz;
double bsx, bsy, bsz;
double csx, csy, csz;
- signed int h, k;
float wh, ht;
+ signed int xh, xk, xl;
+ signed int yh, yk, yl;
+ signed int zh, zk, zl;
+ double xx, xy, xz;
+ double yx, yy, yz;
+ signed int xi, yi;
+ char tmp[256];
+ cairo_text_extents_t size;
+ double cx, cy;
+ const double border = 200.0;
+
+ xh = 1; xk = -1; xl = 0;
+ yh = 0; yk = 0; yl = 1;
+
+ zh = xk*yl - xl*yk;
+ zk = - xh*yl + xl*yh;
+ zl = xh*yk - xk*yh;
+
+ STATUS("Zone axis is %i %i %i\n", zh, zk, zl);
wh = 1024;
ht = 1024;
+ cx = 512.0;
+ cy = 500.0;
surface = cairo_pdf_surface_create("za.pdf", wh, ht);
@@ -112,7 +133,7 @@ static void render_za(UnitCell *cell, double *ref, unsigned int *c,
cairo_fill(dctx);
max_u = 0.0; max_v = 0.0; max_val = 0.0;
- max_res = 0.0; max_h = 0; max_k = 0;
+ max_res = 0.0; max_ux = 0; max_uy = 0;
/* Work out reciprocal lattice spacings and angles for this cut */
if ( cell_get_reciprocal(cell, &asx, &asy, &asz,
@@ -121,31 +142,47 @@ static void render_za(UnitCell *cell, double *ref, unsigned int *c,
ERROR("Couldn't get reciprocal parameters\n");
return;
}
- theta = angle_between(asx, asy, asz, bsx, bsy, bsz);
- as = modulus(asx, asy, asz) / 1e9;
- bs = modulus(bsx, bsy, bsz) / 1e9;
-
- for ( h=-INDMAX; h<INDMAX; h++ ) {
- for ( k=-INDMAX; k<INDMAX; k++ ) {
+ xx = xh*asx + xk*bsx + xl*csx;
+ xy = xh*asy + xk*bsy + xl*csy;
+ xz = xh*asz + xk*bsz + xl*csz;
+ yx = yh*asx + yk*bsx + yl*csx;
+ yy = yh*asy + yk*bsy + yl*csy;
+ yz = yh*asz + yk*bsz + yl*csz;
+ theta = angle_between(xx, xy, xz, yx, yy, yz);
+ as = modulus(xx, xy, xz) / 1e9;
+ bs = modulus(yx, yy, yz) / 1e9;
+
+ for ( xi=-INDMAX; xi<INDMAX; xi++ ) {
+ for ( yi=-INDMAX; yi<INDMAX; yi++ ) {
double u, v, val, res;
int ct;
int nequiv, p;
+ signed int h, k, l;
+
+ h = xi*xh + yi*yh;
+ k = xi*xk + yi*yk;
+ l = xi*xl + yi*yl;
- ct = lookup_count(c, h, k, 0);
- if ( ct < 1 ) continue;
+ /* Do we have this reflection? */
+ if ( find_unique_equiv(items, h, k, l, sym, &h, &k, &l) == 0 ) {
+ continue;
+ }
+
+ /* Reflection in the zone? */
+ if ( zh*h + zk*k + zl*l != 0 ) continue;
switch ( wght ) {
case WGHT_I :
- val = lookup_intensity(ref, h, k, 0);
+ val = lookup_intensity(ref, h, k, l);
break;
case WGHT_SQRTI :
- val = lookup_intensity(ref, h, k, 0);
+ val = lookup_intensity(ref, h, k, l);
val = (val>0.0) ? sqrt(val) : 0.0;
break;
case WGHT_COUNTS :
val = (float)ct;
- val /= (float)num_equivs(h, k, 0, sym);
+ val /= (float)num_equivs(h, k, l, sym);
break;
case WGHT_RAWCOUNTS :
val = (float)ct;
@@ -155,22 +192,28 @@ static void render_za(UnitCell *cell, double *ref, unsigned int *c,
abort();
}
- nequiv = num_equivs(h, k, 0, sym);
+ nequiv = num_equivs(h, k, l, sym);
for ( p=0; p<nequiv; p++ ) {
signed int he, ke, le;
- get_equiv(h, k, 0, &he, &ke, &le, sym, p);
+ signed int ux, uy;
+
+ get_equiv(h, k, l, &he, &ke, &le, sym, p);
+
+ ux = he*xh + ke*xk + le*xl;
+ uy = he*yh + ke*yk + le*yl;
+
+ u = (double)ux*as*sin(theta);
+ v = (double)ux*as*cos(theta) + uy*bs;
- 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(val) > fabs(max_val) ) {
max_val = fabs(val);
}
- if ( fabs(he) > max_h ) max_h = fabs(he);
- if ( fabs(ke) > max_k ) max_k = fabs(ke);
- res = resolution(cell, he, ke, 0);
+ if ( fabs(ux) > max_ux ) max_ux = fabs(ux);
+ if ( fabs(uy) > max_uy ) max_uy = fabs(uy);
+ res = resolution(cell, he, ke, le);
if ( res > max_res ) max_res = res;
}
@@ -189,36 +232,52 @@ static void render_za(UnitCell *cell, double *ref, unsigned int *c,
}
/* Choose whichever scaling factor gives the smallest value */
- scale_u = ((double)(wh/2.0)-50.0) / max_u;
- scale_v = ((double)(ht/2.0)-50.0) / max_v;
+ scale_u = ((double)wh-border) / (2.0*max_u);
+ scale_v = ((double)ht-border) / (2.0*max_v);
scale = (scale_u < scale_v) ? scale_u : scale_v;
- sep_u = as * scale * cos(theta);
- sep_v = bs * scale;
+ sep_u = (double)scale*as*sin(theta);
+ sep_v = (double)scale*as*cos(theta) + scale*bs;
+ /* We are interested in the smaller of the two separations */
max_r = (sep_u < sep_v) ? sep_u : sep_v;
+ max_r /= 2.0; /* Max radius if half the separation */
max_r -= 1.0; /* Add a tiny separation between circles */
+ if ( max_r < 1.0 ) {
+ ERROR("Circle radius is probably too small (%f).\n", max_r);
+ }
- for ( h=-INDMAX; h<INDMAX; h++ ) {
- for ( k=-INDMAX; k<INDMAX; k++ ) {
+ for ( xi=-INDMAX; xi<INDMAX; xi++ ) {
+ for ( yi=-INDMAX; yi<INDMAX; yi++ ) {
double u, v, val;
int ct;
int nequiv, p;
+ signed int h, k, l;
- ct = lookup_count(c, h, k, 0);
- if ( ct < 1 ) continue; /* Must have at least one count */
+ h = xi*xh + yi*yh;
+ k = xi*xk + yi*yk;
+ l = xi*xl + yi*yl;
+
+ /* Do we have this reflection? */
+ if ( find_unique_equiv(items, h, k, l, sym, &h, &k, &l) == 0 ) {
+ continue;
+ }
+
+
+ /* Reflection in the zone? */
+ if ( zh*h + zk*k + zl*l != 0 ) continue;
switch ( wght ) {
case WGHT_I :
- val = lookup_intensity(ref, h, k, 0);
+ val = lookup_intensity(ref, h, k, l);
break;
case WGHT_SQRTI :
- val = lookup_intensity(ref, h, k, 0);
+ val = lookup_intensity(ref, h, k, l);
val = (val>0.0) ? sqrt(val) : 0.0;
break;
case WGHT_COUNTS :
val = (float)ct;
- val /= (float)num_equivs(h, k, 0, sym);
+ val /= (float)num_equivs(h, k, l, sym);
break;
case WGHT_RAWCOUNTS :
val = (float)ct;
@@ -228,19 +287,23 @@ static void render_za(UnitCell *cell, double *ref, unsigned int *c,
abort();
}
- nequiv = num_equivs(h, k, 0, sym);
+ nequiv = num_equivs(h, k, l, sym);
for ( p=0; p<nequiv; p++ ) {
signed int he, ke, le;
+ signed int ux, uy;
float r, g, b;
- get_equiv(h, k, 0, &he, &ke, &le, sym, p);
+ get_equiv(h, k, l, &he, &ke, &le, sym, p);
+
+ ux = he*xh + ke*xk + le*xl;
+ uy = he*yh + ke*yk + le*yl;
- u = (double)he*as*sin(theta);
- v = (double)he*as*cos(theta) + ke*bs;
+ u = (double)ux*as*sin(theta);
+ v = (double)ux*as*cos(theta) + uy*bs;
- cairo_arc(dctx, ((double)wh/2)+u*scale,
- ((double)ht/2)+v*scale, max_r,
- 0, 2*M_PI);
+ cairo_arc(dctx, ((double)cx)+u*scale,
+ ((double)cy)+v*scale,
+ max_r, 0, 2*M_PI);
render_scale(val, max_val/boost, colscale, &r, &g, &b);
cairo_set_source_rgb(dctx, r, g, b);
@@ -253,39 +316,41 @@ static void render_za(UnitCell *cell, double *ref, unsigned int *c,
out:
/* Centre marker */
- cairo_arc(dctx, (double)wh/2,
- (double)ht/2, max_r, 0, 2*M_PI);
+ cairo_arc(dctx, (double)cx,
+ (double)cy, max_r, 0, 2*M_PI);
cairo_set_source_rgb(dctx, 1.0, 0.0, 0.0);
cairo_fill(dctx);
/* Draw indexing lines */
cairo_set_line_width(dctx, 4.0);
- cairo_move_to(dctx, (double)wh/2, (double)ht/2);
- u = (double)max_h*as*sin(theta);
- v = (double)max_h*as*cos(theta) + 0*bs;
- cairo_line_to(dctx, ((double)wh/2)+u*scale,
- ((double)ht/2)+v*scale);
+ cairo_move_to(dctx, (double)cx, (double)cy);
+ u = (double)max_ux*as*sin(theta);
+ v = (double)max_ux*as*cos(theta);
+ cairo_line_to(dctx, cx+u*scale, cy+v*scale);
cairo_set_source_rgb(dctx, 0.0, 1.0, 0.0);
cairo_stroke(dctx);
- cairo_move_to(dctx,((double)wh/2)+u*scale-40.0,
- ((double)ht/2)+v*scale+40.0);
- cairo_set_font_size(dctx, 40.0);
- cairo_show_text(dctx, "h");
+ cairo_set_font_size(dctx, 40.0);
+ snprintf(tmp, 255, "%i%i%i", xh, xk, xl);
+ cairo_text_extents(dctx, tmp, &size);
+
+ cairo_move_to(dctx, cx+u*scale + 20.0, cy+v*scale + size.height/2.0);
+ cairo_show_text(dctx, tmp);
cairo_fill(dctx);
- cairo_move_to(dctx, (double)wh/2, (double)ht/2);
+ snprintf(tmp, 255, "%i%i%i", yh, yk, yl);
+ cairo_text_extents(dctx, tmp, &size);
+
+ cairo_move_to(dctx, (double)cx, (double)cy);
u = 0.0;
- v = (double)max_k*bs;
- cairo_line_to(dctx, ((double)wh/2)+u*scale,
- ((double)ht/2)+v*scale);
+ v = (double)max_uy*bs;
+ cairo_line_to(dctx, cx+u*scale, cy+v*scale);
cairo_set_source_rgb(dctx, 0.0, 1.0, 0.0);
cairo_stroke(dctx);
- cairo_move_to(dctx,((double)wh/2)+u*scale-40.0,
- ((double)ht/2)+v*scale-40.0);
- cairo_set_font_size(dctx, 40.0);
- cairo_show_text(dctx, "k");
+ cairo_move_to(dctx, cx+u*scale - size.width/2.0,
+ cy+v*scale + size.height + 20.0);
+ cairo_show_text(dctx, tmp);
cairo_fill(dctx);
cairo_surface_finish(surface);
@@ -475,7 +540,6 @@ int main(int argc, char *argv[])
ref = new_list_intensity();
cts = new_list_count();
ReflItemList *items = read_reflections(infile, ref, NULL, cts);
- delete_items(items);
if ( ref == NULL ) {
ERROR("Couldn't open file '%s'\n", infile);
return 1;
@@ -485,7 +549,7 @@ int main(int argc, char *argv[])
r = povray_render_animation(cell, ref, cts, nproc);
} else if ( config_zoneaxis ) {
#ifdef HAVE_CAIRO
- render_za(cell, ref, cts, boost, sym, wght, colscale);
+ render_za(cell, items, ref, cts, boost, sym, wght, colscale);
#else
ERROR("This version of CrystFEL was compiled without Cairo");
ERROR(" support, which is required to plot a zone axis");
@@ -497,6 +561,7 @@ int main(int argc, char *argv[])
free(pdb);
free(sym);
+ delete_items(items);
return r;
}