DirAx stuff

git-svn-id: svn://cook.msm.cam.ac.uk:745/diff-tomo/dtr@139 bf6ca9ba-c028-0410-8290-897cf20841d1

1 files changed, 437 insertions, 0 deletions

diff --git a/src/dirax.c b/src/dirax.c
new file mode 100644
index 0000000..4f4d135
--- /dev/null
+++ b/src/dirax.c
@@ -0,0 +1,437 @@
+/*
+ * dirax.c
+ *
+ * Invoke the DirAx auto-indexing program
+ * also: handle DirAx input files
+ *
+ * (c) 2007 Thomas White <taw27@cam.ac.uk>
+ *
+ *  dtr - Diffraction Tomography Reconstruction
+ *
+ */
+
+#ifdef HAVE_CONFIG_H
+#include <config.h>
+#endif
+
+#include <glib.h>
+#include <stdlib.h>
+#include <stdio.h>
+#include <math.h>
+#include <string.h>
+#include <pty.h>
+#include <unistd.h>
+#include <sys/wait.h>
+#include <fcntl.h>
+#include <assert.h>
+#include <sys/ioctl.h>
+#include <termio.h>
+	
+#include "control.h"
+#include "reflections.h"
+#include "utils.h"
+#include "basis.h"
+#include "displaywindow.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;
+	char *copy;
+	
+	copy = strdup(line);
+	for ( i=0; i<strlen(copy); i++ ) {
+		if ( copy[i] == '\r' ) copy[i]='r';
+		if ( copy[i] == '\n' ) copy[i]='n';
+	}
+	printf("DX: DirAx: %s\n", copy);
+	free(copy);
+	
+	/* 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] == ' ')) ) {
+		if ( line[i] == 'R' ) rf = 1;
+		if ( (line[i] == 'D') && rf ) {
+			ctx->dirax_read_cell = 1;
+			if ( ctx->cell ) {
+				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;
+			return;
+		}
+		i++;
+	}
+	
+	/* Parse unit cell vectors as appropriate */
+	if ( ctx->dirax_read_cell == 1 ) {
+		/* First row of unit cell values */
+		float x, y, z;
+		sscanf(line, "%f %f %f", &x, &y, &z);
+		ctx->cell->a.x = x*1e10;  ctx->cell->a.y = y*1e10;  ctx->cell->a.z = z*1e10;
+		ctx->dirax_read_cell++;
+		return;
+	} else if ( ctx->dirax_read_cell == 2 ) {
+		/* First row of unit cell values */
+		float x, y, z;
+		sscanf(line, "%f %f %f", &x, &y, &z);
+		ctx->cell->b.x = x*1e10;  ctx->cell->b.y = y*1e10;  ctx->cell->b.z = z*1e10;
+		ctx->dirax_read_cell++;
+		return;
+	} else if ( ctx->dirax_read_cell == 3 ) {
+		/* First row of unit cell values */
+		float x, y, z;
+		sscanf(line, "%f %f %f", &x, &y, &z);
+		ctx->cell->c.x = x*1e10;  ctx->cell->c.y = y*1e10;  ctx->cell->c.z = z*1e10;
+		printf("DX: Read a reciprocal unit cell\n");
+		displaywindow_update(ctx->dw);
+		ctx->dirax_read_cell = 0;
+		return;
+	}
+	
+	ctx->dirax_read_cell = 0;
+
+}
+
+static void dirax_sendline(const char *line, ControlContext *ctx) {
+	
+	char *copy;
+	int i;
+	
+	write(ctx->dirax_pty, line, strlen(line));
+
+	copy = strdup(line);
+	for ( i=0; i<strlen(copy); i++ ) {
+		if ( copy[i] == '\r' ) copy[i]='\0';
+		if ( copy[i] == '\n' ) copy[i]='\0';
+	}
+	printf("sent '%s'", copy);	/* No newline here */
+	free(copy);
+	
+}
+
+static void dirax_send_next(ControlContext *ctx) {
+
+	switch ( ctx->dirax_step ) {
+	
+		case 0 : {
+			dirax_sendline("read dtr.drx\n", ctx);
+			ctx->dirax_step++;
+			break;
+		}
+		
+		case 1 : {
+			dirax_sendline("go\n", ctx);
+			ctx->dirax_step++;
+			break;
+		}
+		
+		case 2 : {
+			dirax_sendline("cell\n", ctx);
+			ctx->dirax_step++;
+			break;
+		}
+		
+		default: {
+			printf("idle");
+		}
+		
+	}
+	
+	fsync(ctx->dirax_pty);
+	
+}
+
+static gboolean dirax_readable(GIOChannel *dirax, GIOCondition condition, ControlContext *ctx) {
+	
+	int rval;
+	
+	rval = read(ctx->dirax_pty, ctx->dirax_rbuffer+ctx->dirax_rbufpos, ctx->dirax_rbuflen-ctx->dirax_rbufpos);
+	
+	if ( (rval == -1) || (rval == 0) ) {
+	
+		printf("DX: Lost connection to DirAx\n");
+		g_io_channel_shutdown(ctx->dirax, FALSE, NULL);
+		ctx->dirax = NULL;
+		return FALSE;
+		
+	} else {
+	
+		int no_string = 0;
+		
+		ctx->dirax_rbufpos += rval;
+		assert(ctx->dirax_rbufpos <= ctx->dirax_rbuflen);
+		
+		while ( (!no_string) && (ctx->dirax_rbufpos > 0) ) {
+	
+			int i;
+			int block_ready = 0;
+			DirAxInputType type = DIRAX_INPUT_NONE;
+			
+			/* See if there's a full line in the buffer yet */
+			for ( i=0; i<ctx->dirax_rbufpos; i++ ) {	/* Means the last value looked at is roffset-1 */
+			
+				if ( (ctx->dirax_rbuffer[i] == '\n') ) { 
+					block_ready = 1;
+					type = DIRAX_INPUT_LINE;
+					break;
+				}
+				
+				if ( i <= ctx->dirax_rbufpos-7 ) {
+					if ( strncmp(ctx->dirax_rbuffer+i, "Dirax> ", 7) == 0 ) {
+						block_ready = 1;
+						type = DIRAX_INPUT_PROMPT;
+						break;
+					}
+				}
+			
+			}
+			
+			if ( block_ready == 1 ) {
+			
+				switch ( type ) {
+				
+					case DIRAX_INPUT_LINE : {
+					
+						char *block_buffer = NULL;
+						unsigned int new_rbuflen;
+						unsigned int endbit_length;
+						
+						block_buffer = malloc(i+1);
+						memcpy(block_buffer, ctx->dirax_rbuffer, i);
+						block_buffer[i] = '\0';
+						dirax_parseline(block_buffer, ctx);
+						free(block_buffer);
+						
+						/* Now the block's been parsed, it should be forgotten about */
+						endbit_length = i+1;
+						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;
+						if ( new_rbuflen == 0 ) {
+							new_rbuflen = 256;
+						}
+						ctx->dirax_rbuffer = realloc(ctx->dirax_rbuffer, new_rbuflen);
+						ctx->dirax_rbuflen = new_rbuflen;
+						
+						break;
+						
+					}
+					
+					case DIRAX_INPUT_PROMPT : {
+						
+						memmove(ctx->dirax_rbuffer+i, ctx->dirax_rbuffer+i+7, 7);
+						ctx->dirax_rbufpos -= 7;
+						printf("DX: Got prompt: ");
+						dirax_send_next(ctx);
+						printf("\n");
+						
+						break;
+					
+					}
+					
+					default : {
+						printf("DX: Unrecognised input mode (this never happens!)\n");
+						abort();
+					}
+				
+				}
+			
+			} else {
+			
+				if ( ctx->dirax_rbufpos == ctx->dirax_rbuflen ) {
+			
+					/* More buffer space is needed */
+					ctx->dirax_rbuffer = realloc(ctx->dirax_rbuffer, ctx->dirax_rbuflen + 256);
+					ctx->dirax_rbuflen = ctx->dirax_rbuflen + 256;
+					/* The new space gets used at the next read, shortly... */
+			
+				}
+				no_string = 1;
+			
+			}
+			
+		}
+		
+	}
+	
+	return TRUE;
+	
+}
+
+void dirax_invoke(ControlContext *ctx) {
+
+	FILE *fh;
+	Reflection *ref;
+	int pty;
+	unsigned int opts;
+	pid_t dirax_pid;
+	
+	printf("DX: Starting DirAx...\n");
+	
+	fh = fopen("dtr.drx", "w");
+	if ( !fh ) {
+		printf("DX: Couldn't open temporary file dtr.drx\n");
+		return;
+	}
+	fprintf(fh, "%f\n", 0.5);  /* Lie about the wavelength.  DirAx can't handle the truth. */
+	ref = ctx->reflectionlist->reflections;
+	while ( ref ) {
+		fprintf(fh, "%f %f %f %f\n", ref->x/1e10, ref->y/1e10, ref->z/1e10, ref->intensity);
+		ref = ref->next;
+	}
+	fclose(fh);
+	
+	struct termios tty;
+	if ( ioctl(0, TCGETS, &tty) ) {
+		printf("DX: Failed to read pty attributes\n");
+		return;
+	}
+	tty.c_lflag &= ~ECHO;
+	dirax_pid = forkpty(&pty, NULL, &tty, NULL);
+	if ( dirax_pid == -1 ) {
+		printf("DX: Failed to fork.\n");
+		return;
+	}
+	if ( dirax_pid == 0 ) {
+		
+		/* Child process: invoke DirAx */
+		execlp("dirax", "", (char *)NULL);
+		printf("(from the Other Side) Failed to invoke DirAx.\n");
+		_exit(0);
+		
+	}
+	
+	ctx->dirax_pty = pty;
+	ctx->dirax_rbuffer = malloc(256);
+	ctx->dirax_rbuflen = 256;
+	ctx->dirax_rbufpos = 0;
+	
+	/* Set non-blocking */
+	opts = fcntl(pty, F_GETFL);
+	fcntl(pty, F_SETFL, opts | O_NONBLOCK);
+	
+	ctx->dirax_step = 0;
+	ctx->dirax_read_cell = 0;
+	ctx->dirax = g_io_channel_unix_new(pty);
+	g_io_add_watch(ctx->dirax, G_IO_IN | G_IO_HUP, (GIOFunc)dirax_readable, ctx);
+	
+	return;
+	
+}
+
+/* Despite being part of the same module, this has very little to do with invoking DirAx */
+ReflectionList *dirax_load(const char *filename) {
+
+	FILE *fh;
+	char line[256];
+	ReflectionList *list;
+	int lambda_set = 0;
+	
+	fh = fopen(filename, "r");
+	if ( !fh ) {
+		printf("Couldn't open file '%s'\n", filename);
+		return 0;
+	}
+	
+	list = reflectionlist_new();
+	
+	while ( !feof(fh) ) {
+	
+		size_t ptr;
+		float lambda, theta, chib, phib;
+		
+		fgets(line, 255, fh);
+		ptr = skipspace(line);
+		if ( line[ptr] == '!' ) continue;
+		if ( line[ptr] == '\n' ) continue;
+		if ( line[ptr] == '\r' ) continue;
+		if ( sscanf(line+ptr, "%f %f %f\n", &theta, &phib, &chib) == 3 ) {
+			
+			double s, x, y, z;
+			float blah, intensity;
+			
+			/* Try to find an intensity value. Use dummy value if it fails */
+			if ( sscanf(line+ptr, "%f %f %f %f\n", &blah, &blah, &blah, &intensity) != 4 ) {
+				intensity = 1.0;
+			}
+				
+			if ( !lambda_set ) {
+				printf("DX: Wavelength not specified\n");
+				continue;
+			}
+			
+			chib = deg2rad(chib);
+			phib = deg2rad(phib);
+			theta = deg2rad(theta);
+			s = 2*(sin(theta)/lambda);
+			x = -s*cos(chib)*sin(phib);
+			y = +s*cos(chib)*cos(phib);
+			z = +s*sin(chib);
+			reflection_add(list, x, y, z, 1.0, REFLECTION_NORMAL);
+			
+			continue;
+			
+		}
+		if ( sscanf(line+ptr, "%f\n", &lambda) == 1 ) {
+			if ( lambda_set ) {
+				printf("DX: Warning: Found something which looks like a second wavelength\n");
+			}
+			lambda /= 1e10; /* Convert from A to m */
+			lambda_set = 1;
+		}
+		
+	}
+	
+	fclose(fh);
+	
+	return list;
+	
+}
+
+int dirax_is_drxfile(const char *filename) {
+
+	FILE *fh;
+	float lambda;
+	char line[256];
+	
+	fh = fopen(filename, "r");
+	if ( !fh ) {
+		printf("Couldn't open file '%s'\n", filename);
+		return 0;
+	}
+	
+	while ( !feof(fh) ) {
+	
+		size_t ptr;
+		
+		fgets(line, 255, fh);
+		ptr = skipspace(line);
+		if ( line[ptr] == '!' ) continue;
+		if ( line[ptr] == '\n' ) continue;
+		if ( line[ptr] == '\r' ) continue;
+		fscanf(fh, "%f\n", &lambda);
+		fclose(fh);
+		if ( lambda > 0.5 ) {
+			return 1;
+		} else {
+			return 0;
+		}
+		
+	}
+	
+	fclose(fh);
+	
+	return 0;
+	
+}
+