Merge branch 'master' into simple-search

Conflicts:

	src/intensities.c
	src/reproject.c

8 files changed, 646 insertions, 420 deletions

diff --git a/data/displaywindow.ui b/data/displaywindow.ui
index 20a8a95..c59290d 100644
--- a/data/displaywindow.ui
+++ b/data/displaywindow.ui
@@ -5,6 +5,9 @@
 			<menuitem name="savecache" action="SaveCacheAction" />
 			<menuitem name="savehkl" action="SaveHKLAction" />
 			<separator />
+			<menuitem name="savecell" action="SaveCellAction" />
+			<menuitem name="loadcell" action="LoadCellAction" />
+			<separator />
 			<menuitem name="close" action="CloseAction" />
 		</menu>
 		
diff --git a/src/basis.c b/src/basis.c
index 6e7aed8..37fedba 100644
--- a/src/basis.c
+++ b/src/basis.c
@@ -3,7 +3,7 @@
  *
  * Handle basis structures
  *
- * (c) 2007 Thomas White <taw27@cam.ac.uk>
+ * (c) 2007-2009 Thomas White <taw27@cam.ac.uk>
  *
  *  dtr - Diffraction Tomography Reconstruction
  *
@@ -21,135 +21,229 @@
 #include "reflections.h"
 #include "basis.h"
 #include "utils.h"
-
-double basis_efom(ReflectionList *reflectionlist, Basis *basis) {
-
-	int n_indexed, n_counted;
-	Reflection *cur;
-	
-	cur = reflectionlist->reflections;
-	n_indexed = 0;
-	n_counted = 0;
-	while ( cur ) {
-
-		if ( cur->type == REFLECTION_NORMAL ) {
-			
-			/* Can this basis "approximately" account for this reflection? */
-			double det;
-			double a11, a12, a13, a21, a22, a23, a31, a32, a33;
-			double h, k, l;
-			
-			/* Set up the coordinate transform from hkl to xyz */
-			a11 = basis->a.x;  a12 = basis->a.y;  a13 = basis->a.z;
-			a21 = basis->b.x;  a22 = basis->b.y;  a23 = basis->b.z;
-			a31 = basis->c.x;  a32 = basis->c.y;  a33 = basis->c.z;
-			
-			/* Invert the matrix to get hkl from xyz */
-			det = a11*(a22*a33 - a23*a32) - a12*(a21*a33 - a23*a31) + a13*(a21*a32 - a22*a31);
-			h = ((a22*a33-a23*a32)*cur->x + (a23*a31-a21*a33)*cur->y + (a21*a32-a22*a31)*cur->z) / det;
-			k = ((a13*a32-a12*a33)*cur->x + (a11*a33-a13*a31)*cur->y + (a12*a31-a11*a32)*cur->z) / det;
-			l = ((a12*a23-a13*a22)*cur->x + (a13*a21-a11*a23)*cur->y + (a11*a22-a12*a21)*cur->z) / det;
-			
-			/* Calculate the deviations in terms of |a|, |b| and |c| */
-			h = fabs(h);  k = fabs(k);  l = fabs(l);
-			h -= floor(h);  k -= floor(k);  l -= floor(l);
-			if ( h == 1.0 ) h = 0.0;
-			if ( k == 1.0 ) k = 0.0;
-			if ( l == 1.0 ) l = 0.0;
-			
-			/* Define "approximately" here.  Circle in basis space becomes an ellipsoid in reciprocal space */
-			if ( h*h + k*k + l*l <= 0.1*0.1*0.1 ) n_indexed++;
-			
-			n_counted++;
-		
-		}
-		
-		cur = cur->next;
-
-	}
-	
-	return (double)n_indexed / n_counted;
-
+#include "displaywindow.h"
+
+static void basis_print(Basis *cell)
+{
+	printf("%12.8f %12.8f %12.8f\n",
+				cell->a.x/1e9, cell->a.y/1e9, cell->a.z/1e9);
+	printf("%12.8f %12.8f %12.8f\n",
+				cell->b.x/1e9, cell->b.y/1e9, cell->b.z/1e9);
+	printf("%12.8f %12.8f %12.8f\n",
+				cell->c.x/1e9, cell->c.y/1e9, cell->c.z/1e9);
 }
 
-Basis basis_add(Basis u, Basis v) {
-
-	Basis ans;
-	
-	ans.a.x = u.a.x + v.a.x;	ans.a.y = u.a.y + v.a.y;	ans.a.z = u.a.z + v.a.z;
-	ans.b.x = u.b.x + v.b.x;	ans.b.y = u.b.y + v.b.y;	ans.b.z = u.b.z + v.b.z;
-	ans.c.x = u.c.x + v.c.x;	ans.c.y = u.c.y + v.c.y;	ans.c.z = u.c.z + v.c.z;
-	
-	return ans;
-
+static void cell_print(UnitCell *cell)
+{
+	printf("%12.8f %12.8f %12.8f nm\n",
+	       cell->a*1e9, cell->b*1e9, cell->c*1e9);
+	printf("%12.8f %12.8f %12.8f deg\n",
+	       rad2deg(cell->alpha), rad2deg(cell->beta), rad2deg(cell->gamma));
 }
 
-static void basis_print(Basis *cell) {
-
-	printf("%12.8f %12.8f %12.8f\n", cell->a.x/1e9, cell->a.y/1e9, cell->a.z/1e9);
-	printf("%12.8f %12.8f %12.8f\n", cell->b.x/1e9, cell->b.y/1e9, cell->b.z/1e9);
-	printf("%12.8f %12.8f %12.8f\n", cell->c.x/1e9, cell->c.y/1e9, cell->c.z/1e9);	
-
-}
-
-static void cell_print(UnitCell *cell) {
-
-	printf("%12.8f %12.8f %12.8f nm\n", cell->a*1e9, cell->b*1e9, cell->c*1e9);
-	printf("%12.8f %12.8f %12.8f deg\n", rad2deg(cell->alpha), rad2deg(cell->beta), rad2deg(cell->gamma));
-
-}
-
-UnitCell basis_get_cell(Basis *basis) {
-
+UnitCell basis_get_cell(Basis *basis)
+{
 	UnitCell cell;
 	gsl_matrix *m;
 	gsl_matrix *inv;
 	gsl_permutation *perm;
 	double ax, ay, az, bx, by, bz, cx, cy, cz;
 	int s;
-	
-//	basis->a.x = 0.5;  basis->a.y = 0.0;  basis->a.z = 0.0;
-//	basis->b.x = 0.0;  basis->b.y = 0.2;  basis->b.z = 0.0;
-//	basis->c.x = 0.0;  basis->c.y = 0.0;  basis->c.z = 1.0;
+
 	printf("Reciprocal-space cell (nm^-1):\n");
 	basis_print(basis);
-	
+
 	m = gsl_matrix_alloc(3, 3);
-	gsl_matrix_set(m, 0, 0, basis->a.x);  gsl_matrix_set(m, 0, 1, basis->b.x);  gsl_matrix_set(m, 0, 2, basis->c.x);
-	gsl_matrix_set(m, 1, 0, basis->a.y);  gsl_matrix_set(m, 1, 1, basis->b.y);  gsl_matrix_set(m, 1, 2, basis->c.y);
-	gsl_matrix_set(m, 2, 0, basis->a.z);  gsl_matrix_set(m, 2, 1, basis->b.z);  gsl_matrix_set(m, 2, 2, basis->c.z);
-	
+	gsl_matrix_set(m, 0, 0, basis->a.x);
+	gsl_matrix_set(m, 0, 1, basis->b.x);
+	gsl_matrix_set(m, 0, 2, basis->c.x);
+	gsl_matrix_set(m, 1, 0, basis->a.y);
+	gsl_matrix_set(m, 1, 1, basis->b.y);
+	gsl_matrix_set(m, 1, 2, basis->c.y);
+	gsl_matrix_set(m, 2, 0, basis->a.z);
+	gsl_matrix_set(m, 2, 1, basis->b.z);
+	gsl_matrix_set(m, 2, 2, basis->c.z);
+
 	perm = gsl_permutation_alloc(m->size1);
 	inv = gsl_matrix_alloc(m->size1, m->size2);
 	gsl_linalg_LU_decomp(m, perm, &s);
 	gsl_linalg_LU_invert(m, perm, inv);
 	gsl_permutation_free(perm);
 	gsl_matrix_free(m);
-	
+
 	gsl_matrix_transpose(inv);
-	
-	ax = gsl_matrix_get(inv, 0, 0);  bx = gsl_matrix_get(inv, 0, 1);  cx = gsl_matrix_get(inv, 0, 2);
-	ay = gsl_matrix_get(inv, 1, 0);  by = gsl_matrix_get(inv, 1, 1);  cy = gsl_matrix_get(inv, 1, 2);
-	az = gsl_matrix_get(inv, 2, 0);  bz = gsl_matrix_get(inv, 2, 1);  cz = gsl_matrix_get(inv, 2, 2);
-	
+
+	ax = gsl_matrix_get(inv, 0, 0);
+	bx = gsl_matrix_get(inv, 0, 1);
+	cx = gsl_matrix_get(inv, 0, 2);
+	ay = gsl_matrix_get(inv, 1, 0);
+	by = gsl_matrix_get(inv, 1, 1);
+	cy = gsl_matrix_get(inv, 1, 2);
+	az = gsl_matrix_get(inv, 2, 0);
+	bz = gsl_matrix_get(inv, 2, 1);
+	cz = gsl_matrix_get(inv, 2, 2);
+
 	printf("Real-space cell (nm):\n");
 	printf("%12.8f %12.8f %12.8f\n", ax*1e9, ay*1e9, az*1e9);
 	printf("%12.8f %12.8f %12.8f\n", bx*1e9, by*1e9, bz*1e9);
 	printf("%12.8f %12.8f %12.8f\n", cx*1e9, cy*1e9, cz*1e9);
-	
+
 	cell.a = sqrt(ax*ax + ay*ay + az*az);
 	cell.b = sqrt(bx*bx + by*by + bz*bz);
 	cell.c = sqrt(cx*cx + cy*cy + cz*cz);
 	cell.alpha = acos((bx*cx + by*cy + bz*cz)/(cell.b * cell.c));
 	cell.beta = acos((ax*cx + ay*cy + az*cz)/(cell.a * cell.c));
 	cell.gamma = acos((bx*ax + by*ay + bz*az)/(cell.b * cell.a));
-	
+
 	gsl_matrix_free(inv);
-	
-	printf("Cell parameters:\n");	
+
+	printf("Cell parameters:\n");
 	cell_print(&cell);
-	
+
 	return cell;
-	
 }
 
+
+static int basis_do_save(Basis *cell, const char *filename)
+{
+	FILE *fh;
+	UnitCell rcell;
+
+	fh = fopen(filename, "w");
+
+	fprintf(fh, "# DTR unit cell description\n");
+
+	/* A "human-readable" form */
+	rcell = basis_get_cell(cell);
+	fprintf(fh, "# a %12.8f nm\n", rcell.a*1e9);
+	fprintf(fh, "# b %12.8f nm\n", rcell.b*1e9);
+	fprintf(fh, "# c %12.8f nm\n", rcell.c*1e9);
+	fprintf(fh, "# alpha %12.8f deg\n", rad2deg(rcell.alpha));
+	fprintf(fh, "# beta %12.8f deg\n", rad2deg(rcell.beta));
+	fprintf(fh, "# gamma %12.8f deg\n", rad2deg(rcell.gamma));
+
+	/* The useful form */
+	fprintf(fh, "a %12.8f %12.8f %12.8f\n",
+					cell->a.x, cell->a.y, cell->a.z);
+	fprintf(fh, "b %12.8f %12.8f %12.8f\n",
+					cell->b.x, cell->b.y, cell->b.z);
+	fprintf(fh, "c %12.8f %12.8f %12.8f\n",
+					cell->c.x, cell->c.y, cell->c.z);
+
+	fclose(fh);
+
+	return 0;
+}
+
+static gint basis_save_response(GtkWidget *widget, gint response,
+					ControlContext *ctx)
+{
+	if ( response == GTK_RESPONSE_ACCEPT ) {
+		char *filename;
+		filename = gtk_file_chooser_get_filename(
+						GTK_FILE_CHOOSER(widget));
+		if ( basis_do_save(ctx->cell, filename) ) {
+			displaywindow_error("Failed to save unit cell.",
+						ctx->dw);
+		}
+		g_free(filename);
+	}
+
+	gtk_widget_destroy(widget);
+
+	return 0;
+}
+
+void basis_save(ControlContext *ctx)
+{
+	GtkWidget *save;
+
+	save = gtk_file_chooser_dialog_new("Save Unit Cell to File",
+					GTK_WINDOW(ctx->dw->window),
+					GTK_FILE_CHOOSER_ACTION_SAVE,
+					GTK_STOCK_CANCEL, GTK_RESPONSE_CANCEL,
+					GTK_STOCK_SAVE, GTK_RESPONSE_ACCEPT,
+					NULL);
+	g_signal_connect(G_OBJECT(save), "response",
+				G_CALLBACK(basis_save_response), ctx);
+	gtk_widget_show_all(save);
+}
+
+static int basis_do_load(Basis *cell, const char *filename)
+{
+	FILE *fh;
+	float x, y, z;
+	int got_a = 0;
+	int got_b = 0;
+	int got_c = 0;
+
+	fh = fopen(filename, "r");
+
+	while ( !feof(fh) ) {
+
+		char line[256];
+
+		if ( fgets(line, 255, fh) != NULL ) {
+
+			if ( sscanf(line, "a %f %f %f\n", &x, &y, &z) == 3 ) {
+				cell->a.x = x;	cell->a.y = y;	cell->a.z = z;
+				got_a = 1;
+			}
+			if ( sscanf(line,  "b %f %f %f\n", &x, &y, &z) == 3 ) {
+				cell->b.x = x;	cell->b.y = y;	cell->b.z = z;
+				got_b = 1;
+			}
+			if ( sscanf(line, "c %f %f %f\n", &x, &y, &z) == 3 ) {
+				cell->c.x = x;	cell->c.y = y;	cell->c.z = z;
+				got_c = 1;
+			}
+
+		}
+
+	}
+
+	fclose(fh);
+
+	return !(got_a && got_b && got_c);
+}
+
+static gint basis_load_response(GtkWidget *widget, gint response,
+					ControlContext *ctx)
+{
+	if ( response == GTK_RESPONSE_ACCEPT ) {
+		char *filename;
+		filename = gtk_file_chooser_get_filename(
+						GTK_FILE_CHOOSER(widget));
+		if ( ctx->cell ) {
+			free(ctx->cell);
+		}
+		ctx->cell = malloc(sizeof(Basis));
+
+		if ( basis_do_load(ctx->cell, filename) ) {
+			displaywindow_error("Failed to load unit cell.",
+						ctx->dw);
+		} else {
+			displaywindow_update(ctx->dw);
+		}
+		g_free(filename);
+	}
+
+	gtk_widget_destroy(widget);
+
+	return 0;
+}
+
+void basis_load(ControlContext *ctx)
+{
+	GtkWidget *load;
+
+	load = gtk_file_chooser_dialog_new("Load Unit Cell from File",
+					GTK_WINDOW(ctx->dw->window),
+					GTK_FILE_CHOOSER_ACTION_OPEN,
+					GTK_STOCK_CANCEL, GTK_RESPONSE_CANCEL,
+					GTK_STOCK_SAVE, GTK_RESPONSE_ACCEPT,
+					NULL);
+	g_signal_connect(G_OBJECT(load), "response",
+				G_CALLBACK(basis_load_response), ctx);
+	gtk_widget_show_all(load);
+}
diff --git a/src/basis.h b/src/basis.h
index 25eb0ae..2253506 100644
--- a/src/basis.h
+++ b/src/basis.h
@@ -3,12 +3,12 @@
  *
  * Handle basis structures
  *
- * (c) 2007 Thomas White <taw27@cam.ac.uk>
+ * (c) 2007-2009 Thomas White <taw27@cam.ac.uk>
  *
  *  dtr - Diffraction Tomography Reconstruction
  *
  */
- 
+
 #ifndef BASIS_H
 #define BASIS_H
 
@@ -19,11 +19,11 @@
 #include "control.h"
 
 typedef struct {
-	
+
 	double x;
 	double y;
 	double z;
-	
+
 } Vector;
 
 typedef struct basis_struct {
@@ -41,9 +41,8 @@ typedef struct cell_struct {
 	double gamma;
 } UnitCell;
 
-extern double basis_efom(struct reflectionlist_struct *reflectionlist, Basis *basis);
-extern Basis basis_add(Basis u, Basis v);
 extern UnitCell basis_get_cell(Basis *cell);
+extern void basis_save(ControlContext *ctx);
+extern void basis_load(ControlContext *ctx);
 
 #endif	/* BASIS_H */
-
diff --git a/src/dirax.c b/src/dirax.c
index e58983a..5063c60 100644
--- a/src/dirax.c
+++ b/src/dirax.c
@@ -4,7 +4,7 @@
  * Invoke the DirAx auto-indexing program
  * also: handle DirAx input files
  *
- * (c) 2007 Thomas White <taw27@cam.ac.uk>
+ * (c) 2007-2009 Thomas White <taw27@cam.ac.uk>
  *
  *  dtr - Diffraction Tomography Reconstruction
  *
@@ -26,7 +26,7 @@
 #include <assert.h>
 #include <sys/ioctl.h>
 #include <termio.h>
-#include <sgtty.h>	
+#include <sgtty.h>
 
 #include "control.h"
 #include "reflections.h"
@@ -41,11 +41,11 @@ typedef enum {
 	DIRAX_INPUT_PROMPT
 } DirAxInputType;
 
-static void dirax_parseline(const char *line, ControlContext *ctx) {
-	
+static void dirax_parseline(const char *line, ControlContext *ctx)
+{
 	int i, rf;
 	char *copy;
-	
+
 	copy = strdup(line);
 	for ( i=0; i<strlen(copy); i++ ) {
 		if ( copy[i] == '\r' ) copy[i]='r';
@@ -53,15 +53,17 @@ static void dirax_parseline(const char *line, ControlContext *ctx) {
 	}
 	printf("DX: DirAx: %s\n", copy);
 	free(copy);
-	
+
 	if ( strstr(line, "reflections from file") ) {
-		displaywindow_error("DirAx can't understand this data.", ctx->dw);
+		displaywindow_error("DirAx can't understand this data.",
+								ctx->dw);
 		return;
 	}
-	
+
 	/* Is this the first line of a unit cell specification? */
 	rf = 0; i = 0;
-	while ( (i<strlen(line)) && ((line[i] == 'R') || (line[i] == 'D') || (line[i] == ' ')) ) {
+	while ( (i<strlen(line)) && ((line[i] == 'R')
+				|| (line[i] == 'D') || (line[i] == ' ')) ) {
 		if ( line[i] == 'R' ) rf = 1;
 		if ( (line[i] == 'D') && rf ) {
 			ctx->dirax_read_cell = 1;
@@ -69,52 +71,57 @@ static void dirax_parseline(const char *line, ControlContext *ctx) {
 				free(ctx->cell);
 			}
 			ctx->cell = malloc(sizeof(Basis));
-			ctx->cell->a.x = 0.0;  ctx->cell->a.y = 0.0;  ctx->cell->a.z = 0.0;
-			ctx->cell->b.x = 0.0;  ctx->cell->b.y = 0.0;  ctx->cell->b.z = 0.0;
-			ctx->cell->c.x = 0.0;  ctx->cell->c.y = 0.0;  ctx->cell->c.z = 0.0;
+			ctx->cell->a.x = 0.0;  ctx->cell->a.y = 0.0;
+						ctx->cell->a.z = 0.0;
+			ctx->cell->b.x = 0.0;  ctx->cell->b.y = 0.0;
+						ctx->cell->b.z = 0.0;
+			ctx->cell->c.x = 0.0;  ctx->cell->c.y = 0.0;
+						ctx->cell->c.z = 0.0;
 			return;
 		}
 		i++;
 	}
-	
+
 	/* Parse unit cell vectors as appropriate */
 	if ( ctx->dirax_read_cell == 1 ) {
 		/* First row of unit cell values */
 		float x1, x2, x3;
 		sscanf(line, "%f %f %f", &x1, &x2, &x3);
-		ctx->cell->a.x = x1*1e10;  ctx->cell->b.x = x2*1e10;  ctx->cell->c.x = x3*1e10;
+		ctx->cell->a.x = x1*1e10;  ctx->cell->b.x = x2*1e10;
+						ctx->cell->c.x = x3*1e10;
 		ctx->dirax_read_cell++;
 		return;
 	} else if ( ctx->dirax_read_cell == 2 ) {
 		/* First row of unit cell values */
 		float y1, y2, y3;
 		sscanf(line, "%f %f %f", &y1, &y2, &y3);
-		ctx->cell->a.y = y1*1e10;  ctx->cell->b.y = y2*1e10;  ctx->cell->c.y = y3*1e10;
+		ctx->cell->a.y = y1*1e10;  ctx->cell->b.y = y2*1e10;
+						ctx->cell->c.y = y3*1e10;
 		ctx->dirax_read_cell++;
 		return;
 	} else if ( ctx->dirax_read_cell == 3 ) {
 		/* First row of unit cell values */
 		float z1, z2, z3;
 		sscanf(line, "%f %f %f", &z1, &z2, &z3);
-		ctx->cell->a.z = z1*1e10;  ctx->cell->b.z = z2*1e10;  ctx->cell->c.z = z3*1e10;
+		ctx->cell->a.z = z1*1e10;  ctx->cell->b.z = z2*1e10;
+						ctx->cell->c.z = z3*1e10;
 		printf("DX: Read a reciprocal unit cell\n");
 		displaywindow_update(ctx->dw);
 		reproject_lattice_changed(ctx);
 		ctx->dirax_read_cell = 0;
 		return;
 	}
-	
-	ctx->dirax_read_cell = 0;
 
+	ctx->dirax_read_cell = 0;
 }
 
-static void dirax_sendline(const char *line, ControlContext *ctx) {
-	
+static void dirax_sendline(const char *line, ControlContext *ctx)
+{
 	char *copy;
 	int i;
-	
+
 	write(ctx->dirax_pty, line, strlen(line));
-	
+
 	copy = strdup(line);
 	for ( i=0; i<strlen(copy); i++ ) {
 		if ( copy[i] == '\r' ) copy[i]='\0';
@@ -122,236 +129,315 @@ static void dirax_sendline(const char *line, ControlContext *ctx) {
 	}
 	printf("DX: Sent '%s'\n", copy);	/* No newline here */
 	free(copy);
-	
 }
-/* Send a "user" command to DirAx, failing if DirAx is not idle */
-static void dirax_sendline_if_idle(const char *line, ControlContext *ctx) {
 
+/* Send a "user" command to DirAx, failing if DirAx is not idle */
+static void dirax_sendline_if_idle(const char *line, ControlContext *ctx)
+{
 	if ( ctx->dirax_step != 0 ) {
 		printf("DX: DirAx not idle\n");
 		return;
 	}
-	
-	dirax_sendline(line, ctx);
 
+	dirax_sendline(line, ctx);
 }
 
-static void dirax_send_next(ControlContext *ctx) {
-
+static void dirax_send_next(ControlContext *ctx)
+{
 	switch ( ctx->dirax_step ) {
-	
+
 		case 1 : {
 			dirax_sendline("\\echo off\n", ctx);
 			ctx->dirax_step++;
 			break;
 		}
-		
+
 		case 2 : {
 			dirax_sendline("read dtr.drx\n", ctx);
 			ctx->dirax_step++;
 			break;
 		}
-		
+
 		case 3 : {
 			dirax_sendline("dmax 10\n", ctx);
 			ctx->dirax_step++;
 			break;
 		}
-		
+
 		case 4 : {
-			dirax_sendline("indexfit 6\n", ctx);
+			dirax_sendline("indexfit 2\n", ctx);
 			ctx->dirax_step++;
 			break;
 		}
-		
+
 		case 5 : {
 			dirax_sendline("levelfit 300\n", ctx);
 			ctx->dirax_step++;
 			break;
 		}
-		
+
 		case 6 : {
 			dirax_sendline("go\n", ctx);
 			ctx->dirax_step++;
 			break;
 		}
-		
+
 		case 7 : {
 			dirax_sendline("cell\n", ctx);
 			ctx->dirax_step++;
 			break;
 		}
-				
+
 		default: {
 			ctx->dirax_step = 0;
 			printf("DX: Prompt.  DirAx is idle\n");
 		}
-		
+
 	}
-	
 }
 
-static gboolean dirax_readable(GIOChannel *dirax, GIOCondition condition, ControlContext *ctx) {
-	
+static gboolean dirax_readable(GIOChannel *dirax, GIOCondition condition,
+				ControlContext *ctx)
+{
 	int rval;
-	
-	rval = read(ctx->dirax_pty, ctx->dirax_rbuffer+ctx->dirax_rbufpos, ctx->dirax_rbuflen-ctx->dirax_rbufpos);
-	
+
+	rval = read(ctx->dirax_pty, ctx->dirax_rbuffer+ctx->dirax_rbufpos,
+			ctx->dirax_rbuflen-ctx->dirax_rbufpos);
+
 	if ( (rval == -1) || (rval == 0) ) {
-	
+
 		printf("DX: Lost connection to DirAx\n");
 		waitpid(ctx->dirax_pid, NULL, 0);
 		g_io_channel_shutdown(ctx->dirax, FALSE, NULL);
 		ctx->dirax = NULL;
 		displaywindow_update_dirax(ctx, ctx->dw);
 		return FALSE;
-		
+
 	} else {
-	
+
 		int no_string = 0;
-		
+
 		ctx->dirax_rbufpos += rval;
 		assert(ctx->dirax_rbufpos <= ctx->dirax_rbuflen);
-		
+
 		while ( (!no_string) && (ctx->dirax_rbufpos > 0) ) {
-	
+
 			int i;
 			int block_ready = 0;
 			DirAxInputType type = DIRAX_INPUT_NONE;
-			
+
 			/* See if there's a full line in the buffer yet */
-			for ( i=0; i<ctx->dirax_rbufpos-1; i++ ) {	/* Means the last value looked at is rbufpos-2 */
-			
+			for ( i=0; i<ctx->dirax_rbufpos-1; i++ ) {
+				/* Means the last value looked at is rbufpos-2 */
+
 				/* Is there a prompt in the buffer? */
 				if ( i+7 <= ctx->dirax_rbufpos ) {
-					if ( (strncmp(ctx->dirax_rbuffer+i, "Dirax> ", 7) == 0)
-					  || (strncmp(ctx->dirax_rbuffer+i, "PROMPT:", 7) == 0) ) {
+					if ( (strncmp(ctx->dirax_rbuffer+i,
+							"Dirax> ", 7) == 0)
+					  || (strncmp(ctx->dirax_rbuffer+i,
+					  		"PROMPT:", 7) == 0) ) {
 						block_ready = 1;
 						type = DIRAX_INPUT_PROMPT;
 						break;
 					}
 				}
-				
-				if ( (ctx->dirax_rbuffer[i] == '\r') && (ctx->dirax_rbuffer[i+1] == '\n') ) { 
+
+				if ( (ctx->dirax_rbuffer[i] == '\r')
+				  && (ctx->dirax_rbuffer[i+1] == '\n') ) {
 					block_ready = 1;
 					type = DIRAX_INPUT_LINE;
 					break;
 				}
-				
+
 			}
-			
+
 			if ( block_ready ) {
-			
+
 				unsigned int new_rbuflen;
 				unsigned int endbit_length;
-				
+
 				switch ( type ) {
-				
+
 					case DIRAX_INPUT_LINE : {
-					
+
 						char *block_buffer = NULL;
-						
+
 						block_buffer = malloc(i+1);
-						memcpy(block_buffer, ctx->dirax_rbuffer, i);
+						memcpy(block_buffer,
+							ctx->dirax_rbuffer, i);
 						block_buffer[i] = '\0';
-						
+
 						if ( block_buffer[0] == '\r' ) {
-							memmove(block_buffer, block_buffer+1, i);
+							memmove(block_buffer,
+							  block_buffer+1, i);
 						}
-						
-						dirax_parseline(block_buffer, ctx);
+
+						dirax_parseline(block_buffer,
+									ctx);
 						free(block_buffer);
 						endbit_length = i+2;
-						
+
 						break;
-						
+
 					}
-					
+
 					case DIRAX_INPUT_PROMPT : {
-						
+
 						dirax_send_next(ctx);
 						endbit_length = i+7;
 						break;
-					
+
 					}
-					
+
 					default : {
-						printf("DX: Unrecognised input mode (this never happens!)\n");
+						printf(
+			"DX: Unrecognised input mode (this never happens!)\n");
 						abort();
 					}
-				
+
 				}
-			
-				/* Now the block's been parsed, it should be forgotten about */
-				memmove(ctx->dirax_rbuffer, ctx->dirax_rbuffer + endbit_length, ctx->dirax_rbuflen - endbit_length);
-				
+
+				/* Now the block's been parsed, it should be
+				 * forgotten about */
+				memmove(ctx->dirax_rbuffer,
+					ctx->dirax_rbuffer + endbit_length,
+					ctx->dirax_rbuflen - endbit_length);
+
 				/* Subtract the number of bytes removed */
-				ctx->dirax_rbufpos = ctx->dirax_rbufpos - endbit_length;
-				new_rbuflen = ctx->dirax_rbuflen - endbit_length;
+				ctx->dirax_rbufpos = ctx->dirax_rbufpos
+								- endbit_length;
+				new_rbuflen = ctx->dirax_rbuflen
+								- endbit_length;
 				if ( new_rbuflen == 0 ) {
 					new_rbuflen = 256;
 				}
-				ctx->dirax_rbuffer = realloc(ctx->dirax_rbuffer, new_rbuflen);
+				ctx->dirax_rbuffer = realloc(ctx->dirax_rbuffer,
+								new_rbuflen);
 				ctx->dirax_rbuflen = new_rbuflen;
-				
+
 			} else {
-			
-				if ( ctx->dirax_rbufpos == ctx->dirax_rbuflen ) {
-			
+
+				if ( ctx->dirax_rbufpos==ctx->dirax_rbuflen ) {
+
 					/* More buffer space is needed */
-					ctx->dirax_rbuffer = realloc(ctx->dirax_rbuffer, ctx->dirax_rbuflen + 256);
-					ctx->dirax_rbuflen = ctx->dirax_rbuflen + 256;
-					/* The new space gets used at the next read, shortly... */
-			
+					ctx->dirax_rbuffer = realloc(
+						ctx->dirax_rbuffer,
+						ctx->dirax_rbuflen + 256);
+					ctx->dirax_rbuflen = ctx->dirax_rbuflen
+									+ 256;
+					/* The new space gets used at the next
+					 * read, shortly... */
+
 				}
 				no_string = 1;
-			
+
 			}
-			
+
 		}
-		
+
 	}
-	
+
 	return TRUE;
-	
 }
 
-void dirax_stop(ControlContext *ctx) {
+void dirax_stop(ControlContext *ctx)
+{
 	dirax_sendline_if_idle("end\n", ctx);
 }
 
-void dirax_rerun(ControlContext *ctx) {
+void dirax_rerun(ControlContext *ctx)
+{
 	dirax_sendline_if_idle("go\n", ctx);
 	ctx->dirax_step = 7;
 }
 
-void dirax_invoke(ControlContext *ctx) {
+static void dirax_send_random_selection(ReflectionList *r, int n, FILE *fh)
+{
+	char *used;
+	int i;
+
+	used = malloc(n*sizeof(char));
+
+	for ( i=0; i<n; i++ ) {
+		used[i] = '-';
+	}
+
+	i = 0;
+	while ( i < 1000 ) {
+
+		Reflection *ref;
+		int j;
+		long long int ra;
+
+		ra = ((long long int)random() * (long long int)n);
+		ra /= RAND_MAX;
+		if ( used[ra] == 'U' ) {
+			continue;
+		}
+
+
+		/* Dig out the correct reflection. A little faffy
+		 * because of the linked list */
+		ref = r->reflections;
+		for ( j=0; j<ra; j++ ) {
+			ref = ref->next;
+		}
+
+		/* Limit resolution of reflections */
+//		if (	( (ref->x/1e9)*(ref->x/1e9)
+//			+ (ref->y/1e9)*(ref->y/1e9)
+//			+ (ref->z/1e9)*(ref->z/1e9) ) > (20e9)*(20e9) )
+//				continue;
+
+		fprintf(fh, "%10f %10f %10f %8f\n",
+					ref->x/1e10, ref->y/1e10,
+					ref->z/1e10, ref->intensity);
+		used[ra] = 'U';
+
+		i++;
+
+	}
+}
 
+void dirax_invoke(ControlContext *ctx)
+{
 	FILE *fh;
 	Reflection *ref;
 	unsigned int opts;
 	int saved_stderr;
-	
+	int n;
+
 	if ( ctx->dirax ) {
 		dirax_rerun(ctx);
 		return;
 	}
-	
+
 	printf("DX: Starting DirAx...\n");
-	
+
 	fh = fopen("dtr.drx", "w");
 	if ( !fh ) {
 		printf("DX: Couldn't open temporary file dtr.drx\n");
 		return;
 	}
-	fprintf(fh, "%f\n", 0.5);  /* Lie about the wavelength.  DirAx can't handle the truth. */
-	ref = ctx->reflectionlist->reflections;
-	while ( ref ) {
-		fprintf(fh, "%10f %10f %10f %8f\n", ref->x/1e10, ref->y/1e10, ref->z/1e10, ref->intensity);
-		ref = ref->next;
+	fprintf(fh, "%f\n", 0.5);  /* Lie about the wavelength.  */
+
+	n = ctx->reflectionlist->n_reflections;
+	printf("DX: There are %i reflections - ", n);
+	if ( n > 1000 ) {
+		printf("sending a random selection to DirAx\n");
+		dirax_send_random_selection(ctx->reflectionlist, n, fh);
+	} else {
+		printf("sending them all to DirAx\n");
+		ref = ctx->reflectionlist->reflections;
+		while ( ref ) {
+			fprintf(fh, "%10f %10f %10f %8f\n",
+				ref->x/1e10, ref->y/1e10,
+				ref->z/1e10, ref->intensity);
+			ref = ref->next;
+		}
 	}
 	fclose(fh);
-	
+
 	saved_stderr = dup(STDERR_FILENO);
 	ctx->dirax_pid = forkpty(&ctx->dirax_pty, NULL, NULL, NULL);
 	if ( ctx->dirax_pid == -1 ) {
@@ -359,28 +445,28 @@ void dirax_invoke(ControlContext *ctx) {
 		return;
 	}
 	if ( ctx->dirax_pid == 0 ) {
-		
+
 		/* Child process: invoke DirAx */
 		struct termios t;
-		
+
 		/* Turn echo off */
 		tcgetattr(STDIN_FILENO, &t);
 		t.c_lflag &= ~(ECHO | ECHOE | ECHOK | ECHONL);
 		tcsetattr(STDIN_FILENO, TCSANOW, &t);
-		
+
 		/* Reconnect stderr */
 		dup2(saved_stderr, STDERR_FILENO);
-		
+
 		execlp("dirax", "", (char *)NULL);
 		printf("(from the Other Side) Failed to invoke DirAx.\n");
 		_exit(0);
-		
+
 	}
-		
+
 	ctx->dirax_rbuffer = malloc(256);
 	ctx->dirax_rbuflen = 256;
 	ctx->dirax_rbufpos = 0;
-	
+
 	/* Set non-blocking */
 	opts = fcntl(ctx->dirax_pty, F_GETFL);
 	fcntl(ctx->dirax_pty, F_SETFL, opts | O_NONBLOCK);
@@ -388,55 +474,59 @@ void dirax_invoke(ControlContext *ctx) {
 	ctx->dirax_step = 1;	/* This starts the "initialisation" procedure */
 	ctx->dirax_read_cell = 0;
 	ctx->dirax = g_io_channel_unix_new(ctx->dirax_pty);
-	g_io_add_watch(ctx->dirax, G_IO_IN | G_IO_HUP, (GIOFunc)dirax_readable, ctx);
-	
+	g_io_add_watch(ctx->dirax, G_IO_IN | G_IO_HUP, (GIOFunc)dirax_readable,
+									ctx);
+
 	displaywindow_update_dirax(ctx, ctx->dw);
-	
+
 	return;
-	
 }
 
-/* Despite being part of the same module, this has very little to do with invoking DirAx */
-ReflectionList *dirax_load(const char *filename) {
-
+/* Despite being part of the same module, this has very little to do with
+ * invoking DirAx */
+ReflectionList *dirax_load(const char *filename)
+{
 	FILE *fh;
 	char line[256];
 	ReflectionList *list;
 	int lambda_set = 0;
-	
+
 	fh = fopen(filename, "r");
 	if ( !fh ) {
 		printf("Couldn't open file '%s'\n", filename);
 		return 0;
 	}
-	
+
 	list = reflectionlist_new();
-	
+
 	while ( !feof(fh) ) {
-	
+
 		size_t ptr;
 		float lambda, theta, chib, phib;
-		
+
 		fgets(line, 255, fh);
 		ptr = skipspace(line);
 		if ( line[ptr] == '!' ) continue;
 		if ( line[ptr] == '\n' ) continue;
 		if ( line[ptr] == '\r' ) continue;
-		if ( sscanf(line+ptr, "%f %f %f\n", &theta, &phib, &chib) == 3 ) {
-			
+		if ( sscanf(line+ptr, "%f %f %f\n", &theta, &phib,
+								&chib) == 3 ) {
+
 			double s, x, y, z;
 			float blah, intensity;
-			
-			/* Try to find an intensity value. Use dummy value if it fails */
-			if ( sscanf(line+ptr, "%f %f %f %f\n", &blah, &blah, &blah, &intensity) != 4 ) {
+
+			/* Try to find an intensity value. Use dummy value
+			 * if it fails */
+			if ( sscanf(line+ptr, "%f %f %f %f\n", &blah, &blah,
+					&blah, &intensity) != 4 ) {
 				intensity = 1.0;
 			}
-				
+
 			if ( !lambda_set ) {
 				printf("DX: Wavelength not specified\n");
 				continue;
 			}
-			
+
 			chib = deg2rad(chib);
 			phib = deg2rad(phib);
 			theta = deg2rad(theta);
@@ -445,42 +535,43 @@ ReflectionList *dirax_load(const char *filename) {
 			y = +s*cos(chib)*cos(phib);
 			z = +s*sin(chib);
 			reflection_add(list, x, y, z, 1.0, REFLECTION_NORMAL);
-			
+
 			continue;
-			
+
 		}
 		if ( sscanf(line+ptr, "%f\n", &lambda) == 1 ) {
 			if ( lambda_set ) {
-				printf("DX: Warning: Found something which looks like a second wavelength\n");
+				printf("DX: Warning: Found something which "
+					"looks like a second wavelength\n");
 			}
 			lambda /= 1e10; /* Convert from A to m */
 			lambda_set = 1;
 		}
-		
+
 	}
-	
+
 	fclose(fh);
-	
+
 	return list;
-	
-}
 
-int dirax_is_drxfile(const char *filename) {
+}
 
+int dirax_is_drxfile(const char *filename)
+{
 	FILE *fh;
 	float lambda;
 	char line[256];
-	
+
 	fh = fopen(filename, "r");
 	if ( !fh ) {
 		printf("Couldn't open file '%s'\n", filename);
 		return 0;
 	}
-	
+
 	while ( !feof(fh) ) {
-	
+
 		size_t ptr;
-		
+
 		fgets(line, 255, fh);
 		ptr = skipspace(line);
 		if ( line[ptr] == '!' ) continue;
@@ -493,12 +584,10 @@ int dirax_is_drxfile(const char *filename) {
 		} else {
 			return 0;
 		}
-		
+
 	}
-	
+
 	fclose(fh);
-	
+
 	return 0;
-	
 }
-    
diff --git a/src/displaywindow.c b/src/displaywindow.c
index ca3edac..d0c8e3e 100644
--- a/src/displaywindow.c
+++ b/src/displaywindow.c
@@ -422,6 +422,19 @@ static gint displaywindow_savehkl(GtkWidget *widget, DisplayWindow *dw)
 	return 0;
 }
 
+static gint displaywindow_loadcell(GtkWidget *widget, DisplayWindow *dw)
+{
+	basis_load(dw->ctx);
+	return 0;
+}
+
+static gint displaywindow_savecell(GtkWidget *widget, DisplayWindow *dw)
+{
+	basis_save(dw->ctx);
+	return 0;
+}
+
+
 static void displaywindow_addmenubar(DisplayWindow *dw)
 {
 	GtkActionEntry entries[] = {
@@ -429,8 +442,12 @@ static void displaywindow_addmenubar(DisplayWindow *dw)
 	{ "FileAction", NULL, "_File", NULL, NULL, NULL },
 	{ "SaveCacheAction", "filesave", "Save Image Analysis to _Cache",
 			NULL, NULL, G_CALLBACK(displaywindow_savecache) },
-	{ "SaveHKLAction", GTK_STOCK_SAVE, "Save Reflections", NULL, NULL,
+	{ "SaveHKLAction", GTK_STOCK_SAVE, "Save Reflections...", NULL, NULL,
 					G_CALLBACK(displaywindow_savehkl) },
+	{ "LoadCellAction", GTK_STOCK_OPEN, "Load Unit Cell...", NULL, NULL,
+					G_CALLBACK(displaywindow_loadcell) },
+	{ "SaveCellAction", GTK_STOCK_SAVE, "Save Unit Cell...", NULL, NULL,
+					G_CALLBACK(displaywindow_savecell) },
 	{ "CloseAction", GTK_STOCK_QUIT, "_Quit", NULL, NULL,
 					G_CALLBACK(displaywindow_close) },
 
diff --git a/src/intensities.c b/src/intensities.c
index f1bdf19..86f3de1 100644
--- a/src/intensities.c
+++ b/src/intensities.c
@@ -3,7 +3,7 @@
  *
  * Extract integrated intensities by relrod estimation
  *
- * (c) 2007-2008 Thomas White <taw27@cam.ac.uk>
+ * (c) 2007-2009 Thomas White <taw27@cam.ac.uk>
  *
  *  dtr - Diffraction Tomography Reconstruction
  *
@@ -24,8 +24,8 @@
 /* Extract integrated reflection intensities by estimating the spike function
  * based on the observed intensity and the calculated excitation error from
  * the lattice refinement.  Easy.	*/
-void intensities_extract(ControlContext *ctx) {
-
+void intensities_extract(ControlContext *ctx)
+{
 	int i, j;
 	int n_meas, n_dupl, n_notf;
 	double max;
@@ -36,123 +36,111 @@ void intensities_extract(ControlContext *ctx) {
 		reflectionlist_free(ctx->integrated);
 	}
 	ctx->integrated = reflectionlist_new();
-	
+
 	n_meas = 0;
 	n_dupl = 0;
 	n_notf = 0;
 	max = 0;
 	for ( i=0; i<ctx->images->n_images; i++ ) {
-	
+
 		ImageRecord *image;
-		
+
 		image = &ctx->images->images[i];
-		if ( image->rflist == NULL ) image->rflist = reproject_get_reflections(image, ctx->cell_lattice, ctx);
-		
+		if ( image->rflist == NULL )
+			image->rflist = reproject_get_reflections(image,
+							ctx->cell_lattice);
+
 		for ( j=0; j<image->rflist->n_features; j++ ) {
-		
+
 			ImageFeature *feature;
 			signed int h, k, l;
-			
+
 			feature = &image->rflist->features[j];
-			
+
 			h = feature->reflection->h;
 			k = feature->reflection->k;
 			l = feature->reflection->l;
-				
+
 			if ( feature->partner != NULL ) {
-			
+
 				if ( (h!=0) || (k!=0) || (l!=0) ) {
-				
+
 					double intensity;
 					Reflection *ref;
-					
-					/* Perform relrod calculation of doom here.
-					 * TODO: Figure out if this is even possible. */
+
 					intensity = feature->partner->intensity;
-					
-					ref = reflectionlist_find(ctx->integrated, h, k, l);
-					
+
+					ref = reflectionlist_find(
+						      ctx->integrated, h, k, l);
+
 					if ( ref == NULL ) {
-					
+
 						Reflection *new;
-						
-						printf("IN: Adding %3i %3i %3i, intensity=%f\n", h, k, l, intensity);
-								
-						new = reflection_add(ctx->integrated,
-							feature->reflection->x, feature->reflection->y, feature->reflection->z,
-							intensity, REFLECTION_GENERATED);
-						
+
+						new = reflection_add(
+							ctx->integrated,
+							feature->reflection->x,
+							feature->reflection->y,
+							feature->reflection->z,
+							intensity,
+							REFLECTION_GENERATED);
+
 						new->h = h;
 						new->k = k;
 						new->l = l;
-						
-						if ( intensity > max ) max = intensity;
-						
+
+						if ( intensity > max )
+							max = intensity;
+
 						n_meas++;
-					
+
 					} else {
-						
-						printf("IN: Duplicate measurement for %3i %3i %3i - ", h, k, l);
-						
+
 						if ( intensity > ref->intensity ) {
-						
-							printf("stronger.\n");
-						
+
 							ref->x = feature->reflection->x;
 							ref->y = feature->reflection->y;
 							ref->z = feature->reflection->z;
 							ref->intensity = intensity;
-							
-						} else {
-							printf("weaker.\n");
+
 						}
-						
+
 						n_dupl++;
-						
+
 					}
-				
+
 				}
-				
+
 			} else {
 				//printf("IN: %3i %3i %3i not found\n", h, k, l);
 				n_notf++;
 			}
-		
+
 		}
-		
+
 	}
-	
+
 	/* Normalise all reflections to the most intense reflection */
 	reflection = ctx->integrated->reflections;
 	while ( reflection ) {
-	
 		reflection->intensity /= max;
-		
-		/* Test mode */
-	//	if ( (reflection->h == 0) && (reflection->k == -1) && (reflection->l == 0) ) {
-	//		reflection->intensity = 1.0;
-	//	} else {
-	//		reflection->intensity /= 10;
-	//	}
-		
 		reflection = reflection->next;
-		
 	}
-	
+
 	printf("IN: %5i intensities measured\n", n_meas);
 	printf("IN: %5i duplicated measurements\n", n_dupl);
 	printf("IN: %5i predicted reflections not found\n", n_notf);
-
 }
 
-static int intensities_do_save(ReflectionList *integrated, Basis *cell, const char *filename) {
-
+static int intensities_do_save(ReflectionList *integrated, Basis *cell,
+				const char *filename)
+{
 	FILE *fh;
 	Reflection *reflection;
 	UnitCell rcell;
-	
+
 	fh = fopen(filename, "w");
-	
+
 	rcell = basis_get_cell(cell);
 	fprintf(fh, "a %12.8f\n", rcell.a*1e9);
 	fprintf(fh, "b %12.8f\n", rcell.b*1e9);
@@ -160,46 +148,50 @@ static int intensities_do_save(ReflectionList *integrated, Basis *cell, const ch
 	fprintf(fh, "alpha %12.8f\n", rad2deg(rcell.alpha));
 	fprintf(fh, "beta %12.8f\n", rad2deg(rcell.beta));
 	fprintf(fh, "gamma %12.8f\n", rad2deg(rcell.gamma));
-	
+
 	reflection = integrated->reflections;
 	while ( reflection ) {
-		fprintf(fh, "%3i %3i %3i %12.8f\n", reflection->h, reflection->k, reflection->l, reflection->intensity);
+		fprintf(fh, "%3i %3i %3i %12.8f\n",
+			reflection->h, reflection->k, reflection->l,
+			reflection->intensity);
 		reflection = reflection->next;
 	}
 	fclose(fh);
-	
-	return 0;
 
+	return 0;
 }
 
-static gint intensities_save_response(GtkWidget *widget, gint response, ControlContext *ctx) {
-
+static gint intensities_save_response(GtkWidget *widget, gint response,
+					ControlContext *ctx)
+{
 	if ( response == GTK_RESPONSE_ACCEPT ) {
 		char *filename;
-		filename = gtk_file_chooser_get_filename(GTK_FILE_CHOOSER(widget));
-		if ( intensities_do_save(ctx->integrated, ctx->cell, filename) ) {
-			displaywindow_error("Failed to save cache file.", ctx->dw);
+		filename = gtk_file_chooser_get_filename(
+						GTK_FILE_CHOOSER(widget));
+		if ( intensities_do_save(ctx->integrated,
+						ctx->cell, filename) ) {
+			displaywindow_error("Failed to save cache file.",
+						ctx->dw);
 		}
 		g_free(filename);
 	}
-	
+
 	gtk_widget_destroy(widget);
 
 	return 0;
-	
 }
 
-void intensities_save(ControlContext *ctx) {
-
+void intensities_save(ControlContext *ctx)
+{
 	GtkWidget *save;
 
-	save = gtk_file_chooser_dialog_new("Save Reflections to File", GTK_WINDOW(ctx->dw->window),
+	save = gtk_file_chooser_dialog_new("Save Reflections to File",
+					GTK_WINDOW(ctx->dw->window),
 					GTK_FILE_CHOOSER_ACTION_SAVE,
 					GTK_STOCK_CANCEL, GTK_RESPONSE_CANCEL,
 					GTK_STOCK_SAVE, GTK_RESPONSE_ACCEPT,
 					NULL);
-	g_signal_connect(G_OBJECT(save), "response", G_CALLBACK(intensities_save_response), ctx);
+	g_signal_connect(G_OBJECT(save), "response",
+				G_CALLBACK(intensities_save_response), ctx);
 	gtk_widget_show_all(save);
-
 }
-
diff --git a/src/itrans-stat.c b/src/itrans-stat.c
index ee41817..0e03ead 100644
--- a/src/itrans-stat.c
+++ b/src/itrans-stat.c
@@ -441,7 +441,7 @@ static gsl_matrix *itrans_peaksearch_stat_floodfill(gsl_matrix *m,
 					gsl_matrix_set(p, 2, n, val);
 					n++;
 					if ( n == size ) {
-					printf("expanding %i->%i\n", size, size*2);
+
 		p = itrans_peaksearch_stat_matrix_expand(p, size, size*2);
 						size *= 2;
 					}
@@ -450,14 +450,10 @@ static gsl_matrix *itrans_peaksearch_stat_floodfill(gsl_matrix *m,
 			}
 		}
 	}
-	//printf("ff: ending loop, found %d\n",n);
 
 	*count = n;
 	if ( n > 0 ) {
-		//printf("pcheck s1=%d s2=%d\n",p->size1,p->size2);
-		printf("expandingg %i->%i\n", size, n);
 		p = itrans_peaksearch_stat_matrix_expand(p, size, n);
-		//printf("pcheck s1=%d s2=%d\n",p->size1,p->size2);
 	}
 
 	return p;
diff --git a/src/reproject.c b/src/reproject.c
index 48f60da..c41aa3f 100644
--- a/src/reproject.c
+++ b/src/reproject.c
@@ -18,114 +18,151 @@
 #include "imagedisplay.h"
 #include "displaywindow.h"
 #include "image.h"
-#include "mapping.h"
 
 /* Attempt to find partners from the feature list of "image" for each feature in "flist". */
 void reproject_partner_features(ImageFeatureList *rflist, ImageRecord *image) {
 
 	int i;
-	
+
 	for ( i=0; i<rflist->n_features; i++ ) {
-	
+
 		//if ( !reflection_is_easy(rflist->features[i].reflection) ) continue;
-		
+
 		double d = 0.0;
 		ImageFeature *partner;
 		int idx;
-		
-		partner = image_feature_closest(image->features, rflist->features[i].x, rflist->features[i].y,
-						&d, &idx);
-		
+
+		partner = image_feature_closest(image->features, rflist->features[i].x, rflist->features[i].y, &d, &idx);
+
 		if ( (d <= 20.0) && partner ) {
 			rflist->features[i].partner = partner;
 			rflist->features[i].partner_d = d;
 		} else {
 			rflist->features[i].partner = NULL;
 		}
-	
+
 	}
 
 }
 
-ImageFeatureList *reproject_get_reflections(ImageRecord *image, ReflectionList *reflectionlist, ControlContext *ctx) {
+ImageFeatureList *reproject_get_reflections(ImageRecord *image, ReflectionList *reflectionlist) {
 
 	ImageFeatureList *flist;
 	Reflection *reflection;
 	double smax = 0.1e9;
-	double tilt, omega, k;
-	double nx, ny, nz;			/* "normal" vector */
+	double tilt, omega;
+	double nx, ny, nz, nxt, nyt, nzt;	/* "normal" vector (and calculation intermediates) */
 	double kx, ky, kz;			/* Electron wavevector ("normal" times 1/lambda */
-	double ux, uy, uz;			/* "up" vector */
-	double rx, ry, rz;			/* "right" vector */
-	
+	double ux, uy, uz, uxt, uyt, uzt;	/* "up" vector (and calculation intermediates) */
+	//int done_debug = 0;
+
 	flist = image_feature_list_new();
-	
+
 	tilt = image->tilt;
 	omega = image->omega;
-	k = 1.0/image->lambda;
-	
-	/* Calculate the (normalised) incident electron wavevector */
-	mapping_rotate(0.0, 0.0, 1.0, &nx, &ny, &nz, omega, tilt);
+
+	/* Calculate the (normalised) incident electron wavevector
+	 *	n is rotated "into" the reconstruction, so only one omega step. */
+	nxt = 0.0;  nyt = 0.0;  nzt = 1.0;
+	nx = nxt;  ny = cos(tilt)*nyt + sin(tilt)*nzt;  nz = -sin(tilt)*nyt + cos(tilt)*nzt;
+	nxt = nx;  nyt = ny;  nzt = nz;
+	nx = nxt*cos(-omega) + nyt*sin(-omega);  ny = -nxt*sin(-omega) + nyt*cos(-omega);  nz = nzt;
 	kx = nx / image->lambda;
 	ky = ny / image->lambda;
 	kz = nz / image->lambda;	/* This is the centre of the Ewald sphere */
-	
-	/* Determine where "up" is */
-	mapping_rotate(0.0, 1.0, 0.0, &ux, &uy, &uz, omega, tilt);
-	
-	/* Determine where "right" is */
-	mapping_rotate(1.0, 0.0, 0.0, &rx, &ry, &rz, omega, tilt);
-	
+	//reflection_add(ctx->reflectionlist, kx, ky, kz, 1, REFLECTION_VECTOR_MARKER_1);
+
+	/* Determine where "up" is
+	 *	See above. */
+	uxt = 0.0;  uyt = 1.0;  uzt = 0.0;
+	ux = uxt;  uy = cos(tilt)*uyt + sin(tilt)*uzt;  uz = -sin(tilt)*uyt + cos(tilt)*uzt;
+	uxt = ux;  uyt = uy;  uzt = uz;
+	ux = uxt*cos(-omega) + uyt*-sin(omega);  uy = -uxt*sin(omega) + uyt*cos(omega);  uz = uzt;
+	//reflection_add(ctx->reflectionlist, ux*50e9, uy*50e9, uz*50e9, 1, REFLECTION_VECTOR_MARKER_2);
+
 	reflection = reflectionlist->reflections;
 	do {
-	
+
 		double xl, yl, zl;
 		double a, b, c;
-		double s1, s2, s, t;
-		double g_sq, gn;
-		
+		double A1, A2, s1, s2, s;
+
 		/* Get the coordinates of the reciprocal lattice point */
 		xl = reflection->x;
 		yl = reflection->y;
 		zl = reflection->z;
-		g_sq = modulus_squared(xl, yl, zl);
-		gn = xl*nx + yl*ny + zl*nz;
-		
+
 		/* Next, solve the relrod equation to calculate the excitation error */
 		a = 1.0;
-		b = 2.0*(gn - k);
-		c = -2.0*gn*k + g_sq;
-		t = -0.5*(b + sign(b)*sqrt(b*b - 4.0*a*c));
-		s1 = t/a;
-		s2 = c/t;
+		b = -2.0*(xl*nx + yl*ny + zl*nz);
+		c = xl*xl + yl*yl + zl*zl - 1.0/(image->lambda*image->lambda);
+		A1 = (-b + sqrt(b*b-4.0*a*c))/(2.0*a);
+		A2 = (-b - sqrt(b*b-4.0*a*c))/(2.0*a);
+		s1 = 1.0/image->lambda - A1;
+		s2 = 1.0/image->lambda - A2;
 		if ( fabs(s1) < fabs(s2) ) s = s1; else s = s2;
-		
+
 		/* Skip this reflection if s is large */
 		if ( fabs(s) <= smax ) {
-		
+
 			double xi, yi, zi;
 			double gx, gy, gz;
+			double cx, cy, cz;
 			double theta;
 			double x, y;
-			
+			double rx, ry, rz;
+
 			/* Determine the intersection point */
 			xi = xl + s*nx;  yi = yl + s*ny;  zi = zl + s*nz;
-			
+
 			/* Calculate Bragg angle */
 			gx = xi - kx;
 			gy = yi - ky;
 			gz = zi - kz;	/* This is the vector from the centre of the sphere to the intersection */
 			theta = angle_between(-kx, -ky, -kz, gx, gy, gz);
-			
+
 			if ( theta > 0.0 ) {
-			
-				double dx, dy, psi;
-				
-				/* Calculate azimuth of point in image (anticlockwise from +x) */
-				dx = xi*rx + yi*ry + zi*rz;
-				dy = xi*ux + yi*uy + zi*uz;
-				psi = atan2(dy, dx);
-				
+
+				double psi, disc;
+
+				//reflection_add(ctx->reflectionlist, xl, yl, zl, 1, REFLECTION_GENERATED);
+				//reflection_add(ctx->reflectionlist, xi, yi, zi, 1, REFLECTION_MARKER);
+
+				/* Calculate azimuth of point in image (clockwise from "up", will be changed later) */
+				cx = yi*nz-zi*ny;  cy = nx*zi-nz*xi;  cz = ny*xi-nx*yi;  /* c = i x n */
+				psi = angle_between(cx, cy, cz, ux, uy, uz);
+
+				/* Finally, resolve the +/- Pi ambiguity from the previous step */
+				rx = cy*nz-cz*ny;  ry = nx*cz-nz*cx;  rz = ny*cx-nx*cy;  /* r = [i x n] x n */
+				disc = angle_between(rx, ry, rz, ux, uy, uz);
+			//	if ( (i==20) && !done_debug ) {
+			//		reflection_add(ctx->reflectionlist, xi, yi, zi, 1, REFLECTION_VECTOR_MARKER_3);
+			//		reflection_add(ctx->reflectionlist, cx, cy, cz, 1, REFLECTION_VECTOR_MARKER_4);
+			//		reflection_add(ctx->reflectionlist, rx, ry, rz, 1, REFLECTION_VECTOR_MARKER_4);
+			//		printf("psi=%f deg, disc=%f deg\n", rad2deg(psi), rad2deg(disc));
+			//	}
+				if ( (psi >= M_PI_2) && (disc >= M_PI_2) ) {
+					psi -= M_PI_2;		/* Case 1 */
+			//		if ( (i==20) && !done_debug ) printf("case 1\n");
+				} else if ( (psi >= M_PI_2) && (disc < M_PI_2) ) {
+					psi = 3*M_PI_2 - psi;	/* Case 2 */
+			//		if ( (i==20) && !done_debug ) printf("case 2\n");
+				} else if ( (psi < M_PI_2) && (disc < M_PI_2) ) {
+					psi = 3*M_PI_2 - psi;	/* Case 3 */
+			//		if ( (i==20) && !done_debug ) printf("case 3\n");
+				} else if ( (psi < M_PI_2) && (disc >= M_PI_2) ) {
+					psi = 3*M_PI_2 + psi;	/* Case 4 */
+			//		if ( (i==20) && !done_debug ) printf("case 4\n");
+				}
+
+			//	if ( (i==20) && !done_debug ) printf("final psi=%f clockwise from 'up'\n", rad2deg(psi));
+				psi = M_PI_2 - psi;	/* Anticlockwise from "+x" instead of clockwise from "up" */
+			//	if ( (i==20) && !done_debug ) printf("final psi=%f anticlockwise from +x\n", rad2deg(psi));
+
+				psi += omega;
+			//	if ( (i==20) && !done_debug ) printf("final psi=%f anticlockwise from +tilt axis\n", rad2deg(psi));
+			//	if ( (i==20) && !done_debug ) done_debug = 1;
+
 				/* Calculate image coordinates from polar representation */
 				if ( image->fmode == FORMULATION_CLEN ) {
 					x = image->camera_len*tan(theta)*cos(psi);
@@ -138,39 +175,39 @@ ImageFeatureList *reproject_get_reflections(ImageRecord *image, ReflectionList *
 					x /= image->pixel_size;
 					y /= image->pixel_size;
 				} else {
-					fprintf(stderr, "Unrecognised formulation mode in reproject_get_reflections\n");
+					fprintf(stderr, "Unrecognised formulation mode in reproject_get_reflections()\n");
 					return NULL;
 				}
-				
+
 				x += image->x_centre;
 				y += image->y_centre;
-				
+
 				/* Sanity check */
 				if ( (x>=0) && (x<image->width) && (y>=0) && (y<image->height) ) {
-				
-					/* Record the reflection
-					 *  Intensity should be multiplied by relrod spike function except
-					 *  reprojected reflections aren't used quantitatively for anything
-					 */
-					image_add_feature_reflection(flist, x, y, image,
-								     reflection->intensity, reflection);
-					
+
+					/* Record the reflection */
+					image_add_feature_reflection(flist, x, y, image, reflection->intensity, reflection);
+					/* Intensity should be multiplied by relrod spike function except
+					 * reprojected reflections aren't used quantitatively for anything */
+
+					//printf("Reflection at %f, %f\n", x, y);
+
 				} /* else it's outside the picture somewhere */
-				
+
 			} /* else this is the central beam so don't worry about it */
-		
+
 		}
-		
+
 		reflection = reflection->next;
-	
+
 	} while ( reflection );
-	
+
 	/* Partner features only if the image has a feature list.  This allows the test
 	 *	program to use this function to generate simulated data. */
 	if ( image->features != NULL ) {
 		reproject_partner_features(flist, image);
 	}
-	
+
 	return flist;
 
 }
@@ -179,21 +216,21 @@ ImageFeatureList *reproject_get_reflections(ImageRecord *image, ReflectionList *
 void reproject_cell_to_lattice(ControlContext *ctx) {
 
 	int i;
-	
+
 	if ( ctx->cell_lattice ) {
 		reflection_list_from_new_cell(ctx->cell_lattice, ctx->cell);
 	} else {
 		ctx->cell_lattice = reflection_list_from_cell(ctx->cell);
 	}
-	
+
 	displaywindow_enable_cell_functions(ctx->dw, TRUE);
-	
+
 	/* Invalidate all the reprojected feature lists */
 	for ( i=0; i<ctx->images->n_images; i++ ) {
 		image_feature_list_free(ctx->images->images[i].rflist);
 		ctx->images->images[i].rflist = NULL;
 	}
-	
+
 }
 
 /* Notify that ctx->cell has changed (also need to call displaywindow_update) */
@@ -203,4 +240,3 @@ void reproject_lattice_changed(ControlContext *ctx) {
 	displaywindow_update_imagestack(ctx->dw);
 
 }
-