Simplify render_hkl --zone-axis

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 ) {