aboutsummaryrefslogtreecommitdiff
path: root/src
diff options
context:
space:
mode:
authorThomas White <taw@physics.org>2011-08-02 15:09:34 +0200
committerThomas White <taw@physics.org>2012-02-22 15:27:35 +0100
commitee25b4c412c3759158eae3153170cb7940aee143 (patch)
tree13545931778c8f780e773e0f5e8af44d4dd7cddd /src
parent6a30e47ab458cd89f8aa3ed9b48af719cb936988 (diff)
Simplify render_hkl --zone-axis
Diffstat (limited to 'src')
-rw-r--r--src/reflist-utils.c61
-rw-r--r--src/reflist-utils.h5
-rw-r--r--src/render_hkl.c87
3 files changed, 95 insertions, 58 deletions
diff --git a/src/reflist-utils.c b/src/reflist-utils.c
index af43b1ff..d6272217 100644
--- a/src/reflist-utils.c
+++ b/src/reflist-utils.c
@@ -414,3 +414,64 @@ RefList *asymmetric_indices(RefList *in, const SymOpList *sym)
return new;
}
+
+
+/**
+ * resolution_limits:
+ * @list: A %RefList
+ * @cell: A %UnitCell
+ * @rmin: Place to store the minimum 1/d value
+ * @rmax: Place to store the maximum 1/d value
+ *
+ * This function calculates the minimum and maximum values of 1/d, where
+ * 2dsin(theta) = wavelength. The answers are in m^-1.
+ **/
+void resolution_limits(RefList *list, UnitCell *cell,
+ double *rmin, double *rmax)
+{
+ Reflection *refl;
+ RefListIterator *iter;
+
+ *rmin = INFINITY;
+ *rmax = 0.0;
+
+ for ( refl = first_refl(list, &iter);
+ refl != NULL;
+ refl = next_refl(refl, iter) )
+ {
+ double r;
+ signed int h, k, l;
+
+ get_indices(refl, &h, &k, &l);
+ r = 2.0 * resolution(cell, h, k, l);
+
+ if ( r > *rmax ) *rmax = r;
+ if ( r < *rmin ) *rmin = r;
+ }
+}
+
+
+/**
+ * max_intensity:
+ * @list: A %RefList
+ *
+ * Returns: The maximum intensity in @list.
+ **/
+double max_intensity(RefList *list)
+{
+ Reflection *refl;
+ RefListIterator *iter;
+ double max;
+
+ max = -INFINITY;
+
+ for ( refl = first_refl(list, &iter);
+ refl != NULL;
+ refl = next_refl(refl, iter) )
+ {
+ double val = get_intensity(refl);
+ if ( val > max ) max = val;
+ }
+
+ return max;
+}
diff --git a/src/reflist-utils.h b/src/reflist-utils.h
index 59aa6172..d14e5f8e 100644
--- a/src/reflist-utils.h
+++ b/src/reflist-utils.h
@@ -44,4 +44,9 @@ extern int find_equiv_in_list(RefList *list, signed int h, signed int k,
extern RefList *asymmetric_indices(RefList *in, const SymOpList *sym);
+extern void resolution_limits(RefList *list, UnitCell *cell,
+ double *rmin, double *rmax);
+
+extern double max_intensity(RefList *list);
+
#endif /* REFLIST_UTILS_H */
diff --git a/src/render_hkl.c b/src/render_hkl.c
index 18d6bb81..50450a94 100644
--- a/src/render_hkl.c
+++ b/src/render_hkl.c
@@ -82,27 +82,19 @@ static void show_help(const char *s)
#ifdef HAVE_CAIRO
-static void draw_circles(signed int xh, signed int xk, signed int xl,
- signed int yh, signed int yk, signed int yl,
+static void draw_circles(double xh, double xk, double xl,
+ double yh, double yk, double yl,
signed int zh, signed int zk, signed int zl,
RefList *list, const SymOpList *sym,
cairo_t *dctx, int wght, double boost, int colscale,
UnitCell *cell, double radius, double theta,
double as, double bs, double cx, double cy,
- double scale,
- signed int *max_ux, signed int *max_uy,
- double *max_val, double *max_u, double *max_v,
- double *max_res)
+ double scale, double *max_val)
{
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 */
@@ -110,9 +102,8 @@ static void draw_circles(signed int xh, signed int xk, signed int xl,
refl != NULL;
refl = next_refl(refl, iter) ) {
- double u, v, val, res;
+ double u, v, val;
signed int ha, ka, la;
- int xi, yi;
int i, n;
get_indices(refl, &ha, &ka, &la);
@@ -123,6 +114,7 @@ static void draw_circles(signed int xh, signed int xk, signed int xl,
for ( i=0; i<n; i++ ) {
signed int h, k, l;
+ double xi, yi;
get_equiv(sym, m, i, ha, ka, la, &h, &k, &l);
@@ -171,10 +163,6 @@ static void draw_circles(signed int xh, signed int xk, signed int xl,
} else {
- /* Find max vectors in plane for scaling */
- if ( fabs(u) > fabs(*max_u) ) *max_u = fabs(u);
- if ( fabs(v) > fabs(*max_v) ) *max_v = fabs(v);
-
/* Find max value for colour scale */
if ( !isnan(val) && !isinf(val)
&& (fabs(val) > fabs(*max_val)) )
@@ -182,16 +170,6 @@ static void draw_circles(signed int xh, signed int xk, signed int xl,
*max_val = fabs(val);
}
- /* Find max indices */
- if ( (yi==0) && (fabs(xi) > *max_ux) )
- *max_ux = fabs(xi);
- if ( (xi==0) && (fabs(yi) > *max_uy) )
- *max_uy = fabs(yi);
-
- /* Find max resolution */
- res = resolution(cell, h, k, l);
- if ( res > *max_res ) *max_res = res;
-
}
}
@@ -261,11 +239,10 @@ static void render_za(UnitCell *cell, RefList *list,
{
cairo_surface_t *surface;
cairo_t *dctx;
- double max_u, max_v, max_res, max_val;
- double scale_u, scale_v, scale;
+ double max_val;
+ double scale1, scale2, scale;
double sep_u, sep_v, max_r;
double u, v;
- signed int max_ux, max_uy;
double as, bs, theta;
double asx, asy, asz;
double bsx, bsy, bsz;
@@ -279,6 +256,7 @@ static void render_za(UnitCell *cell, RefList *list,
double cx, cy;
const double border = 200.0;
int png;
+ double rmin, rmax;
/* Vector product to determine the zone axis. */
zh = xk*yl - xl*yk;
@@ -304,19 +282,16 @@ static void render_za(UnitCell *cell, RefList *list,
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;
+ as = modulus(xx, xy, xz);
+ bs = modulus(yx, yy, yz);
- scale = 1.0;
- draw_circles(xh, xk, xl, yh, yk, yl, zh, zk, zl,
- list, sym, NULL, wght, boost, colscale, cell,
- 0.0, theta, as, bs, 0.0, 0.0, scale,
- &max_ux, &max_uy, &max_val, &max_u, &max_v, &max_res);
+ resolution_limits(list, cell, &rmin, &rmax);
+ printf("Resolution limits: 1/d = %.2f - %.2f nm^-1"
+ " (d = %.2f - %.2f A)\n",
+ rmin/1e9, rmax/1e9, (1.0/rmin)/1e-10, (1.0/rmax)/1e-10);
- max_res /= 1e9;
- printf("Maximum resolution is 1/d = %5.3f nm^-1, d = %5.3f nm\n",
- max_res*2.0, 1.0/(max_res*2.0));
+ max_val = max_intensity(list);
if ( max_val <= 0.0 ) {
STATUS("Couldn't find max value.\n");
return;
@@ -327,23 +302,18 @@ static void render_za(UnitCell *cell, RefList *list,
max_val = scale_top;
}
- /* Choose whichever scaling factor gives the smallest value */
- 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;
+ scale1 = ((double)wh-border) / (2.0*rmax);
+ scale2 = ((double)ht-border) / (2.0*rmax);
+ scale = (scale1 < scale2) ? scale1 : scale2;
+ /* Work out the spot radius */
sep_u = scale*as;
sep_v = 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 is 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);
- }
-
- if ( outfile == NULL ) outfile = "za.pdf";
+ /* Create surface */
if ( strcmp(outfile+strlen(outfile)-4, ".png") == 0 ) {
png = 1;
surface = cairo_image_surface_create(CAIRO_FORMAT_ARGB32,
@@ -382,7 +352,7 @@ static void render_za(UnitCell *cell, RefList *list,
draw_circles(xh, xk, xl, yh, yk, yl, zh, zk, zl,
list, sym, dctx, wght, boost, colscale, cell,
max_r, theta, as, bs, cx, cy, scale,
- NULL, NULL, &max_val, NULL, NULL, NULL);
+ &max_val);
/* Centre marker */
cairo_arc(dctx, (double)cx,
@@ -394,8 +364,8 @@ static void render_za(UnitCell *cell, RefList *list,
cairo_set_line_cap(dctx, CAIRO_LINE_CAP_ROUND);
cairo_set_line_width(dctx, 4.0);
cairo_move_to(dctx, (double)cx, (double)cy);
- u = (2.0+max_ux)*as*sin(theta);
- v = (2.0+max_ux)*as*cos(theta);
+ u = rmax*sin(theta);
+ v = rmax*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);
@@ -412,14 +382,12 @@ static void render_za(UnitCell *cell, RefList *list,
cairo_text_extents(dctx, tmp, &size);
cairo_move_to(dctx, (double)cx, (double)cy);
- u = 0.0;
- v = (2.0+max_uy)*bs;
- cairo_line_to(dctx, cx+u*scale, cy+v*scale);
+ cairo_line_to(dctx, cx, cy+rmax*scale);
cairo_set_source_rgb(dctx, 0.0, 1.0, 0.0);
cairo_stroke(dctx);
- cairo_move_to(dctx, cx+u*scale - size.width/2.0,
- cy+v*scale + size.height + 20.0);
+ cairo_move_to(dctx, cx - size.width/2.0,
+ cy+rmax*scale + size.height + 20.0);
render_overlined_indices(dctx, yh, yk, yl);
cairo_fill(dctx);
@@ -684,6 +652,9 @@ int main(int argc, char *argv[])
weighting = strdup("I");
}
+ if ( outfile == NULL ) outfile = "za.pdf";
+
+
if ( strcmp(weighting, "I") == 0 ) {
wght = WGHT_I;
} else if ( strcmp(weighting, "sqrtI") == 0 ) {