Merge branch 'master' of ssh://git-weiss@jade.msm.cam.ac.uk/srv/git/dtr

Conflicts:

	src/dirax.c

6 files changed, 339 insertions, 262 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 79b20ea..6e97d21 100644
--- a/src/dirax.c
+++ b/src/dirax.c
@@ -35,14 +35,12 @@
 #include "displaywindow.h"
 #include "reproject.h"
 
-
 typedef enum {
 	DIRAX_INPUT_NONE,
 	DIRAX_INPUT_LINE,
 	DIRAX_INPUT_PROMPT
 } DirAxInputType;
 
-
 static void dirax_parseline(const char *line, ControlContext *ctx)
 {
 	int i, rf;
@@ -58,14 +56,14 @@ static void dirax_parseline(const char *line, ControlContext *ctx)
 
 	if ( strstr(line, "reflections from file") ) {
 		displaywindow_error("DirAx can't understand this data.",
-					ctx->dw);
+								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] == ' ')) ) {
+				|| (line[i] == 'D') || (line[i] == ' ')) ) {
 		if ( line[i] == 'R' ) rf = 1;
 		if ( (line[i] == 'D') && rf ) {
 			ctx->dirax_read_cell = 1;
@@ -73,15 +71,12 @@ 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++;
@@ -92,27 +87,24 @@ static void dirax_parseline(const char *line, ControlContext *ctx)
 		/* 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);
@@ -123,7 +115,6 @@ static void dirax_parseline(const char *line, ControlContext *ctx)
 	ctx->dirax_read_cell = 0;
 }
 
-
 static void dirax_sendline(const char *line, ControlContext *ctx)
 {
 	char *copy;
@@ -140,7 +131,6 @@ static void dirax_sendline(const char *line, ControlContext *ctx)
 	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)
 {
@@ -152,7 +142,6 @@ static void dirax_sendline_if_idle(const char *line, ControlContext *ctx)
 	dirax_sendline(line, ctx);
 }
 
-
 static void dirax_send_next(ControlContext *ctx)
 {
 	switch ( ctx->dirax_step ) {
@@ -182,7 +171,7 @@ static void dirax_send_next(ControlContext *ctx)
 		}
 
 		case 5 : {
-			dirax_sendline("levelfit 200\n", ctx);
+			dirax_sendline("levelfit 300\n", ctx);
 			ctx->dirax_step++;
 			break;
 		}
@@ -207,7 +196,6 @@ static void dirax_send_next(ControlContext *ctx)
 	}
 }
 
-
 static gboolean dirax_readable(GIOChannel *dirax, GIOCondition condition,
 				ControlContext *ctx)
 {
@@ -240,7 +228,7 @@ static gboolean dirax_readable(GIOChannel *dirax, GIOCondition condition,
 
 			/* 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 */
+				/* Means the last value looked at is rbufpos-2 */
 
 				/* Is there a prompt in the buffer? */
 				if ( i+7 <= ctx->dirax_rbufpos ) {
@@ -255,7 +243,7 @@ static gboolean dirax_readable(GIOChannel *dirax, GIOCondition condition,
 				}
 
 				if ( (ctx->dirax_rbuffer[i] == '\r')
-					&& (ctx->dirax_rbuffer[i+1] == '\n') ) {
+				  && (ctx->dirax_rbuffer[i+1] == '\n') ) {
 					block_ready = 1;
 					type = DIRAX_INPUT_LINE;
 					break;
@@ -281,11 +269,11 @@ static gboolean dirax_readable(GIOChannel *dirax, GIOCondition condition,
 
 						if ( block_buffer[0] == '\r' ) {
 							memmove(block_buffer,
-							 block_buffer+1, i);
+							  block_buffer+1, i);
 						}
 
 						dirax_parseline(block_buffer,
-								ctx);
+									ctx);
 						free(block_buffer);
 						endbit_length = i+2;
 
@@ -302,28 +290,29 @@ static gboolean dirax_readable(GIOChannel *dirax, GIOCondition condition,
 					}
 
 					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;
+								- 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);
+								new_rbuflen);
 				ctx->dirax_rbuflen = new_rbuflen;
 
 			} else {
@@ -332,10 +321,10 @@ static gboolean dirax_readable(GIOChannel *dirax, GIOCondition condition,
 
 					/* More buffer space is needed */
 					ctx->dirax_rbuffer = realloc(
-					ctx->dirax_rbuffer,
-					ctx->dirax_rbuflen + 256);
+						ctx->dirax_rbuffer,
+						ctx->dirax_rbuflen + 256);
 					ctx->dirax_rbuflen = ctx->dirax_rbuflen
-							+ 256;
+									+ 256;
 					/* The new space gets used at the next
 					 * read, shortly... */
 
@@ -349,23 +338,19 @@ static gboolean dirax_readable(GIOChannel *dirax, GIOCondition condition,
 	}
 
 	return TRUE;
-
 }
 
-
 void dirax_stop(ControlContext *ctx)
 {
 	dirax_sendline_if_idle("end\n", ctx);
 }
 
-
 void dirax_rerun(ControlContext *ctx)
 {
 	dirax_sendline_if_idle("go\n", ctx);
 	ctx->dirax_step = 7;
 }
 
-
 static void dirax_send_random_selection(ReflectionList *r, int n, FILE *fh)
 {
 	char *used;
@@ -414,14 +399,12 @@ static void dirax_send_random_selection(ReflectionList *r, int n, FILE *fh)
 	}
 }
 
-
 void dirax_invoke(ControlContext *ctx)
 {
 	FILE *fh;
 	Reflection *ref;
 	unsigned int opts;
 	int saved_stderr;
-	int n;
 
 	if ( ctx->dirax ) {
 		dirax_rerun(ctx);
@@ -435,23 +418,12 @@ void dirax_invoke(ControlContext *ctx)
 		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. */
-
-	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;
-		}
+	fprintf(fh, "%f\n", 0.5);  /* Lie about the wavelength.  */
+	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);
 
@@ -492,11 +464,12 @@ void dirax_invoke(ControlContext *ctx)
 	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);
+									ctx);
 
 	displaywindow_update_dirax(ctx, ctx->dw);
-}
 
+	return;
+}
 
 /* Despite being part of the same module, this has very little to do with
  * invoking DirAx */
@@ -525,14 +498,14 @@ ReflectionList *dirax_load(const char *filename)
 		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 */
+			/* 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;
@@ -569,12 +542,11 @@ ReflectionList *dirax_load(const char *filename)
 	fclose(fh);
 
 	return list;
-}
 
+}
 
 int dirax_is_drxfile(const char *filename)
 {
-
 	FILE *fh;
 	float lambda;
 	char line[256];
diff --git a/src/displaywindow.c b/src/displaywindow.c
index 6488b34..1ffd95e 100644
--- a/src/displaywindow.c
+++ b/src/displaywindow.c
@@ -414,6 +414,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[] = {
@@ -421,8 +434,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 a2e19c3..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);
-		
+		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);
-
 }
-