Merge branch 'master' into simple-search

9 files changed, 662 insertions, 533 deletions

diff --git a/src/dirax.c b/src/dirax.c
index 5063c60..b1d421e 100644
--- a/src/dirax.c
+++ b/src/dirax.c
@@ -171,7 +171,7 @@ static void dirax_send_next(ControlContext *ctx)
 		}
 
 		case 5 : {
-			dirax_sendline("levelfit 300\n", ctx);
+			dirax_sendline("levelfit 200\n", ctx);
 			ctx->dirax_step++;
 			break;
 		}
diff --git a/src/glbits.c b/src/glbits.c
index ba03adf..86b710e 100644
--- a/src/glbits.c
+++ b/src/glbits.c
@@ -35,7 +35,7 @@ static GLuint glbits_load_shader(const char *filename, GLenum type) {
 	FILE *fh;
 	int l;
 	GLint status;
-	
+
 	fh = fopen(filename, "r");
 	if ( fh == NULL ) {
 		fprintf(stderr, "Couldn't load shader '%s'\n", filename);
@@ -48,7 +48,7 @@ static GLuint glbits_load_shader(const char *filename, GLenum type) {
 	shader = glCreateShader(type);
 	glShaderSource(shader, 1, &source, NULL);
 	glCompileShader(shader);
-	glGetShaderiv(shader, GL_COMPILE_STATUS, &status); 
+	glGetShaderiv(shader, GL_COMPILE_STATUS, &status);
 	if ( status == GL_FALSE ) {
 		glGetShaderInfoLog(shader, 4095, &l, text);
 		if ( l > 0 ) {
@@ -57,13 +57,13 @@ static GLuint glbits_load_shader(const char *filename, GLenum type) {
 			printf("Shader compilation failed.\n");
 		}
 	}
-	
+
 	return shader;
-	
+
 }
 
 static void glbits_load_shaders(DisplayWindow *dw) {
-	
+
 	/* Lighting-per-fragment */
 	dw->gl_vshader_lightpp = glbits_load_shader(DATADIR"/dtr/light-pp.vert", GL_VERTEX_SHADER);
 	dw->gl_fshader_lightpp = glbits_load_shader(DATADIR"/dtr/light-pp.frag", GL_FRAGMENT_SHADER);
@@ -71,7 +71,7 @@ static void glbits_load_shaders(DisplayWindow *dw) {
 	glAttachShader(dw->gl_program_lightpp, dw->gl_vshader_lightpp);
 	glAttachShader(dw->gl_program_lightpp, dw->gl_fshader_lightpp);
 	glLinkProgram(dw->gl_program_lightpp);
-	
+
 }
 
 static void glbits_delete_shaders(DisplayWindow *dw) {
@@ -93,7 +93,7 @@ static void glbits_delete_shaders(DisplayWindow *dw) {
 	normals[3*i + 1] = yv;					\
 	normals[3*i + 2] = zv;					\
 	i++;
-	
+
 #define DRAW_BLOB							\
 	double step = M_PI/(double)BLOB_BITS;				\
 	int is, js;							\
@@ -121,13 +121,13 @@ static void glbits_delete_shaders(DisplayWindow *dw) {
 			ADD_VERTEX					\
 		}							\
 	}
-	
+
 #define DRAW_POINTER_LINE												\
 	glBegin(GL_LINES);												\
 		glVertex3f(1.0, 0.0, 0.0);										\
 		glVertex3f(0.0, 0.0, 0.0);										\
 	glEnd();
-	
+
 #define DRAW_POINTER_HEAD												\
 	glPushMatrix();													\
 	for ( pointer_head_face = 1; pointer_head_face <= 4; pointer_head_face++ ) {					\
@@ -162,9 +162,9 @@ void glbits_prepare(DisplayWindow *dw) {
 	ControlContext *ctx;
 	GLfloat *vertices;
 	GLfloat *normals;
-	
+
 	ctx = dw->ctx;
-	
+
 	/* "Measured" reflections */
 	if ( dw->gl_use_buffers ) {
 		glGenBuffers(1, &dw->gl_ref_vertex_buffer);
@@ -193,7 +193,7 @@ void glbits_prepare(DisplayWindow *dw) {
 				i++;
 			}
 			reflection = reflection->next;
-		};	
+		};
 		if ( dw->gl_use_buffers ) {
 			glBindBuffer(GL_ARRAY_BUFFER, dw->gl_ref_vertex_buffer);
 			glBufferData(GL_ARRAY_BUFFER, 3*dw->gl_ref_num_vertices*sizeof(GLfloat), vertices, GL_STATIC_DRAW);
@@ -206,7 +206,7 @@ void glbits_prepare(DisplayWindow *dw) {
 			dw->gl_ref_normal_array = normals;
 		}
 	}
-	
+
 	/* Marker "reflections" */
 	if ( dw->gl_use_buffers ) {
 		glGenBuffers(1, &dw->gl_marker_vertex_buffer);
@@ -230,7 +230,7 @@ void glbits_prepare(DisplayWindow *dw) {
 				DRAW_BLOB
 			}
 			reflection = reflection->next;
-		};	
+		};
 		if ( dw->gl_use_buffers ) {
 			glBindBuffer(GL_ARRAY_BUFFER, dw->gl_marker_vertex_buffer);
 			glBufferData(GL_ARRAY_BUFFER, 3*dw->gl_marker_num_vertices*sizeof(GLfloat), vertices, GL_STATIC_DRAW);
@@ -243,7 +243,7 @@ void glbits_prepare(DisplayWindow *dw) {
 			dw->gl_marker_normal_array = normals;
 		}
 	}
-	
+
 	/* Generated reflections */
 	if ( dw->gl_use_buffers ) {
 		glGenBuffers(1, &dw->gl_gen_vertex_buffer);
@@ -266,7 +266,7 @@ void glbits_prepare(DisplayWindow *dw) {
 				double size = 5.0 * log(1+0.1*reflection->intensity);
 				DRAW_BLOB
 				reflection = reflection->next;
-			};	
+			};
 			if ( dw->gl_use_buffers ) {
 				glBindBuffer(GL_ARRAY_BUFFER, dw->gl_gen_vertex_buffer);
 				glBufferData(GL_ARRAY_BUFFER, 3*dw->gl_gen_num_vertices*sizeof(GLfloat), vertices, GL_STATIC_DRAW);
@@ -283,25 +283,25 @@ void glbits_prepare(DisplayWindow *dw) {
 	} else {
 		dw->gl_gen_num_vertices = 0;
 	}
-	
+
 	/* Indexing lines */
 	glLineWidth(2.0);
 	if ( ctx->cell && dw->lines ) {
-		
+
 		int max_ind;
 		signed int h, k, l;
-		
+
 		max_ind = 10;
 		dw->gl_line_num_vertices = 3*2*((2*1+1)*(2*1+1));
-		
+
 		if ( dw->gl_use_buffers ) {
 			glGenBuffers(1, &dw->gl_line_vertex_buffer);
 		}
 		reflection = ctx->reflectionlist->reflections;
 		vertices = malloc(3*dw->gl_line_num_vertices*sizeof(GLfloat));
-		
+
 		i=0;
-		
+
 		/* Lines parallel to a */
 		for ( k=-1; k<=1; k++ ) {
 			for ( l=-1; l<=1; l++ ) {
@@ -341,7 +341,7 @@ void glbits_prepare(DisplayWindow *dw) {
 				i++;
 			}
 		}
-		
+
 		if ( dw->gl_use_buffers ) {
 			glBindBuffer(GL_ARRAY_BUFFER, dw->gl_line_vertex_buffer);
 			glBufferData(GL_ARRAY_BUFFER, 3*dw->gl_line_num_vertices*sizeof(GLfloat), vertices, GL_STATIC_DRAW);
@@ -349,12 +349,12 @@ void glbits_prepare(DisplayWindow *dw) {
 		} else {
 			dw->gl_line_vertex_array = vertices;
 		}
-		
+
 	}
-	
+
 	dw->gl_list_id = glGenLists(1);
 	glNewList(dw->gl_list_id, GL_COMPILE);
-	
+
 	#if 0
 	GLUquadric *quad;
 	quad = gluNewQuadric();
@@ -364,7 +364,7 @@ void glbits_prepare(DisplayWindow *dw) {
 	glMaterialf(GL_FRONT_AND_BACK, GL_SHININESS, 10.0);
 	gluSphere(quad, 10, 32, 32);
 	#endif
-	
+
 	/* Bounding cube: 100 nm^-1 side length */
 	if ( dw->cube ) {
 		glMaterialfv(GL_FRONT_AND_BACK, GL_AMBIENT_AND_DIFFUSE, black);
@@ -380,7 +380,7 @@ void glbits_prepare(DisplayWindow *dw) {
 			glVertex3f(-50.0, -50.0, 50.0);
 			glNormal3f(50.0, -50.0, 50.0);
 			glVertex3f(50.0, -50.0, 50.0);
-			
+
 		glEnd();
 		glBegin(GL_LINE_LOOP);
 			glNormal3f(50.0, 50.0, -50.0);
@@ -411,7 +411,7 @@ void glbits_prepare(DisplayWindow *dw) {
 			glVertex3f(50.0, -50.0, -50.0);
 		glEnd();
 	}
-	
+
 	/* x, y, z pointers */
 	int pointer_head_face;
 	glPushMatrix();
@@ -427,7 +427,7 @@ void glbits_prepare(DisplayWindow *dw) {
 	glMaterialfv(GL_FRONT_AND_BACK, GL_AMBIENT_AND_DIFFUSE, red);
 	DRAW_POINTER_HEAD
 	glPopMatrix();
-	
+
 	glPushMatrix();
 	glRotatef(90.0, 0.0, 0.0, 1.0);
 	glScalef(10.0, 1.0, 1.0);
@@ -442,7 +442,7 @@ void glbits_prepare(DisplayWindow *dw) {
 	glMaterialfv(GL_FRONT_AND_BACK, GL_AMBIENT_AND_DIFFUSE, green);
 	DRAW_POINTER_HEAD
 	glPopMatrix();
-	
+
 	glPushMatrix();
 	glRotatef(-90.0, 0.0, 1.0, 0.0);
 	glScalef(10.0, 1.0, 1.0);
@@ -457,13 +457,13 @@ void glbits_prepare(DisplayWindow *dw) {
 	glMaterialfv(GL_FRONT_AND_BACK, GL_AMBIENT_AND_DIFFUSE, bblue);
 	DRAW_POINTER_HEAD
 	glPopMatrix();
-	
+
 	/* Plot the other reflections */
 	reflection = ctx->reflectionlist->reflections;
 	while ( reflection != NULL ) {
-	
+
 		if ( reflection->type == REFLECTION_VECTOR_MARKER_1 ) {
-		
+
 			glMaterialfv(GL_FRONT_AND_BACK, GL_EMISSION, black);
 			glMaterialfv(GL_FRONT_AND_BACK, GL_AMBIENT_AND_DIFFUSE, red);
 			glMaterialfv(GL_FRONT_AND_BACK, GL_SPECULAR, black);
@@ -474,11 +474,11 @@ void glbits_prepare(DisplayWindow *dw) {
 				glNormal3f(reflection->x/1e9, reflection->y/1e9, reflection->z/1e9);
 				glVertex3f(reflection->x/1e9, reflection->y/1e9, reflection->z/1e9);
 			glEnd();
-			
+
 		}
-		
+
 		if ( reflection->type == REFLECTION_VECTOR_MARKER_2 ) {
-		
+
 			glMaterialfv(GL_FRONT_AND_BACK, GL_EMISSION, black);
 			glMaterialfv(GL_FRONT_AND_BACK, GL_AMBIENT_AND_DIFFUSE, green);
 			glMaterialfv(GL_FRONT_AND_BACK, GL_SPECULAR, black);
@@ -489,11 +489,11 @@ void glbits_prepare(DisplayWindow *dw) {
 				glNormal3f(reflection->x/1e9, reflection->y/1e9, reflection->z/1e9);
 				glVertex3f(reflection->x/1e9, reflection->y/1e9, reflection->z/1e9);
 			glEnd();
-			
+
 		}
-		
+
 		if ( reflection->type == REFLECTION_VECTOR_MARKER_3 ) {
-		
+
 			glMaterialfv(GL_FRONT_AND_BACK, GL_EMISSION, black);
 			glMaterialfv(GL_FRONT_AND_BACK, GL_AMBIENT_AND_DIFFUSE, blue);
 			glMaterialfv(GL_FRONT_AND_BACK, GL_SPECULAR, black);
@@ -504,89 +504,89 @@ void glbits_prepare(DisplayWindow *dw) {
 				glNormal3f(reflection->x/1e9, reflection->y/1e9, reflection->z/1e9);
 				glVertex3f(reflection->x/1e9, reflection->y/1e9, reflection->z/1e9);
 			glEnd();
-			
+
 		}
-		
+
 		reflection = reflection->next;
-		
+
 	};
-	
+
 	/* Draw the reciprocal unit cell if one is available */
 	if ( ctx->cell && !dw->lines ) {
-	
+
 		glBegin(GL_LINES);
-			
+
 			glMaterialfv(GL_FRONT_AND_BACK, GL_AMBIENT_AND_DIFFUSE, black);
 			glMaterialfv(GL_FRONT_AND_BACK, GL_SPECULAR, black);
 			glMaterialf(GL_FRONT_AND_BACK, GL_SHININESS, 0.0);
 			glNormal3f(1.0, 0.0, 0.0);
-			
+
 			glMaterialfv(GL_FRONT_AND_BACK, GL_EMISSION, red);
 			glVertex3f(0.0, 0.0, 0.0);
 			glVertex3f(ctx->cell->a.x/1e9, ctx->cell->a.y/1e9, ctx->cell->a.z/1e9);
-			
+
 			glMaterialfv(GL_FRONT_AND_BACK, GL_EMISSION, green);
 			glVertex3f(0.0, 0.0, 0.0);
 			glVertex3f(ctx->cell->b.x/1e9, ctx->cell->b.y/1e9, ctx->cell->b.z/1e9);
-			
+
 			glMaterialfv(GL_FRONT_AND_BACK, GL_EMISSION, bblue);
 			glVertex3f(0.0, 0.0, 0.0);
 			glVertex3f(ctx->cell->c.x/1e9, ctx->cell->c.y/1e9, ctx->cell->c.z/1e9);
-			
+
 			glMaterialfv(GL_FRONT_AND_BACK, GL_EMISSION, purple);
 			glVertex3f(ctx->cell->a.x/1e9, ctx->cell->a.y/1e9, ctx->cell->a.z/1e9);
 			glVertex3f(ctx->cell->a.x/1e9 + ctx->cell->b.x/1e9,
 					ctx->cell->a.y/1e9 + ctx->cell->b.y/1e9,
 					ctx->cell->a.z/1e9 + ctx->cell->b.z/1e9);
-			
+
 			glVertex3f(ctx->cell->a.x/1e9 + ctx->cell->b.x/1e9,
 					ctx->cell->a.y/1e9 + ctx->cell->b.y/1e9,
 					ctx->cell->a.z/1e9 + ctx->cell->b.z/1e9);
 			glVertex3f(ctx->cell->b.x/1e9, ctx->cell->b.y/1e9, ctx->cell->b.z/1e9);
-			
+
 			glVertex3f(ctx->cell->c.x/1e9, ctx->cell->c.y/1e9, ctx->cell->c.z/1e9);
 			glVertex3f(ctx->cell->c.x/1e9 + ctx->cell->a.x/1e9,
 					ctx->cell->c.y/1e9 + ctx->cell->a.y/1e9,
 					ctx->cell->c.z/1e9 + ctx->cell->a.z/1e9);
-			
+
 			glVertex3f(ctx->cell->c.x/1e9 + ctx->cell->a.x/1e9,
 					ctx->cell->c.y/1e9 + ctx->cell->a.y/1e9,
 					ctx->cell->c.z/1e9 + ctx->cell->a.z/1e9);
 			glVertex3f(ctx->cell->c.x/1e9 + ctx->cell->a.x/1e9 + ctx->cell->b.x/1e9,
 					ctx->cell->c.y/1e9 + ctx->cell->a.y/1e9 + ctx->cell->b.y/1e9,
 					ctx->cell->c.z/1e9 + ctx->cell->a.z/1e9 + ctx->cell->b.z/1e9);
-			
+
 			glVertex3f(ctx->cell->c.x/1e9 + ctx->cell->a.x/1e9 + ctx->cell->b.x/1e9,
 					ctx->cell->c.y/1e9 + ctx->cell->a.y/1e9 + ctx->cell->b.y/1e9,
 					ctx->cell->c.z/1e9 + ctx->cell->a.z/1e9 + ctx->cell->b.z/1e9);
 			glVertex3f(ctx->cell->c.x/1e9 + ctx->cell->b.x/1e9,
 					ctx->cell->c.y/1e9 + ctx->cell->b.y/1e9,
 					ctx->cell->c.z/1e9 + ctx->cell->b.z/1e9);
-			
+
 			glVertex3f(ctx->cell->c.x/1e9 + ctx->cell->b.x/1e9,
 					ctx->cell->c.y/1e9 + ctx->cell->b.y/1e9,
 					ctx->cell->c.z/1e9 + ctx->cell->b.z/1e9);
 			glVertex3f(ctx->cell->c.x/1e9, ctx->cell->c.y/1e9, ctx->cell->c.z/1e9);
-			
+
 			glVertex3f(ctx->cell->c.x/1e9 + ctx->cell->b.x/1e9,
 					ctx->cell->c.y/1e9 + ctx->cell->b.y/1e9,
 					ctx->cell->c.z/1e9 + ctx->cell->b.z/1e9);
 			glVertex3f(ctx->cell->b.x/1e9, ctx->cell->b.y/1e9, ctx->cell->b.z/1e9);
-			
+
 			glVertex3f(ctx->cell->c.x/1e9 + ctx->cell->a.x/1e9 + ctx->cell->b.x/1e9,
 					ctx->cell->c.y/1e9 + ctx->cell->a.y/1e9 + ctx->cell->b.y/1e9,
 					ctx->cell->c.z/1e9 + ctx->cell->a.z/1e9 + ctx->cell->b.z/1e9);
 			glVertex3f(ctx->cell->a.x/1e9 + ctx->cell->b.x/1e9,
 					ctx->cell->a.y/1e9 + ctx->cell->b.y/1e9,
 					ctx->cell->a.z/1e9 + ctx->cell->b.z/1e9);
-			
+
 			glVertex3f(ctx->cell->c.x/1e9 + ctx->cell->a.x/1e9,
 					ctx->cell->c.y/1e9 + ctx->cell->a.y/1e9,
 					ctx->cell->c.z/1e9 + ctx->cell->a.z/1e9);
 			glVertex3f(ctx->cell->a.x/1e9, ctx->cell->a.y/1e9, ctx->cell->a.z/1e9);
-						
+
 		glEnd();
-		
+
 	}
 
 	/* Tilt axis */
@@ -611,22 +611,22 @@ void glbits_prepare(DisplayWindow *dw) {
 		DRAW_POINTER_HEAD
 		glPopMatrix();
 	}
-	
+
 	/* Zero plane (must be drawn last for transparency to work) */
-	glBegin(GL_QUADS);
-		glMaterialfv(GL_FRONT_AND_BACK, GL_EMISSION, black);
-		glMaterialfv(GL_FRONT_AND_BACK, GL_AMBIENT_AND_DIFFUSE, glass);
-		glMaterialfv(GL_FRONT_AND_BACK, GL_SPECULAR, black);
-		glMaterialf(GL_FRONT_AND_BACK, GL_SHININESS, 0.0);
-		glNormal3f(0.0, 0.0, 1.0);
-		glVertex3f(50.0, 50.0, 0.0);
-		glVertex3f(50.0, -50.0, 0.0);
-		glVertex3f(-50.0, -50.0, 0.0);
-		glVertex3f(-50.0, 50.0, 0.0);
-	glEnd();
-	
+//	glBegin(GL_QUADS);
+//		glMaterialfv(GL_FRONT_AND_BACK, GL_EMISSION, black);
+//		glMaterialfv(GL_FRONT_AND_BACK, GL_AMBIENT_AND_DIFFUSE, glass);
+//		glMaterialfv(GL_FRONT_AND_BACK, GL_SPECULAR, black);
+//		glMaterialf(GL_FRONT_AND_BACK, GL_SHININESS, 0.0);
+//		glNormal3f(0.0, 0.0, 1.0);
+//		glVertex3f(50.0, 50.0, 0.0);
+//		glVertex3f(50.0, -50.0, 0.0);
+//		glVertex3f(-50.0, -50.0, 0.0);
+//		glVertex3f(-50.0, 50.0, 0.0);
+//	glEnd();
+
 	glEndList();
-	
+
 	//printf("DW: Vertex counts: meas:%i, mark:%i, gen:%i\n", dw->gl_ref_num_vertices, dw->gl_marker_num_vertices, dw->gl_gen_num_vertices);
 
 }
@@ -634,7 +634,7 @@ void glbits_prepare(DisplayWindow *dw) {
 void glbits_set_ortho(DisplayWindow *dw, GLfloat w, GLfloat h) {
 
 	GLfloat aspect = w/h;
-	
+
 	glMatrixMode(GL_PROJECTION);
 	glLoadIdentity();
 	glOrtho(-aspect*(dw->distance/2.0), aspect*(dw->distance/2.0), -(dw->distance/2.0), (dw->distance/2.0), 0.001, 400.0);
@@ -651,8 +651,8 @@ void glbits_set_perspective(DisplayWindow *dw, GLfloat w, GLfloat h) {
 
 }
 
-gint glbits_expose(GtkWidget *widget, GdkEventExpose *event, DisplayWindow *dw) {
-	
+gint glbits_expose(GtkWidget *widget, GdkEventExpose *event, DisplayWindow *dw)
+{
 	GdkGLContext *glcontext = gtk_widget_get_gl_context(widget);
 	GdkGLDrawable *gldrawable = gtk_widget_get_gl_drawable(widget);
 	float m[4][4];
@@ -664,57 +664,57 @@ gint glbits_expose(GtkWidget *widget, GdkEventExpose *event, DisplayWindow *dw)
 	GLfloat light0_diffuse[]  = { 0.8, 0.8, 0.8, 1.0 };
 	GLfloat light0_specular[] = { 0.8, 0.8, 0.8, 1.0 };
 	GLfloat grey[] = { 0.6, 0.6, 0.6, 1.0 };
-	
+	GLfloat bg_top[] = { 1.0, 1.0, 1.0, 1.0 };
+	GLfloat bg_bot[] = { 1.0, 1.0, 1.0, 1.0 };
+	GLfloat w = dw->drawing_area->allocation.width;
+	GLfloat h = dw->drawing_area->allocation.height;
+	GLfloat aspect = w/h;
+
 	if ( !gdk_gl_drawable_gl_begin(gldrawable, glcontext) ) {
 		return 0;
 	}
-	
+
 	glClear(GL_COLOR_BUFFER_BIT | GL_DEPTH_BUFFER_BIT);
-	
+
 	if ( dw->background ) {
-		
-		GLfloat w = dw->drawing_area->allocation.width;
-		GLfloat h = dw->drawing_area->allocation.height;
-		GLfloat aspect = w/h;
-		GLfloat bg_top[] = { 0.0, 0.3, 1.0, 1.0 };
-		GLfloat bg_bot[] = { 0.0, 0.7, 1.0, 1.0 };
-		
-		glMatrixMode(GL_MODELVIEW);
-		glLoadIdentity();
-		glTranslatef(0.0, 0.0, -20.0);
-		/* Draw the background (this is done before setting up rotations) */
-		/* Set up "private" projection matrix */
-		glMatrixMode(GL_PROJECTION);
-		glPushMatrix();
-		glLoadIdentity();
-		glOrtho(-aspect*3.0, aspect*3.0, -3.0, 3.0, 0.001, 21.0);
-		/* Draw background plane */
-		glMaterialfv(GL_FRONT, GL_AMBIENT_AND_DIFFUSE, black);
-		glMaterialfv(GL_FRONT, GL_SPECULAR, black);
-		glMaterialf(GL_FRONT, GL_SHININESS, 0.0);
-		glBegin(GL_QUADS);
-			glMaterialfv(GL_FRONT_AND_BACK, GL_EMISSION, bg_bot);
-			glVertex3f(-3.0*aspect, -3.0, 0.0);
-			glVertex3f(+3.0*aspect, -3.0, 0.0);
-			glMaterialfv(GL_FRONT_AND_BACK, GL_EMISSION, bg_top);
-			glVertex3f(+3.0*aspect, +3.0, 0.0);
-			glVertex3f(-3.0*aspect, +3.0, 0.0);
-		glEnd();
-		/* Restore the old projection matrix */
-		glPopMatrix();
-		glClear(GL_DEPTH_BUFFER_BIT);	/* Background does not count for depth test purposes */
-		
+		bg_top[0] = 0.0; bg_top[1] = 0.3; bg_top[2] = 1.0;
+		bg_bot[0] = 0.0; bg_bot[1] = 0.7; bg_bot[2] = 1.0;
 	}
-	
+
+	glMatrixMode(GL_MODELVIEW);
+	glLoadIdentity();
+	glTranslatef(0.0, 0.0, -20.0);
+	/* Draw the background (this is done before setting up rotations) */
+	/* Set up "private" projection matrix */
+	glMatrixMode(GL_PROJECTION);
+	glPushMatrix();
+	glLoadIdentity();
+	glOrtho(-aspect*3.0, aspect*3.0, -3.0, 3.0, 0.001, 21.0);
+	/* Draw background plane */
+	glMaterialfv(GL_FRONT, GL_AMBIENT_AND_DIFFUSE, black);
+	glMaterialfv(GL_FRONT, GL_SPECULAR, black);
+	glMaterialf(GL_FRONT, GL_SHININESS, 0.0);
+	glBegin(GL_QUADS);
+		glMaterialfv(GL_FRONT_AND_BACK, GL_EMISSION, bg_bot);
+		glVertex3f(-3.0*aspect, -3.0, 0.0);
+		glVertex3f(+3.0*aspect, -3.0, 0.0);
+		glMaterialfv(GL_FRONT_AND_BACK, GL_EMISSION, bg_top);
+		glVertex3f(+3.0*aspect, +3.0, 0.0);
+		glVertex3f(-3.0*aspect, +3.0, 0.0);
+	glEnd();
+	/* Restore the old projection matrix */
+	glPopMatrix();
+	glClear(GL_DEPTH_BUFFER_BIT);
+
 	glMatrixMode(GL_MODELVIEW);
 	glLoadIdentity();
-	
+
 	/* Set up lighting */
 	glEnable(GL_LIGHT0);
 	glLightfv(GL_LIGHT0, GL_POSITION, light0_position);
 	glLightfv(GL_LIGHT0, GL_DIFFUSE, light0_diffuse);
 	glLightfv(GL_LIGHT0, GL_SPECULAR, light0_specular);
-	
+
 	/* The z component of this makes no difference if the projection is orthographic,
 	 *	but controls zoom when perspective is used */
 	glTranslatef(0.0, 0.0, -400.0);
@@ -726,24 +726,24 @@ gint glbits_expose(GtkWidget *widget, GdkEventExpose *event, DisplayWindow *dw)
 	glPushClientAttrib(GL_CLIENT_VERTEX_ARRAY_BIT);
 	glEnableClientState(GL_VERTEX_ARRAY);
 	glEnableClientState(GL_NORMAL_ARRAY);
-			
+
 	if ( dw->mode == DW_MAPPED ) {
-		
+
 		/* Draw the "measured" reflections */
 		if ( dw->gl_ref_num_vertices ) {
-			
+
 			GLfloat att[] = {1.0, 1.0, 0.0};
-			
+
 			glPointParameterfv(GL_POINT_DISTANCE_ATTENUATION, att);
 			glPointSize(20.0);
 			glEnable(GL_POINT_SMOOTH);
-			
+
 			glMaterialfv(GL_FRONT_AND_BACK, GL_EMISSION, black);
 			glMaterialfv(GL_FRONT_AND_BACK, GL_AMBIENT_AND_DIFFUSE, black);
 			glMaterialfv(GL_FRONT_AND_BACK, GL_SPECULAR, black);
 			glColor3f(0.0, 0.0, 0.0);
 			glMaterialf(GL_FRONT_AND_BACK, GL_SHININESS, 0.0);
-			
+
 			if ( dw->gl_use_buffers ) {
 				glBindBuffer(GL_ARRAY_BUFFER, dw->gl_ref_vertex_buffer);
 				glVertexPointer(3, GL_FLOAT, 0, NULL);
@@ -756,20 +756,20 @@ gint glbits_expose(GtkWidget *widget, GdkEventExpose *event, DisplayWindow *dw)
 				glNormalPointer(GL_FLOAT, 0, dw->gl_ref_normal_array);
 				glDrawArrays(GL_POINTS, 0, dw->gl_ref_num_vertices);
 			}
-			
+
 			glDisable(GL_POINT_SMOOTH);
-			
+
 		}
-		
+
 		/* Draw marker points */
 		if ( dw->gl_marker_num_vertices ) {
-		
+
 			glMaterialfv(GL_FRONT_AND_BACK, GL_AMBIENT_AND_DIFFUSE, blue);
 			glMaterialfv(GL_FRONT_AND_BACK, GL_EMISSION, black);
 			glMaterialfv(GL_FRONT_AND_BACK, GL_SPECULAR, white);
 			glMaterialf(GL_FRONT_AND_BACK, GL_SHININESS, 100.0);
 			glColor3f(0.0, 0.0, 1.0);
-			
+
 			if ( dw->gl_use_buffers ) {
 				glBindBuffer(GL_ARRAY_BUFFER, dw->gl_marker_vertex_buffer);
 				glVertexPointer(3, GL_FLOAT, 0, NULL);
@@ -782,20 +782,20 @@ gint glbits_expose(GtkWidget *widget, GdkEventExpose *event, DisplayWindow *dw)
 				glNormalPointer(GL_FLOAT, 0, dw->gl_marker_normal_array);
 				glDrawArrays(GL_QUADS, 0, dw->gl_marker_num_vertices);
 			}
-			
+
 		}
-		
+
 	} else {
 
 		/* Draw generated reflections */
 		if ( dw->gl_gen_num_vertices ) {
-			
+
 			glMaterialfv(GL_FRONT_AND_BACK, GL_EMISSION, black);
 			glMaterialfv(GL_FRONT_AND_BACK, GL_AMBIENT_AND_DIFFUSE, gold);
 			glMaterialfv(GL_FRONT_AND_BACK, GL_SPECULAR, white);
 			glMaterialf(GL_FRONT_AND_BACK, GL_SHININESS, 100.0);
 			glColor3f(0.7, 0.7, 0.0);
-			
+
 			if ( dw->gl_use_shaders ) glUseProgram(dw->gl_program_lightpp);
 			if ( dw->gl_use_buffers ) {
 				glBindBuffer(GL_ARRAY_BUFFER, dw->gl_gen_vertex_buffer);
@@ -810,19 +810,19 @@ gint glbits_expose(GtkWidget *widget, GdkEventExpose *event, DisplayWindow *dw)
 				glDrawArrays(GL_QUADS, 0, dw->gl_gen_num_vertices);
 			}
 			if ( dw->gl_use_shaders ) glUseProgram(0);
-			
+
 		}
-	
+
 	}
 	glDisable(GL_NORMAL_ARRAY);
-	
+
 	/* Draw indexing lines */
 	if ( dw->lines && dw->gl_line_num_vertices ) {
-	
+
 		glMaterialfv(GL_FRONT_AND_BACK, GL_EMISSION, grey);
 		glMaterialfv(GL_FRONT_AND_BACK, GL_AMBIENT_AND_DIFFUSE, black);
 		glMaterialfv(GL_FRONT_AND_BACK, GL_SPECULAR, black);
-		
+
 		if ( dw->gl_use_buffers ) {
 			glBindBuffer(GL_ARRAY_BUFFER, dw->gl_line_vertex_buffer);
 			glVertexPointer(3, GL_FLOAT, 0, NULL);
@@ -833,18 +833,18 @@ gint glbits_expose(GtkWidget *widget, GdkEventExpose *event, DisplayWindow *dw)
 			glDrawArrays(GL_LINES, 0, dw->gl_line_num_vertices);
 		}
 	}
-	
+
 	glPopClientAttrib();
-	
+
 	/* Draw everything else */
 	glCallList(dw->gl_list_id);
-	
+
 	if ( gdk_gl_drawable_is_double_buffered(gldrawable) ) {
 		gdk_gl_drawable_swap_buffers(gldrawable);
 	} else {
 		glFlush();
 	}
-	
+
 	gdk_gl_drawable_gl_end(gldrawable);
 
 	return TRUE;
@@ -857,29 +857,29 @@ gboolean glbits_configure(GtkWidget *widget, GdkEventConfigure *event, DisplayWi
 	GdkGLDrawable *gldrawable = gtk_widget_get_gl_drawable(widget);
 	GLfloat w = widget->allocation.width;
 	GLfloat h = widget->allocation.height;
-	
+
 	/* Set viewport */
 	if ( !gdk_gl_drawable_gl_begin(gldrawable, glcontext) ) {
 		return FALSE;
 	}
 	glViewport(0, 0, w, h);
-	
+
 	glEnable(GL_LIGHTING);
 	glEnable(GL_DEPTH_TEST);
 	glEnable(GL_BLEND);
 	glBlendFunc(GL_SRC_ALPHA, GL_ONE_MINUS_SRC_ALPHA);
 	glShadeModel(GL_SMOOTH);
 	glEnable(GL_LINE_SMOOTH);
-	
+
 	/* Nudge the projection matrix routines to preserve the aspect ratio */
 	if ( dw->view == DW_ORTHO ) {
 		glbits_set_ortho(dw, w, h);
 	} else {
 		glbits_set_perspective(dw, w, h);
 	}
-	
+
 	gdk_gl_drawable_gl_end(gldrawable);
-	
+
 	return FALSE;
 
 }
@@ -904,7 +904,7 @@ void glbits_free_resources(DisplayWindow *dw) {
 		if ( dw->gl_line_vertex_array != NULL ) free(dw->gl_line_vertex_array);
 	}
 	glDeleteLists(dw->gl_list_id, 1);
-	
+
 }
 
 void glbits_final_free_resources(DisplayWindow *dw) {
@@ -923,22 +923,21 @@ gint glbits_realise(GtkWidget *widget, DisplayWindow *dw) {
 	GdkGLDrawable *gldrawable = gtk_widget_get_gl_drawable(widget);
 	GLfloat w = widget->allocation.width;
 	GLfloat h = widget->allocation.height;
-	
+
 	if ( !gdk_gl_drawable_gl_begin(gldrawable, glcontext) ) {
 		return 0;
 	}
-	
+
 	/* This has to be done once an OpenGL context has been created */
 	GLenum glew_err = glewInit();
 	if (glew_err != GLEW_OK) fprintf(stderr, "GLEW initialisation error: %s\n", glewGetErrorString(glew_err));
-		
+
 	glbits_set_ortho(dw, w, h);
 	glbits_first_prepare(dw);
 	dw->realised = 1;
-	
+
 	gdk_gl_drawable_gl_end(gldrawable);
 
 	return 0;
 
 }
-
diff --git a/src/imagedisplay.c b/src/imagedisplay.c
index 8915c41..958622f 100644
--- a/src/imagedisplay.c
+++ b/src/imagedisplay.c
@@ -34,10 +34,10 @@ static void imagedisplay_rescale(ImageDisplay *imagedisplay, unsigned int v_w, u
 
 	unsigned int w, h;
 	float aspect_image, aspect_window;
-	
+
 	w = imagedisplay->imagerecord.width;
 	h = imagedisplay->imagerecord.height;
-	
+
 	/* Preserve aspect ratio */
 	aspect_image = (float)w/h;
 	aspect_window = (float)v_w/v_h;
@@ -46,16 +46,16 @@ static void imagedisplay_rescale(ImageDisplay *imagedisplay, unsigned int v_w, u
 	} else {
 		v_h = v_w/aspect_image;
 	}
-	
+
 	if ( imagedisplay->pixbuf_scaled ) {
 		g_object_unref(imagedisplay->pixbuf_scaled);
 	}
-	
+
 	/* Create the scaled pixbuf from the 8-bit display data */
 	imagedisplay->pixbuf_scaled = gdk_pixbuf_scale_simple(imagedisplay->pixbuf, v_w, v_h, GDK_INTERP_BILINEAR);
 	imagedisplay->view_width = v_w;
 	imagedisplay->view_height = v_h;
-	
+
 }
 
 static gboolean imagedisplay_configure_event(GtkWidget *widget, GdkEventConfigure *event, ImageDisplay *imagedisplay) {
@@ -63,7 +63,7 @@ static gboolean imagedisplay_configure_event(GtkWidget *widget, GdkEventConfigur
 	imagedisplay->drawingarea_width = event->width;
 	imagedisplay->drawingarea_height = event->height;
 	imagedisplay_rescale(imagedisplay, event->width, event->height);
-	
+
 	return FALSE;
 
 }
@@ -75,14 +75,14 @@ void imagedisplay_put_data(ImageDisplay *imagedisplay, ImageRecord imagerecord)
 	int min, max;
 	double c, scale;
 	int b;
-	
+
 	h = imagerecord.height;
 	w = imagerecord.width;
-	
+
 	if ( imagedisplay->pixbuf ) {
 		g_object_unref(imagedisplay->pixbuf);
 	}
-	
+
 	min = 2<<15;  max = 0;
 	b = (h+w)/20;	/* Number of pixels of ignored border */
 	for ( y=b; y<h-b; y++ ) {
@@ -93,33 +93,33 @@ void imagedisplay_put_data(ImageDisplay *imagedisplay, ImageRecord imagerecord)
 			if ( val < min ) min = val;
 		}
 	}
-	
+
 	c = 0.1;
 	scale = 255.0 / log(1+c*(max-min));
-	
+
 	/* Turn 16-bit image data into 8-bit display data */
 	imagedisplay->data = malloc(3*w*h);
 	for ( y=0; y<h; y++ ) {
 		for ( x=0; x<w; x++ ) {
-		
+
 			uint16_t val16, val8;
-			
+
 			val16 = imagerecord.image[x+w*y];
 			val8 = scale * log(1+c*(val16-min));
-			
+
 			imagedisplay->data[3*( x+w*(h-1-y) )] = val8;
 			imagedisplay->data[3*( x+w*(h-1-y) )+1] = val8;
 			imagedisplay->data[3*( x+w*(h-1-y) )+2] = val8;
-			
+
 		}
 	}
-	
+
 	memcpy(&imagedisplay->imagerecord, &imagerecord, sizeof(ImageRecord));
-	
+
 	/* Create the pixbuf from the 8-bit display data */
 	imagedisplay->pixbuf = gdk_pixbuf_new_from_data(imagedisplay->data, GDK_COLORSPACE_RGB, FALSE, 8, w, h, w*3,
 				(GdkPixbufDestroyNotify)imagedisplay_free_data, imagedisplay);
-	
+
 	if ( imagedisplay->realised ) {
 		imagedisplay_force_redraw(imagedisplay);
 	}
@@ -127,9 +127,9 @@ void imagedisplay_put_data(ImageDisplay *imagedisplay, ImageRecord imagerecord)
 }
 
 void imagedisplay_clear_marks(ImageDisplay *imagedisplay) {
-	
+
 	ImageDisplayMark *cur;
-	
+
 	cur = imagedisplay->marks;
 	while ( cur ) {
 		ImageDisplayMark *next = cur->next;
@@ -137,28 +137,28 @@ void imagedisplay_clear_marks(ImageDisplay *imagedisplay) {
 		cur = next;
 	}
 	imagedisplay->marks = NULL;
-	
+
 }
 
 static void imagedisplay_destroyed(GtkWidget *widget, ImageDisplay *imagedisplay) {
 
 	imagedisplay_clear_marks(imagedisplay);
-	
+
 	if ( imagedisplay->flags & IMAGEDISPLAY_QUIT_IF_CLOSED ) {
 		gtk_exit(0);
 	}
-	
+
 	g_object_unref(G_OBJECT(imagedisplay->gc_centre));
 	g_object_unref(G_OBJECT(imagedisplay->gc_tiltaxis));
 	g_object_unref(G_OBJECT(imagedisplay->gc_marks_1));
 	g_object_unref(G_OBJECT(imagedisplay->gc_marks_2));
-	
+
 	if ( imagedisplay->flags & IMAGEDISPLAY_FREE ) {
 		free(imagedisplay->imagerecord.image);
 	}
 
 	free(imagedisplay);
-	
+
 }
 
 void imagedisplay_close(ImageDisplay *imagedisplay) {
@@ -176,11 +176,11 @@ static void imagedisplay_add_scalebar(ImageDisplay *imagedisplay, GtkWidget *dra
 	PangoRectangle rect;
 	int bwidth, bheight;
 	int view_height = imagedisplay->view_height;
-	
+
 	sb = mapping_scale_bar_length(&imagedisplay->imagerecord);
 	layout = gtk_widget_create_pango_layout(drawingarea, "1 nm^-1");
 	pango_layout_get_pixel_extents(layout, &rect, NULL);
-	
+
 	bwidth = (sb*scale)+20;
 	bheight = rect.height+30;
 	if ( rect.width > bwidth ) bwidth = rect.width+20;
@@ -200,20 +200,20 @@ static void imagedisplay_add_scalebar(ImageDisplay *imagedisplay, GtkWidget *dra
 						xoffs+(x2), yoffs+imagedisplay->view_height-1-(y2)))
 
 static gboolean imagedisplay_redraw(GtkWidget *drawingarea, GdkEventExpose *event, ImageDisplay *imagedisplay) {
-	
+
 	double scale, xoffs, yoffs;
 	ImageDisplayMark *cur;
 	double max;
-	
+
 	xoffs = ((double)imagedisplay->drawingarea_width - imagedisplay->view_width) / 2;
 	yoffs = ((double)imagedisplay->drawingarea_height - imagedisplay->view_height) / 2;
 	scale = (double)imagedisplay->view_width/imagedisplay->imagerecord.width;
-	
+
 	gdk_draw_pixbuf(drawingarea->window, drawingarea->style->bg_gc[GTK_WIDGET_STATE(drawingarea)],
 			imagedisplay->pixbuf_scaled,
 			0, 0, xoffs, yoffs, imagedisplay->view_width, imagedisplay->view_height,
 			GDK_RGB_DITHER_NONE, 0, 0);
-	
+//return FALSE;
 	if ( imagedisplay->flags & IMAGEDISPLAY_SHOW_TILT_AXIS ) {
 		/* This is nasty, but works */
 		imagedisplay_draw_line(imagedisplay->gc_tiltaxis,
@@ -229,12 +229,12 @@ static gboolean imagedisplay_redraw(GtkWidget *drawingarea, GdkEventExpose *even
 			(imagedisplay->imagerecord.y_centre - imagedisplay->imagerecord.width
 					* tan(imagedisplay->imagerecord.omega)) * scale);
 	}
-	
+
 	/* Add scale bar */
 	if ( imagedisplay->flags & IMAGEDISPLAY_SCALE_BAR ) {
 		imagedisplay_add_scalebar(imagedisplay, drawingarea, scale, xoffs, yoffs);
 	}
-	
+
 	/* NB This calls the function above, which sorts out stuff */
 	if ( imagedisplay->flags & IMAGEDISPLAY_SHOW_CENTRE ) {
 		imagedisplay_draw_line(imagedisplay->gc_centre,
@@ -253,23 +253,23 @@ static gboolean imagedisplay_redraw(GtkWidget *drawingarea, GdkEventExpose *even
 	cur = imagedisplay->marks;
 	max = 0.0;
 	while ( cur ) {
-	
+
 		if ( cur->weight < 0.0 ) {
 			//printf("ID: Warning: ImageDisplayMark with negative weight\n");
 			cur = cur->next;
 			continue;
 		}
-		
+
 		if ( log(1+0.1*cur->weight) > max ) max = log(1+0.1*cur->weight);
 		cur = cur->next;
-		
+
 	}
-	
+
 	cur = imagedisplay->marks;
 	while ( cur ) {
-	
+
 		GdkGC *gc;
-		
+
 		switch ( cur->type ) {
 			case IMAGEDISPLAY_MARK_CIRCLE_1 : gc = imagedisplay->gc_marks_1; break;
 			case IMAGEDISPLAY_MARK_CIRCLE_2 : gc = imagedisplay->gc_marks_2; break;
@@ -278,39 +278,39 @@ static gboolean imagedisplay_redraw(GtkWidget *drawingarea, GdkEventExpose *even
 			case IMAGEDISPLAY_MARK_LINE_2 : gc = imagedisplay->gc_marks_2; break;
 			default : gc = imagedisplay->gc_marks_1; break;
 		}
-		
+
 		if ( (cur->type == IMAGEDISPLAY_MARK_CIRCLE_1)
 		  || (cur->type == IMAGEDISPLAY_MARK_CIRCLE_2)
 		  || (cur->type == IMAGEDISPLAY_MARK_CIRCLE_3) ) {
-		
+
 			double r;
-		
+
 			if ( cur->weight < 0.0 ) {
 				cur = cur->next;
 				continue;
 			}
-			
+
 			if ( cur->type == IMAGEDISPLAY_MARK_CIRCLE_1 ) {
 				r = 10.0;
 			} else {
-				r = 20.0 * (log(1+0.1*cur->weight)/max);
+				r = 10.0 * (log(1+0.1*cur->weight)/max);
 			}
-			
+
 			gdk_draw_arc(drawingarea->window, gc, FALSE,
 						xoffs + cur->x*scale - r,
 						yoffs + imagedisplay->view_height-1-cur->y*scale - r,
 						2*r, 2*r, 0, 64*360);
-			
+
 		} else if ( (cur->type == IMAGEDISPLAY_MARK_LINE_1) || (cur->type == IMAGEDISPLAY_MARK_LINE_2) ) {
-		
+
 			gdk_draw_line(drawingarea->window, gc,
 						xoffs + cur->x*scale,
 						yoffs + imagedisplay->view_height-1-cur->y*scale,
 						xoffs + cur->x2*scale,
 						yoffs + imagedisplay->view_height-1-cur->y2*scale );
-		
+
 		}
-		
+
 		cur = cur->next;
 	}
 
@@ -321,11 +321,11 @@ static gboolean imagedisplay_redraw(GtkWidget *drawingarea, GdkEventExpose *even
 static gint imagedisplay_realize(GtkWidget *widget, ImageDisplay *imagedisplay) {
 
 	GdkColor colour;
-	
+
 	imagedisplay->gc_centre = gdk_gc_new(imagedisplay->drawingarea->window);
 	gdk_color_parse("yellow", &colour);
 	gdk_gc_set_rgb_fg_color(imagedisplay->gc_centre, &colour);
-	
+
 	imagedisplay->gc_tiltaxis = gdk_gc_new(imagedisplay->drawingarea->window);
 	gdk_color_parse("#6600dd", &colour);
 	gdk_gc_set_rgb_fg_color(imagedisplay->gc_tiltaxis, &colour);
@@ -350,14 +350,14 @@ static gint imagedisplay_realize(GtkWidget *widget, ImageDisplay *imagedisplay)
 	imagedisplay->realised = TRUE;
 
 	return 0;
-	
+
 }
 
 ImageDisplay *imagedisplay_new_nowindow(ImageRecord imagerecord, ImageDisplayFlags flags, const char *message,
 					GCallback mouse_click_func, gpointer callback_data) {
 
 	ImageDisplay *imagedisplay;
-	
+
 	imagedisplay = malloc(sizeof(ImageDisplay));
 	imagedisplay->imagerecord = imagerecord;
 	imagedisplay->view_width = 512;
@@ -370,27 +370,27 @@ ImageDisplay *imagedisplay_new_nowindow(ImageRecord imagerecord, ImageDisplayFla
 	imagedisplay->pixbuf_scaled = NULL;
 	imagedisplay->realised = FALSE;
 	imagedisplay->window = NULL;
-	
+
 	imagedisplay_put_data(imagedisplay, imagerecord);
-	
+
 	imagedisplay->vbox = gtk_vbox_new(FALSE, 0);
-	
+
 	if ( message ) {
 		GtkWidget *label;
 		label = gtk_label_new(message);
 		gtk_misc_set_alignment(GTK_MISC(label), 0, 0.5);
 		gtk_box_pack_start(GTK_BOX(imagedisplay->vbox), label, FALSE, TRUE, 3);
 	}
-	
+
 	imagedisplay->drawingarea = gtk_drawing_area_new();
 	gtk_box_pack_start(GTK_BOX(imagedisplay->vbox), imagedisplay->drawingarea, TRUE, TRUE, 0);
-	
+
 	if ( imagedisplay->mouse_click_func ) {
 		gtk_widget_add_events(GTK_WIDGET(imagedisplay->drawingarea), GDK_BUTTON_PRESS_MASK);
 		g_signal_connect(GTK_OBJECT(imagedisplay->drawingarea), "button-press-event",
 					G_CALLBACK(imagedisplay->mouse_click_func), callback_data);
 	}
-	
+
 	g_signal_connect(GTK_OBJECT(imagedisplay->drawingarea), "realize",
 				G_CALLBACK(imagedisplay_realize), imagedisplay);
 	g_signal_connect(GTK_OBJECT(imagedisplay->drawingarea), "destroy",
@@ -399,7 +399,7 @@ ImageDisplay *imagedisplay_new_nowindow(ImageRecord imagerecord, ImageDisplayFla
 				G_CALLBACK(imagedisplay_configure_event), imagedisplay);
 	g_signal_connect(GTK_OBJECT(imagedisplay->drawingarea), "expose-event",
 				G_CALLBACK(imagedisplay_redraw), imagedisplay);
-	
+
 	return imagedisplay;
 
 }
@@ -411,24 +411,24 @@ ImageDisplay *imagedisplay_open_with_message(ImageRecord imagerecord, const char
 
 	ImageDisplay *imagedisplay;
 	GdkGeometry geom;
-	
+
 	imagedisplay = imagedisplay_new_nowindow(imagerecord, flags, message, mouse_click_func, callback_data);
-	
+
 	imagedisplay->window = gtk_window_new(GTK_WINDOW_TOPLEVEL);
 	gtk_container_add(GTK_CONTAINER(imagedisplay->window), imagedisplay->vbox);
 	imagedisplay->title = strdup(title);
 	gtk_window_set_title(GTK_WINDOW(imagedisplay->window), imagedisplay->title);
-	
+
 	geom.min_width = 128;	geom.min_height = 128;
 	gtk_window_set_geometry_hints(GTK_WINDOW(imagedisplay->window), GTK_WIDGET(imagedisplay->drawingarea),
 					&geom, GDK_HINT_MIN_SIZE);
-	
+
 	gtk_window_set_default_size(GTK_WINDOW(imagedisplay->window), 512, 512);
-	
+
 	gtk_widget_show_all(imagedisplay->window);
-	
+
 	return imagedisplay;
-	
+
 }
 
 ImageDisplay *imagedisplay_open(ImageRecord image, const char *title, ImageDisplayFlags flags) {
@@ -438,13 +438,13 @@ ImageDisplay *imagedisplay_open(ImageRecord image, const char *title, ImageDispl
 void imagedisplay_add_mark(ImageDisplay *imagedisplay, double x, double y, ImageDisplayMarkType type, double weight) {
 
 	ImageDisplayMark *new;
-	
+
 	new = malloc(sizeof(ImageDisplayMark));
 	new->x = x;  new->y = y;
 	new->type = type;
 	new->weight = weight;
 	new->next = NULL;
-	
+
 	if ( !imagedisplay->marks ) {
 		imagedisplay->marks = new;
 	} else {
@@ -461,14 +461,14 @@ void imagedisplay_add_line(ImageDisplay *imagedisplay, double x1, double y1,
 				double x2, double y2, ImageDisplayMarkType type) {
 
 	ImageDisplayMark *new;
-	
+
 	new = malloc(sizeof(ImageDisplayMark));
 	new->x = x1;  new->y = y1;
 	new->x2 = x2;  new->y2 = y2;
 	new->type = type;
 	new->weight = 1.0;	/* This field makes little sense for a line */
 	new->next = NULL;
-	
+
 	if ( !imagedisplay->marks ) {
 		imagedisplay->marks = new;
 	} else {
@@ -486,4 +486,3 @@ void imagedisplay_force_redraw(ImageDisplay *imagedisplay) {
 	gtk_widget_queue_draw_area(imagedisplay->drawingarea, 0, 0, imagedisplay->drawingarea_width,
 								    imagedisplay->drawingarea_height);
 }
-
diff --git a/src/intensities.c b/src/intensities.c
index e5ff1e2..88de5f2 100644
--- a/src/intensities.c
+++ b/src/intensities.c
@@ -184,6 +184,7 @@ static int intensities_do_save(ReflectionList *integrated, Basis *cell,
 	fh = fopen(filename, "w");
 
 	rcell = basis_get_cell(cell);
+	fprintf(fh, "#INT This file contains reflection intensities\n");
 	fprintf(fh, "a %12.8f\n", rcell.a*1e9);
 	fprintf(fh, "b %12.8f\n", rcell.b*1e9);
 	fprintf(fh, "c %12.8f\n", rcell.c*1e9);
diff --git a/src/itrans.c b/src/itrans.c
index 1018a66..fb5f649 100644
--- a/src/itrans.c
+++ b/src/itrans.c
@@ -3,8 +3,8 @@
  *
  * Parameterise features in an image for reconstruction
  *
- * (c) 2007 Thomas White <taw27@cam.ac.uk>
- *	    Gordon Ball <gfb21@cam.ac.uk>
+ * (c) 2007-2009 Thomas White <taw27@cam.ac.uk>
+ * (c) 2007      Gordon Ball <gfb21@cam.ac.uk>
  *
  *  dtr - Diffraction Tomography Reconstruction
  *
@@ -21,83 +21,107 @@
 #include "itrans-zaefferer.h"
 #include "itrans-stat.h"
 
+/* Return an empty feature list */
+ImageFeatureList *itrans_peaksearch_none(ImageRecord *imagerecord)
+{
+	return image_feature_list_new();
+}
+
 /* Return a feature list for the given image */
-ImageFeatureList *itrans_process_image(ImageRecord *image, PeakSearchMode psmode) {
-	
+ImageFeatureList *itrans_process_image(ImageRecord *image,
+						PeakSearchMode psmode)
+{
 	ImageFeatureList *list;
 
 	switch ( psmode ) {
-		case PEAKSEARCH_THRESHOLD :		{ list = itrans_peaksearch_threshold(image); break; }
-		case PEAKSEARCH_ADAPTIVE_THRESHOLD :	{ list = itrans_peaksearch_adaptive_threshold(image); break; }
-		case PEAKSEARCH_ZAEFFERER :		{ list = itrans_peaksearch_zaefferer(image); break; }
-		case PEAKSEARCH_STAT :			{ list = itrans_peaksearch_stat(image); break; }
+		case PEAKSEARCH_NONE : {
+			list = itrans_peaksearch_none(image);
+			break;
+		}
+		case PEAKSEARCH_THRESHOLD : {
+			list = itrans_peaksearch_threshold(image);
+			break;
+		}
+		case PEAKSEARCH_ADAPTIVE_THRESHOLD : {
+			list = itrans_peaksearch_adaptive_threshold(image);
+			break;
+		}
+		case PEAKSEARCH_ZAEFFERER : {
+			list = itrans_peaksearch_zaefferer(image);
+			break;
+		}
+		case PEAKSEARCH_STAT : {
+			list = itrans_peaksearch_stat(image);
+			break;
+		}
 		default: list = NULL;
 	}
-	
+
 	return list;
-	
 }
 
 /* Quantification radius */
 #define R 20
 
 /* Quantify each feature in the image's feature list */
-void itrans_quantify_features(ImageRecord *image) {
-
+void itrans_quantify_features(ImageRecord *image)
+{
 	int i;
-	
+
 	for ( i=0; i<image->features->n_features; i++ ) {
 
 		double theta;
 		double background, total;
 		int n;
 		int xi, yi;
-		
+
 		background = 0.0;
 		n = 0;
 		for ( theta=0; theta<2*M_PI; theta+=2*M_PI/10 ) {
-		
+
 			double x, y;
 			long int xd, yd;
-			
+
 			x = image->features->features[i].x + R*sin(theta);
 			y = image->features->features[i].y + R*cos(theta);
 			xd = rint(x);
 			yd = rint(y);
-			
-			if ( (xd>=0) && (yd>=0) && (xd<image->width) && (yd<image->height) ) {
-				background += image->image[xd + yd*image->width];
+
+			if ( (xd>=0) && (yd>=0)
+			  && (xd<image->width) && (yd<image->height) ) {
+				background += image->image[xd+yd*image->width];
 				n++;
 			}
 
 		}
 		background /= n;
-		
+
 		total = 0.0;
 		for ( xi=-20; xi<=20; xi++ ) {
 			for ( yi=-20; yi<=20; yi++ ) {
-			
+
 				double x, y;
 				long int xd, yd;
-				
+
 				if ( xi*xi + yi*yi > 20 ) continue;
-				
+
 				x = image->features->features[i].x + xi;
 				y = image->features->features[i].y + yi;
-				
+
 				xd = rint(x);
 				yd = rint(y);
-				
-				if ( (xd>=0) && (yd>=0) && (xd<image->width) && (yd<image->height) ) {
-					total += image->image[xd + yd*image->width] - background;
+
+				if ( (xd>=0) && (yd>=0)
+				  && (xd<image->width) && (yd<image->height) ) {
+					total +=
+						image->image[xd+yd*image->width]
+								- background;
 				}
-			
+
 			}
 		}
-		
+
 		image->features->features[i].intensity = total;
 
 	}
-
 }
-
diff --git a/src/main.c b/src/main.c
index 494dc46..a70afaf 100644
--- a/src/main.c
+++ b/src/main.c
@@ -3,7 +3,7 @@
  *
  * The Top Level Source File
  *
- * (c) 2007-2008 Thomas White <taw27@cam.ac.uk>
+ * (c) 2007-2009 Thomas White <taw27@cam.ac.uk>
  * (c) 2007      Gordon Ball <gfb21@cam.ac.uk>
  *
  *  dtr - Diffraction Tomography Reconstruction
@@ -36,108 +36,123 @@
 #include "dirax.h"
 #include "itrans.h"
 
-void main_do_reconstruction(ControlContext *ctx) {
-
+void main_do_reconstruction(ControlContext *ctx)
+{
 	if ( ctx->inputfiletype != INPUT_DRX ) {
-		
+
 		int val = 0;
-		
+
 		/* Initial centering */
-		prealign_sum_stack(ctx->images, ctx->have_centres, ctx->sum_stack);
+		prealign_sum_stack(ctx->images, ctx->have_centres,
+								ctx->sum_stack);
 		if ( ctx->finecentering ) {
 			prealign_fine_centering(ctx->images, ctx->sum_stack);
 		}
-		
+
 		if ( ctx->cache_filename == NULL ) {
-			
+
 			int i;
-		
+
 			/* Find all the features */
 			printf("MA: Analysing images..."); fflush(stdout);
 			for ( i=0; i<ctx->images->n_images; i++ ) {
-				ctx->images->images[i].features = itrans_process_image(&ctx->images->images[i],
-										       ctx->psmode);
-				itrans_quantify_features(&ctx->images->images[i]);
+				ctx->images->images[i].features =
+				   itrans_process_image(&ctx->images->images[i],
+				   ctx->psmode);
+				itrans_quantify_features(
+						       &ctx->images->images[i]);
 			}
 			printf("done.\n");
-		
+
 		} else {
-		
+
 			int i;
-			
-			printf("MA: Loading previous image analysis from '%s'\n", ctx->cache_filename);
-			if ( cache_load(ctx->images, ctx->cache_filename) ) val = 1;
-			
-			/* Quantify all the features (unnecessary if the file already contains intensities,
-			 *	but many of mine don't, and it doesn't hurt. */
+
+			printf("MA: Loading previous image analysis from %s\n",
+							ctx->cache_filename);
+			if ( cache_load(ctx->images, ctx->cache_filename) )
+									val = 1;
+
+			/* Quantify all the features (unnecessary if the file
+			 *  already contains intensities,
+			 *  but many of mine don't, and it doesn't hurt. */
 			printf("MA: Quantifying features..."); fflush(stdout);
 			for ( i=0; i<ctx->images->n_images; i++ ) {
-				itrans_quantify_features(&ctx->images->images[i]);
+				itrans_quantify_features(
+						       &ctx->images->images[i]);
 			}
 			printf("done.\n");
-		
+
 		}
-		
+
 		if ( ctx->finecentering ) {
 			prealign_feature_centering(ctx->images);
 		}
-		
+
 		if ( !val ) mapping_map_features(ctx);
-	
+
 	} /* else has already been created by dirax_load() */
-	
+
 	if ( ctx->reflectionlist ) {
 		ctx->dw = displaywindow_open(ctx);
 	} else {
 		fprintf(stderr, "Reconstruction failed.\n");
 		gtk_exit(0);
 	}
-	
 }
 
-static gint main_method_window_response(GtkWidget *method_window, gint response, ControlContext *ctx) {
-
+static gint main_method_window_response(GtkWidget *method_window, gint response,
+					ControlContext *ctx)
+{
 	if ( response == GTK_RESPONSE_OK ) {
-	
+
 		int val = 0;
-		
-		switch ( gtk_combo_box_get_active(GTK_COMBO_BOX(ctx->combo_peaksearch)) ) {
-			case 0 : ctx->psmode = PEAKSEARCH_THRESHOLD; break;
-			case 1 : ctx->psmode = PEAKSEARCH_ADAPTIVE_THRESHOLD; break;
-			case 2 : ctx->psmode = PEAKSEARCH_ZAEFFERER; break;
-			case 3 : ctx->psmode = PEAKSEARCH_STAT; break;
-			case 4 : ctx->psmode = PEAKSEARCH_CACHED; break;
-			default: ctx->psmode = PEAKSEARCH_NONE; break;	/* Happens when reading from a cache file */
+
+		switch ( gtk_combo_box_get_active(
+				GTK_COMBO_BOX(ctx->combo_peaksearch)) ) {
+			case 0 : ctx->psmode = PEAKSEARCH_NONE; break;
+			case 1 : ctx->psmode = PEAKSEARCH_THRESHOLD; break;
+			case 2 : ctx->psmode = PEAKSEARCH_ADAPTIVE_THRESHOLD;
+									break;
+			case 3 : ctx->psmode = PEAKSEARCH_ZAEFFERER; break;
+			case 4 : ctx->psmode = PEAKSEARCH_STAT; break;
+			case 5 : ctx->psmode = PEAKSEARCH_CACHED; break;
+			/* Happens when reading from a cache file */
+			default: ctx->psmode = PEAKSEARCH_NONE; break;
 		}
-		
-		if ( gtk_toggle_button_get_active(GTK_TOGGLE_BUTTON(ctx->checkbox_prealign)) ) {
+
+		if ( gtk_toggle_button_get_active(
+				GTK_TOGGLE_BUTTON(ctx->checkbox_prealign)) ) {
 			ctx->prealign = TRUE;
 		} else {
 			ctx->prealign = FALSE;
 		}
-		
-		if ( gtk_toggle_button_get_active(GTK_TOGGLE_BUTTON(ctx->checkbox_finecentering)) ) {
+
+		if ( gtk_toggle_button_get_active(
+			GTK_TOGGLE_BUTTON(ctx->checkbox_finecentering)) ) {
 			ctx->finecentering = TRUE;
 		} else {
 			ctx->finecentering = FALSE;
 		}
-		
-		if ( gtk_toggle_button_get_active(GTK_TOGGLE_BUTTON(ctx->checkbox_sumstack)) ) {
+
+		if ( gtk_toggle_button_get_active(
+				GTK_TOGGLE_BUTTON(ctx->checkbox_sumstack)) ) {
 			ctx->sum_stack = TRUE;
 		} else {
 			ctx->sum_stack = FALSE;
 		}
-		
+
 		if ( ctx->psmode == PEAKSEARCH_CACHED ) {
-			ctx->cache_filename = gtk_file_chooser_get_filename(GTK_FILE_CHOOSER(ctx->cache_file_selector));
+			ctx->cache_filename = gtk_file_chooser_get_filename(
+				GTK_FILE_CHOOSER(ctx->cache_file_selector));
 			if ( !ctx->cache_filename ) {
 				return 1;
 			}
 		}
-		
+
 		gtk_widget_destroy(method_window);
 		while ( gtk_events_pending() ) gtk_main_iteration();
-		
+
 		/* Load the input */
 		if ( ctx->inputfiletype == INPUT_QDRP ) {
 			val = qdrp_read(ctx);
@@ -147,139 +162,170 @@ static gint main_method_window_response(GtkWidget *method_window, gint response,
 			ctx->reflectionlist = dirax_load(ctx->filename);
 			if ( !ctx->reflectionlist ) val = 1;
 		}
-		
+
 		if ( val ) {
 			fprintf(stderr, "Reconstruction failed.\n");
 			gtk_exit(0);
 		}
-		
+
 		if ( ctx->prealign ) {
-			prealign_do_series(ctx);	/* this will eventually call main_do_reconstruction() */
+			/* this will eventually call main_do_reconstruction() */
+			prealign_do_series(ctx);
 		} else {
 			main_do_reconstruction(ctx);
 		}
-		
+
 	} else {
 		gtk_exit(0);
 	}
-	
-	return 0;
 
+	return 0;
 }
 
-static gint main_peaksearch_changed(GtkWidget *method_window, ControlContext *ctx) {
-
-	if ( gtk_combo_box_get_active(GTK_COMBO_BOX(ctx->combo_peaksearch)) == 4 ) {
-		gtk_widget_set_sensitive(GTK_WIDGET(ctx->cache_file_selector), TRUE);
+static gint main_peaksearch_changed(GtkWidget *method_window,
+					ControlContext *ctx)
+{
+	if ( gtk_combo_box_get_active(
+				GTK_COMBO_BOX(ctx->combo_peaksearch)) == 5 ) {
+		gtk_widget_set_sensitive(GTK_WIDGET(ctx->cache_file_selector),
+						TRUE);
 	} else {
-		gtk_widget_set_sensitive(GTK_WIDGET(ctx->cache_file_selector), FALSE);
+		gtk_widget_set_sensitive(GTK_WIDGET(ctx->cache_file_selector),
+						FALSE);
 	}
 
 	return 0;
-
 }
 
-void main_method_dialog_open(ControlContext *ctx) {
-
+void main_method_dialog_open(ControlContext *ctx)
+{
 	GtkWidget *method_window;
 	GtkWidget *vbox;
 	GtkWidget *hbox;
 	GtkWidget *table;
 	GtkWidget *peaksearch_label;
 	GtkWidget *cache_file_selector_label;
-	
-	method_window = gtk_dialog_new_with_buttons("Reconstruction Parameters", NULL,
-		GTK_DIALOG_DESTROY_WITH_PARENT,	GTK_STOCK_CANCEL, GTK_RESPONSE_CLOSE,
-						GTK_STOCK_OK, GTK_RESPONSE_OK, NULL);
+
+	method_window = gtk_dialog_new_with_buttons("Reconstruction Parameters",
+					NULL,
+					GTK_DIALOG_DESTROY_WITH_PARENT,
+					GTK_STOCK_CANCEL, GTK_RESPONSE_CLOSE,
+					GTK_STOCK_OK, GTK_RESPONSE_OK, NULL);
 	gtk_window_set_default_size(GTK_WINDOW(method_window), 400, -1);
-	
+
 	vbox = gtk_vbox_new(FALSE, 0);
 	hbox = gtk_hbox_new(TRUE, 0);
-	gtk_box_pack_start(GTK_BOX(GTK_DIALOG(method_window)->vbox), GTK_WIDGET(hbox), FALSE, FALSE, 7);
+	gtk_box_pack_start(GTK_BOX(GTK_DIALOG(method_window)->vbox),
+					GTK_WIDGET(hbox), FALSE, FALSE, 7);
 	gtk_box_pack_start(GTK_BOX(hbox), GTK_WIDGET(vbox), FALSE, TRUE, 10);
-	
+
 	table = gtk_table_new(5, 2, FALSE);
 	gtk_table_set_row_spacings(GTK_TABLE(table), 5);
-	
+
 	peaksearch_label = gtk_label_new("Peak Search: ");
-	gtk_table_attach_defaults(GTK_TABLE(table), peaksearch_label, 1, 2, 1, 2);
+	gtk_table_attach_defaults(GTK_TABLE(table), peaksearch_label,
+								1, 2, 1, 2);
 	gtk_misc_set_alignment(GTK_MISC(peaksearch_label), 1, 0.5);
 	ctx->combo_peaksearch = gtk_combo_box_new_text();
-	gtk_combo_box_append_text(GTK_COMBO_BOX(ctx->combo_peaksearch), "Simple Thresholding");
-	gtk_combo_box_append_text(GTK_COMBO_BOX(ctx->combo_peaksearch), "Adaptive Thresholding");
-	gtk_combo_box_append_text(GTK_COMBO_BOX(ctx->combo_peaksearch), "Zaefferer Gradient Search");
-	gtk_combo_box_append_text(GTK_COMBO_BOX(ctx->combo_peaksearch), "Iterative Statistical Analysis");
-	gtk_combo_box_append_text(GTK_COMBO_BOX(ctx->combo_peaksearch), "Get From Cache File");
-	gtk_combo_box_set_active(GTK_COMBO_BOX(ctx->combo_peaksearch), 2);
-	gtk_table_attach_defaults(GTK_TABLE(table), ctx->combo_peaksearch, 2, 3, 1, 2);
-	g_signal_connect(G_OBJECT(ctx->combo_peaksearch), "changed", G_CALLBACK(main_peaksearch_changed), ctx);
-	
+	gtk_combo_box_append_text(GTK_COMBO_BOX(ctx->combo_peaksearch),
+					"None");
+	gtk_combo_box_append_text(GTK_COMBO_BOX(ctx->combo_peaksearch),
+					"Simple Thresholding");
+	gtk_combo_box_append_text(GTK_COMBO_BOX(ctx->combo_peaksearch),
+					"Adaptive Thresholding");
+	gtk_combo_box_append_text(GTK_COMBO_BOX(ctx->combo_peaksearch),
+					"Zaefferer Gradient Search");
+	gtk_combo_box_append_text(GTK_COMBO_BOX(ctx->combo_peaksearch),
+					"Iterative Statistical Analysis");
+	gtk_combo_box_append_text(GTK_COMBO_BOX(ctx->combo_peaksearch),
+					"Get From Cache File");
+	gtk_combo_box_set_active(GTK_COMBO_BOX(ctx->combo_peaksearch), 3);
+	gtk_table_attach_defaults(GTK_TABLE(table), ctx->combo_peaksearch,
+								2, 3, 1, 2);
+	g_signal_connect(G_OBJECT(ctx->combo_peaksearch), "changed",
+				G_CALLBACK(main_peaksearch_changed), ctx);
+
 	cache_file_selector_label = gtk_label_new("Cache File to Load: ");
-	gtk_table_attach_defaults(GTK_TABLE(table), cache_file_selector_label, 1, 2, 2, 3);
+	gtk_table_attach_defaults(GTK_TABLE(table), cache_file_selector_label,
+								1, 2, 2, 3);
 	gtk_misc_set_alignment(GTK_MISC(cache_file_selector_label), 1, 0.5);
-	ctx->cache_file_selector = gtk_file_chooser_button_new("Select Cache File to Load",
-							       GTK_FILE_CHOOSER_ACTION_OPEN);
-	gtk_table_attach_defaults(GTK_TABLE(table), ctx->cache_file_selector, 2, 3, 2, 3);
+	ctx->cache_file_selector = gtk_file_chooser_button_new(
+						"Select Cache File to Load",
+						GTK_FILE_CHOOSER_ACTION_OPEN);
+	gtk_table_attach_defaults(GTK_TABLE(table), ctx->cache_file_selector,
+						2, 3, 2, 3);
 	gtk_widget_set_sensitive(GTK_WIDGET(ctx->cache_file_selector), FALSE);
-	
-	ctx->checkbox_prealign = gtk_check_button_new_with_label("Manually pre-align the image stack");
-	gtk_table_attach_defaults(GTK_TABLE(table), ctx->checkbox_prealign, 1, 3, 3, 4);
-	
-	ctx->checkbox_finecentering = gtk_check_button_new_with_label("Perform fine pattern centering");
-	gtk_table_attach_defaults(GTK_TABLE(table), ctx->checkbox_finecentering, 1, 3, 4, 5);
-	gtk_toggle_button_set_active(GTK_TOGGLE_BUTTON(ctx->checkbox_finecentering), TRUE);
-	
-	ctx->checkbox_sumstack = gtk_check_button_new_with_label("Show summed image stack");
-	gtk_table_attach_defaults(GTK_TABLE(table), ctx->checkbox_sumstack, 1, 3, 5, 6);
-	
+
+	ctx->checkbox_prealign = gtk_check_button_new_with_label(
+					"Manually pre-align the image stack");
+	gtk_table_attach_defaults(GTK_TABLE(table), ctx->checkbox_prealign,
+					1, 3, 3, 4);
+
+	ctx->checkbox_finecentering = gtk_check_button_new_with_label(
+					"Perform fine pattern centering");
+	gtk_table_attach_defaults(GTK_TABLE(table), ctx->checkbox_finecentering,
+					1, 3, 4, 5);
+	gtk_toggle_button_set_active(
+			GTK_TOGGLE_BUTTON(ctx->checkbox_finecentering), TRUE);
+
+	ctx->checkbox_sumstack = gtk_check_button_new_with_label(
+					"Show summed image stack");
+	gtk_table_attach_defaults(GTK_TABLE(table), ctx->checkbox_sumstack,
+					1, 3, 5, 6);
+
 	if ( ctx->inputfiletype == INPUT_DRX ) {
-		gtk_combo_box_append_text(GTK_COMBO_BOX(ctx->combo_peaksearch), "3D coordinates from DirAx file");
-		gtk_widget_set_sensitive(GTK_WIDGET(ctx->combo_peaksearch), FALSE);
-		gtk_combo_box_set_active(GTK_COMBO_BOX(ctx->combo_peaksearch), 5);
+		gtk_combo_box_append_text(GTK_COMBO_BOX(ctx->combo_peaksearch),
+					"3D coordinates from DirAx file");
+		gtk_widget_set_sensitive(GTK_WIDGET(ctx->combo_peaksearch),
+						FALSE);
+		gtk_combo_box_set_active(GTK_COMBO_BOX(ctx->combo_peaksearch),
+									6);
 	}
-	
+
 	gtk_box_pack_start(GTK_BOX(vbox), GTK_WIDGET(table), TRUE, TRUE, 5);
-	
-	g_signal_connect(G_OBJECT(method_window), "response", G_CALLBACK(main_method_window_response), ctx);
-	
-	gtk_widget_show_all(method_window);
 
-}
+	g_signal_connect(G_OBJECT(method_window), "response",
+				G_CALLBACK(main_method_window_response), ctx);
 
+	gtk_widget_show_all(method_window);
+}
 
-int main(int argc, char *argv[]) {
 
+int main(int argc, char *argv[])
+{
 	ControlContext *ctx;
 	struct stat stat_buffer;
 	FILE *fh;
 
 	gtk_init(&argc, &argv);
-	
+
 	if ( gtk_gl_init_check(&argc, &argv) == FALSE ) {
 		fprintf(stderr, "Could not initialise gtkglext\n");
 		return 1;
 	}
-	
+
 	if ( gdk_gl_query_extension() == FALSE ) {
 		fprintf(stderr, "OpenGL not supported\n");
 		return 1;
 	}
-	
+
 	if ( argc != 2 ) {
-		fprintf(stderr, "Syntax: %s [<MRC file> | qdrp.rc | reflect.cache]\n", argv[0]);
+		fprintf(stderr,
+			"Syntax: %s [<MRC file> | qdrp.rc | reflect.cache]\n",
+			argv[0]);
 		return 1;
 	}
-	
+
 	ctx = control_ctx_new();
 	ctx->filename = strdup(argv[1]);
-	
+
 	fh = fopen(ctx->filename, "r");
 	if ( !fh ) {
 		printf("Couldn't open file '%s'\n", ctx->filename);
 		return 1;
 	}
 	fclose(fh);
-	
+
 	if ( qdrp_is_qdrprc(ctx->filename) ) {
 		printf("QDRP input file detected.\n");
 		ctx->inputfiletype = INPUT_QDRP;
@@ -293,20 +339,18 @@ int main(int argc, char *argv[]) {
 		fprintf(stderr, "Unrecognised input file type\n");
 		return 1;
 	}
-	
+
 	if ( stat(argv[1], &stat_buffer) == -1 ) {
 		fprintf(stderr, "File '%s' not found\n", argv[1]);
 		return 1;
 	}
-	
+
 	main_method_dialog_open(ctx);
-		
+
 	gtk_main();
-	
+
 	free(ctx->filename);
 	free(ctx);
-	
+
 	return 0;
-	
 }
-
diff --git a/src/mapping.c b/src/mapping.c
index 16dbe34..4ffc57d 100644
--- a/src/mapping.c
+++ b/src/mapping.c
@@ -25,29 +25,29 @@ void mapping_rotate(double x, double y, double z, double *ddx, double *ddy, doub
 
 	double nx, ny, nz;
 	double x_temp, y_temp, z_temp;
-	
+
 	/* First: rotate image clockwise until tilt axis is aligned horizontally. */
 	nx = x*cos(omega) + y*sin(omega);
 	ny = -x*sin(omega) + y*cos(omega);
 	nz = z;
-	
+
 	/* Now, tilt about the x-axis ANTICLOCKWISE around +x, i.e. the "wrong" way.
 		This is because the crystal is rotated in the experiment, not the Ewald sphere. */
-	x_temp = nx; y_temp = ny; z_temp = nz;	
+	x_temp = nx; y_temp = ny; z_temp = nz;
 	nx = x_temp;
 	ny = cos(tilt)*y_temp + sin(tilt)*z_temp;
 	nz = -sin(tilt)*y_temp + cos(tilt)*z_temp;
-	
+
 	/* Finally, reverse the omega rotation to restore the location of the image in 3D space */
-	x_temp = nx; y_temp = ny; z_temp = nz;	
+	x_temp = nx; y_temp = ny; z_temp = nz;
 	nx = x_temp*cos(-omega) + y_temp*sin(-omega);
 	ny = -x_temp*sin(-omega) + y_temp*cos(-omega);
 	nz = z_temp;
-	
+
 	*ddx = nx;
 	*ddy = ny;
 	*ddz = nz;
-	
+
 }
 
 int mapping_map_to_space(ImageFeature *refl, double *ddx, double *ddy, double *ddz, double *twotheta) {
@@ -55,22 +55,22 @@ int mapping_map_to_space(ImageFeature *refl, double *ddx, double *ddy, double *d
 	/* "Input" space */
 	double d, x, y;
 	ImageRecord *imagerecord;
-	
+
 	/* Angular description of reflection */
 	double theta, psi, k;
-	
+
 	/* Reciprocal space */
 	double tilt;
 	double omega;
-	
+
 	double x_temp, y_temp, z_temp;
-	
+
 	imagerecord = refl->parent;
 	x = refl->x - imagerecord->x_centre;	y = refl->y - imagerecord->y_centre;
 	tilt = imagerecord->tilt;
 	omega = imagerecord->omega;
 	k = 1/imagerecord->lambda;
-	
+
 	/* Calculate an angular description of the reflection */
 	if ( imagerecord->fmode == FORMULATION_CLEN ) {
 		x /= imagerecord->resolution;
@@ -87,15 +87,15 @@ int mapping_map_to_space(ImageFeature *refl, double *ddx, double *ddy, double *d
 		return -1;
 	}
 	psi = atan2(y, x);
-	
+
 	x_temp = k*sin(theta)*cos(psi);
 	y_temp = k*sin(theta)*sin(psi);
 	z_temp = k - k*cos(theta);
-	
+
 	mapping_rotate(x_temp, y_temp, z_temp, ddx, ddy, ddz, omega, tilt);
-	
+
 	*twotheta = theta;	/* Radians.  I've used the "wrong" nomenclature above */
-	
+
 	return 0;
 
 }
@@ -105,17 +105,17 @@ int mapping_scale(ImageFeature *refl, double *ddx, double *ddy) {
 	double x, y;
 	ImageRecord *imagerecord;
 	double k;
-	
+
 	imagerecord = refl->parent;
 	x = refl->x - imagerecord->x_centre;
 	y = refl->y - imagerecord->y_centre;
 	k = 1/imagerecord->lambda;
-	
+
 	if ( imagerecord->fmode == FORMULATION_CLEN ) {
 		x /= imagerecord->resolution;
 		y /= imagerecord->resolution;	/* Convert pixels to metres */
 		*ddx = x * k / imagerecord->camera_len;
-		*ddy = y * k / imagerecord->camera_len;	
+		*ddy = y * k / imagerecord->camera_len;
 	} else if (imagerecord->fmode == FORMULATION_PIXELSIZE ) {
 		*ddx = x * imagerecord->pixel_size;
 		*ddy = y * imagerecord->pixel_size;	/* Convert pixels to metres^-1 */
@@ -123,9 +123,9 @@ int mapping_scale(ImageFeature *refl, double *ddx, double *ddy) {
 		fprintf(stderr, "Unrecognised formulation mode in mapping_scale.\n");
 		return -1;
 	}
-	
+
 	return 0;
-	
+
 }
 
 /* Return the length of a 1 nm^-1 scale bar in the given image (in pixels)
@@ -137,7 +137,7 @@ double mapping_scale_bar_length(ImageRecord *image) {
 		case FORMULATION_CLEN : return 1.0e9*image->resolution*image->camera_len*image->lambda;
 		default : fprintf(stderr, "Unrecognised formulation mode in mapping_scale_bar_length.\n");
 	}
-	
+
 	return 0.0;
 
 }
@@ -145,29 +145,31 @@ double mapping_scale_bar_length(ImageRecord *image) {
 void mapping_map_features(ControlContext *ctx) {
 
 	int i;
-	
+
 	/* Create reflection list for measured reflections */
 	if ( ctx->reflectionlist ) reflectionlist_free(ctx->reflectionlist);
 	ctx->reflectionlist = reflectionlist_new();
-	
+
 	printf("MP: Mapping to 3D..."); fflush(stdout);
 	for ( i=0; i<ctx->images->n_images; i++ ) {
-	
+
 		int j;
-		
+
 		/* Iterate over the features in this image */
 		for ( j=0; j<ctx->images->images[i].features->n_features; j++ ) {
-		
+
 			double nx, ny, nz, twotheta;
-			
+
 			if ( !mapping_map_to_space(&ctx->images->images[i].features->features[j],
 											&nx, &ny, &nz, &twotheta) ) {
 				reflection_add(ctx->reflectionlist, nx, ny, nz,
 					ctx->images->images[i].features->features[j].intensity, REFLECTION_NORMAL);
+			} else {
+				printf("Couldn't map\n");
 			}
-			
+
 		}
-		
+
 	}
 	printf("done.\n");
 
@@ -182,12 +184,11 @@ void mapping_adjust_axis(ControlContext *ctx, double offset) {
 							rad2deg(ctx->images->images[i].omega+offset));
 		ctx->images->images[i].omega += offset;
 	}
-	
+
 	mapping_map_features(ctx);
 	if ( ctx->dw ) {
 		displaywindow_update_imagestack(ctx->dw);
 		displaywindow_update(ctx->dw);
 	}
-	
-}
 
+}
diff --git a/src/reflections.c b/src/reflections.c
index a589231..ce5f7a9 100644
--- a/src/reflections.c
+++ b/src/reflections.c
@@ -32,9 +32,9 @@ static void reflectionlist_init(ReflectionList *reflectionlist) {
 ReflectionList *reflectionlist_new() {
 
 	ReflectionList *reflectionlist = malloc(sizeof(ReflectionList));
-	
+
 	reflectionlist_init(reflectionlist);
-	
+
 	return reflectionlist;
 
 }
@@ -44,10 +44,10 @@ void reflectionlist_clear_markers(ReflectionList *reflectionlist) {
 	Reflection *reflection = reflectionlist->reflections;
 	Reflection *prev = NULL;
 	int del = 0;
-	
+
 	 while ( reflection ) {
 		Reflection *next = reflection->next;
-		
+
 		if ( (reflection->type == REFLECTION_MARKER) || (reflection->type == REFLECTION_GENERATED)
 		 || (reflection->type == REFLECTION_VECTOR_MARKER_1) || (reflection->type == REFLECTION_VECTOR_MARKER_2)
 		 || (reflection->type == REFLECTION_VECTOR_MARKER_3) ) {
@@ -61,14 +61,14 @@ void reflectionlist_clear_markers(ReflectionList *reflectionlist) {
 		} else {
 			prev = reflection;
 		}
-		
+
 		reflection = next;
-	
+
 	};
-	
+
 	reflectionlist->n_reflections -= del;
 	reflectionlist->last_reflection = prev;
-	
+
 }
 void reflectionlist_clear(ReflectionList *reflectionlist) {
 
@@ -78,7 +78,7 @@ void reflectionlist_clear(ReflectionList *reflectionlist) {
 		free(reflection);
 		reflection = next;
 	};
-	
+
 	reflectionlist_init(reflectionlist);
 
 }
@@ -92,27 +92,32 @@ Reflection *reflection_add(ReflectionList *reflectionlist, double x, double y, d
 
 	Reflection *new_reflection;
 	Reflection *nearest;
-	
+
 	if ( reflectionlist->list_capped ) return NULL;
-	
+
 	if ( reflectionlist->n_reflections > 1e7 ) {
 		fprintf(stderr, "More than 10 million reflections on list.  I think this is silly.\n");
 		fprintf(stderr, "No further reflections will be stored.  Go and fix the peak detection.\n");
 		reflectionlist->list_capped = 1;
 	}
-	
-	nearest = reflectionlist_find_nearest_type(reflectionlist, x, y, z, type);
-	if ( nearest && distance3d(x, y, z, nearest->x, nearest->y, nearest->z) < 0.1e9 ) return NULL;
-	
+
+//	nearest = reflectionlist_find_nearest_type(reflectionlist, x, y, z, type);
+//	if ( nearest && distance3d(x, y, z, nearest->x, nearest->y, nearest->z) < 0.1e9 ) {
+//	printf("Too close\n");
+//	return NULL;
+//}
 	new_reflection = malloc(sizeof(Reflection));
 	new_reflection->next = NULL;
 	new_reflection->x = x;
 	new_reflection->y = y;
 	new_reflection->z = z;
+	new_reflection->h = 999;
+	new_reflection->k = 333;
+	new_reflection->l = 111;
 	new_reflection->intensity = intensity;
 	new_reflection->type = type;
 	new_reflection->found = 0;
-	
+
 	if ( reflectionlist->last_reflection ) {
 		reflectionlist->last_reflection->next = new_reflection;
 		reflectionlist->last_reflection = new_reflection;
@@ -121,16 +126,16 @@ Reflection *reflection_add(ReflectionList *reflectionlist, double x, double y, d
 		reflectionlist->last_reflection = new_reflection;
 	}
 	reflectionlist->n_reflections++;
-	
+
 	return new_reflection;
-	
+
 }
 
 double reflectionlist_largest_g(ReflectionList *reflectionlist) {
 
 	double max = 0.0;
 	Reflection *reflection;
-	
+
 	reflection = reflectionlist->reflections;
 	while ( reflection ) {
 		if ( reflection->type == REFLECTION_NORMAL ) {
@@ -140,7 +145,7 @@ double reflectionlist_largest_g(ReflectionList *reflectionlist) {
 		}
 		reflection = reflection->next;
 	};
-	
+
 	return max;
 
 }
@@ -150,7 +155,7 @@ Reflection *reflectionlist_find_nearest(ReflectionList *reflectionlist, double x
 	double max = +INFINITY;
 	Reflection *reflection;
 	Reflection *best = NULL;
-	
+
 	reflection = reflectionlist->reflections;
 	while ( reflection ) {
 		if ( reflection->type == REFLECTION_NORMAL ) {
@@ -163,7 +168,7 @@ Reflection *reflectionlist_find_nearest(ReflectionList *reflectionlist, double x
 		}
 		reflection = reflection->next;
 	};
-	
+
 	return best;
 
 }
@@ -174,7 +179,7 @@ Reflection *reflectionlist_find_nearest_longer_unknown(ReflectionList *reflectio
 	double max = +INFINITY;
 	Reflection *reflection;
 	Reflection *best = NULL;
-	
+
 	reflection = reflectionlist->reflections;
 	while ( reflection ) {
 		if ( (reflection->type == REFLECTION_NORMAL) && (!reflection->found) ) {
@@ -187,7 +192,7 @@ Reflection *reflectionlist_find_nearest_longer_unknown(ReflectionList *reflectio
 		}
 		reflection = reflection->next;
 	};
-	
+
 	return best;
 
 }
@@ -198,7 +203,7 @@ Reflection *reflectionlist_find_nearest_type(ReflectionList *reflectionlist, dou
 	double max = +INFINITY;
 	Reflection *reflection;
 	Reflection *best = NULL;
-	
+
 	reflection = reflectionlist->reflections;
 	while ( reflection ) {
 		if ( reflection->type == type ) {
@@ -211,7 +216,7 @@ Reflection *reflectionlist_find_nearest_type(ReflectionList *reflectionlist, dou
 		}
 		reflection = reflection->next;
 	};
-	
+
 	return best;
 
 }
@@ -223,9 +228,9 @@ ReflectionList *reflection_list_from_cell(Basis *basis) {
 	double max_res;
 	signed int h, k, l;
 	int max_order_a, max_order_b, max_order_c;
-	
+
 	ordered = reflectionlist_new();
-	
+
 	max_res = 21e9;
 	do {
 		max_order_a = max_res/modulus(basis->a.x, basis->a.y, basis->a.z);
@@ -234,36 +239,35 @@ ReflectionList *reflection_list_from_cell(Basis *basis) {
 		max_res -= 1e9;
 	} while ( (max_order_a * max_order_b * max_order_c * 8) > 1e4 );
 	printf("Selected maximum resolution %8.5f nm^-1\n", max_res/1e9);
-	
+
 	for ( h=-max_order_a; h<=max_order_a; h++ ) {
 		for ( k=-max_order_b; k<=max_order_b; k++ ) {
 			for ( l=-max_order_c; l<=max_order_c; l++ ) {
 
 				double x, y, z;
-				
-				/* Test mode */
-				//if ( h != 0 ) continue;
-				
+
 				x = h*basis->a.x + k*basis->b.x + l*basis->c.x;
 				y = h*basis->a.y + k*basis->b.y + l*basis->c.y;
 				z = h*basis->a.z + k*basis->b.z + l*basis->c.z;
-				
+
 				if ( ( x*x + y*y + z*z ) <= max_res*max_res ) {
-					Reflection *ref;
-					ref = reflection_add(ordered, x, y, z, 1.0, REFLECTION_GENERATED);
-					if ( ref ) {	/* Check it's not NULL */
-						ref->h = h;	ref->k = k;	ref->l = l;
-					}
 					if ( (h!=0) || (k!=0) || (l!=0) ) {
-					//	reflection_add(ctx->reflectionlist, x, y, z, 1.0, REFLECTION_GENERATED);
-						reflection_add(ordered, x, y, z, 1.0, REFLECTION_GENERATED);
+						Reflection *ref;
+						ref = reflection_add(ordered,
+							x, y, z, 1.0,
+							REFLECTION_GENERATED);
+						if ( ref ) {
+							ref->h = h;
+							ref->k = k;
+							ref->l = l;
+						}
 					}
 				}
-				
+
 			}
 		}
 	}
-	
+
 	return ordered;
 
 }
@@ -271,19 +275,19 @@ ReflectionList *reflection_list_from_cell(Basis *basis) {
 void reflection_list_from_new_cell(ReflectionList *ordered, Basis *basis) {
 
 	Reflection *ref;
-	
+
 	ref = ordered->reflections;
-	
+
 	while ( ref ) {
-	
+
 		signed int h, k, l;
-		
+
 		h = ref->h;	k = ref->k;	l = ref->l;
-				
+
 		ref->x = h*basis->a.x + k*basis->b.x + l*basis->c.x;
 		ref->y = h*basis->a.y + k*basis->b.y + l*basis->c.y;
 		ref->z = h*basis->a.z + k*basis->b.z + l*basis->c.z;
-		
+
 		ref = ref->next;
 
 	}
@@ -296,7 +300,7 @@ int reflection_is_easy(Reflection *reflection) {
 	if ( reflection->h ) return !(reflection->k || reflection->l);
 	if ( reflection->k ) return !(reflection->h || reflection->l);
 	if ( reflection->l ) return !(reflection->h || reflection->k);
-	
+
 	return 0;	/* 000 */
 
 }
@@ -304,7 +308,7 @@ int reflection_is_easy(Reflection *reflection) {
 Reflection *reflectionlist_find(ReflectionList *reflectionlist, signed int h, signed int k, signed int l) {
 
 	Reflection *reflection;
-	
+
 	reflection = reflectionlist->reflections;
 	while ( reflection ) {
 		if ( (reflection->h==h) && (reflection->k==k) && (reflection->l==l) ) {
@@ -312,8 +316,7 @@ Reflection *reflectionlist_find(ReflectionList *reflectionlist, signed int h, si
 		}
 		reflection = reflection->next;
 	};
-	
+
 	return NULL;
 
 }
-
diff --git a/src/reproject.c b/src/reproject.c
index df85f4f..e0199af 100644
--- a/src/reproject.c
+++ b/src/reproject.c
@@ -18,6 +18,7 @@
 #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) {
@@ -37,6 +38,8 @@ void reproject_partner_features(ImageFeatureList *rflist, ImageRecord *image) {
 		if ( (d <= 20.0) && partner ) {
 			rflist->features[i].partner = partner;
 			rflist->features[i].partner_d = d;
+			image->features->features[idx].partner =
+				&rflist->features[i];
 		} else {
 			rflist->features[i].partner = NULL;
 		}
@@ -45,16 +48,19 @@ void reproject_partner_features(ImageFeatureList *rflist, ImageRecord *image) {
 
 }
 
-ImageFeatureList *reproject_get_reflections(ImageRecord *image, ReflectionList *reflectionlist) {
-
+ImageFeatureList *reproject_get_reflections(ImageRecord *image,
+						ReflectionList *reflectionlist)
+{
 	ImageFeatureList *flist;
 	Reflection *reflection;
-	double smax = 0.4e9;
+	double smax = 0.2e9;
 	double tilt, omega, k;
-	double nx, ny, nz;			/* "normal" vector */
-	double kx, ky, kz;			/* Electron wavevector ("normal" times 1/lambda */
-	double ux, uy, uz, uxt, uyt, uzt;	/* "up" vector (and calculation intermediates) */
-	//int done_debug = 0;
+	double nx, ny, nz; /* "normal" vector */
+	double kx, ky, kz; /* Electron wavevector ("normal" times 1/lambda */
+	double ux, uy, uz; /* "up" vector */
+	double rx, ry, rz; /* "right" vector */
+	int nrr = 0;
+	int n = 0;
 
 	flist = image_feature_list_new();
 
@@ -67,15 +73,12 @@ ImageFeatureList *reproject_get_reflections(ImageRecord *image, ReflectionList *
 	kx = nx / image->lambda;
 	ky = ny / image->lambda;
 	kz = nz / image->lambda;	/* This is the centre of the Ewald sphere */
-	//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);
+	/* 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 = reflectionlist->reflections;
 	do {
@@ -92,7 +95,8 @@ ImageFeatureList *reproject_get_reflections(ImageRecord *image, ReflectionList *
 		g_sq = modulus_squared(xl, yl, zl);
 		gn = xl*nx + yl*ny + zl*nz;
 
-		/* Next, solve the relrod equation to calculate the excitation error */
+		/* 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;
@@ -106,10 +110,8 @@ ImageFeatureList *reproject_get_reflections(ImageRecord *image, ReflectionList *
 
 			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;
@@ -117,64 +119,36 @@ ImageFeatureList *reproject_get_reflections(ImageRecord *image, ReflectionList *
 			/* 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 */
+			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 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");
-				}
+				double dx, dy, psi;
 
-			//	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));
+				/* 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);
 
-				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 */
+				/* Get image coordinates from polar
+				 * representation */
 				if ( image->fmode == FORMULATION_CLEN ) {
 					x = image->camera_len*tan(theta)*cos(psi);
 					y = image->camera_len*tan(theta)*sin(psi);
 					x *= image->resolution;
 					y *= image->resolution;
-				} else if ( image->fmode == FORMULATION_PIXELSIZE ) {
+				} else if ( image->fmode==FORMULATION_PIXELSIZE ) {
 					x = tan(theta)*cos(psi) / image->lambda;
 					y = tan(theta)*sin(psi) / image->lambda;
 					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;
 				}
 
@@ -182,31 +156,45 @@ ImageFeatureList *reproject_get_reflections(ImageRecord *image, ReflectionList *
 				y += image->y_centre;
 
 				/* Sanity check */
-				if ( (x>=0) && (x<image->width) && (y>=0) && (y<image->height) ) {
-
-					/* 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);
+				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);
+					printf("adding %i %i %i\n",
+						reflection->h,
+						reflection->k,
+						reflection->l);
+					n++;
 
 				} /* else it's outside the picture somewhere */
 
-			} /* else this is the central beam so don't worry about it */
+			} /* else this is the central beam */
 
 		}
 
 		reflection = reflection->next;
 
+		nrr++;
+
 	} 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. */
+	/* Partner features only if the image has a feature list.  This allows
+	 *  the test programs to use this function to generate simulated data
+	 */
 	if ( image->features != NULL ) {
 		reproject_partner_features(flist, image);
 	}
 
+	printf("processed %i, found %i\n", nrr, n);
+
 	return flist;
 
 }
@@ -226,8 +214,20 @@ void reproject_cell_to_lattice(ControlContext *ctx) {
 
 	/* 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;
+
+		ImageRecord *image;
+		int j;
+
+		image = &ctx->images->images[i];
+		if ( image->rflist != NULL ) {
+			image_feature_list_free(image->rflist);
+			image->rflist = NULL;
+		}
+
+		for ( j=0; j<image->features->n_features; j++ ) {
+			image->features->features[i].partner = NULL;
+		}
+
 	}
 
 }
@@ -235,7 +235,65 @@ void reproject_cell_to_lattice(ControlContext *ctx) {
 /* Notify that ctx->cell has changed (also need to call displaywindow_update) */
 void reproject_lattice_changed(ControlContext *ctx) {
 
+	int i;
+	int total_reprojected = 0;
+	int total_measured = 0;
+	int partnered_reprojected = 0;
+	int partnered_measured = 0;
+
 	reproject_cell_to_lattice(ctx);
-	displaywindow_update_imagestack(ctx->dw);
+	printf("%i reflections\n", ctx->cell_lattice->n_reflections);
 
+	for ( i=0; i<ctx->images->n_images; i++ ) {
+
+		ImageRecord *image;
+		int j;
+
+		image = &ctx->images->images[i];
+
+		/* Perform relrod projection */
+		image->rflist = reproject_get_reflections(image,
+							ctx->cell_lattice);
+
+		/* Loop over reprojected reflections */
+		for ( j=0; j<image->rflist->n_features; j++ ) {
+
+			double d;
+			ImageFeature f;
+			ImageFeature *p;
+
+			f = image->rflist->features[j];
+			p = image->rflist->features[j].partner;
+
+			if ( p != NULL ) {
+				d = distance(f.x, f.y, p->x, p->y);
+				partnered_reprojected++;
+			}
+		}
+		total_reprojected += image->rflist->n_features;
+
+		/* Loop over measured reflections */
+		for ( j=0; j<image->features->n_features; j++ ) {
+
+			double d;
+			ImageFeature f;
+			ImageFeature *p;
+
+			f = image->features->features[j];
+			p = image->features->features[j].partner;
+
+			if ( p != NULL ) {
+				d = distance(f.x, f.y, p->x, p->y);
+				partnered_measured++;
+			}
+		}
+		total_measured += image->features->n_features;
+	}
+	printf("%i images\n", ctx->images->n_images);
+	printf("%i measured reflections\n", total_measured);
+	printf("%i reprojected reflections\n", total_reprojected);
+	printf("%i partnered measured reflections\n", partnered_measured);
+	printf("%i partnered reprojected reflections\n", partnered_reprojected);
+
+	displaywindow_update_imagestack(ctx->dw);
 }