From a74ca2086dd001f243ea72f552ea73e5d0d07587 Mon Sep 17 00:00:00 2001 From: Thomas White Date: Wed, 18 Mar 2009 15:11:36 +0000 Subject: Linewrap, copyright and brackets in dirax.c --- src/dirax.c | 338 +++++++++++++++++++++++++++++++++--------------------------- 1 file changed, 188 insertions(+), 150 deletions(-) diff --git a/src/dirax.c b/src/dirax.c index e58983a..439d03b 100644 --- a/src/dirax.c +++ b/src/dirax.c @@ -4,7 +4,7 @@ * Invoke the DirAx auto-indexing program * also: handle DirAx input files * - * (c) 2007 Thomas White + * (c) 2007-2009 Thomas White * * dtr - Diffraction Tomography Reconstruction * @@ -26,7 +26,7 @@ #include #include #include -#include +#include #include "control.h" #include "reflections.h" @@ -41,11 +41,11 @@ typedef enum { DIRAX_INPUT_PROMPT } DirAxInputType; -static void dirax_parseline(const char *line, ControlContext *ctx) { - +static void dirax_parseline(const char *line, ControlContext *ctx) +{ int i, rf; char *copy; - + copy = strdup(line); for ( i=0; idw); + displaywindow_error("DirAx can't understand this data.", + ctx->dw); return; } - + /* Is this the first line of a unit cell specification? */ rf = 0; i = 0; - while ( (idirax_read_cell = 1; @@ -69,52 +71,63 @@ static void dirax_parseline(const char *line, ControlContext *ctx) { free(ctx->cell); } ctx->cell = malloc(sizeof(Basis)); - ctx->cell->a.x = 0.0; ctx->cell->a.y = 0.0; ctx->cell->a.z = 0.0; - ctx->cell->b.x = 0.0; ctx->cell->b.y = 0.0; ctx->cell->b.z = 0.0; - ctx->cell->c.x = 0.0; ctx->cell->c.y = 0.0; ctx->cell->c.z = 0.0; + ctx->cell->a.x = 0.0; + ctx->cell->a.y = 0.0; + ctx->cell->a.z = 0.0; + ctx->cell->b.x = 0.0; + ctx->cell->b.y = 0.0; + ctx->cell->b.z = 0.0; + ctx->cell->c.x = 0.0; + ctx->cell->c.y = 0.0; + ctx->cell->c.z = 0.0; return; } i++; } - + /* Parse unit cell vectors as appropriate */ if ( ctx->dirax_read_cell == 1 ) { /* First row of unit cell values */ float x1, x2, x3; sscanf(line, "%f %f %f", &x1, &x2, &x3); - ctx->cell->a.x = x1*1e10; ctx->cell->b.x = x2*1e10; ctx->cell->c.x = x3*1e10; + ctx->cell->a.x = x1*1e10; + ctx->cell->b.x = x2*1e10; + ctx->cell->c.x = x3*1e10; ctx->dirax_read_cell++; return; } else if ( ctx->dirax_read_cell == 2 ) { /* First row of unit cell values */ float y1, y2, y3; sscanf(line, "%f %f %f", &y1, &y2, &y3); - ctx->cell->a.y = y1*1e10; ctx->cell->b.y = y2*1e10; ctx->cell->c.y = y3*1e10; + ctx->cell->a.y = y1*1e10; + ctx->cell->b.y = y2*1e10; + ctx->cell->c.y = y3*1e10; ctx->dirax_read_cell++; return; } else if ( ctx->dirax_read_cell == 3 ) { /* First row of unit cell values */ float z1, z2, z3; sscanf(line, "%f %f %f", &z1, &z2, &z3); - ctx->cell->a.z = z1*1e10; ctx->cell->b.z = z2*1e10; ctx->cell->c.z = z3*1e10; + ctx->cell->a.z = z1*1e10; + ctx->cell->b.z = z2*1e10; + ctx->cell->c.z = z3*1e10; printf("DX: Read a reciprocal unit cell\n"); displaywindow_update(ctx->dw); reproject_lattice_changed(ctx); ctx->dirax_read_cell = 0; return; } - - ctx->dirax_read_cell = 0; + ctx->dirax_read_cell = 0; } -static void dirax_sendline(const char *line, ControlContext *ctx) { - +static void dirax_sendline(const char *line, ControlContext *ctx) +{ char *copy; int i; - + write(ctx->dirax_pty, line, strlen(line)); - + copy = strdup(line); for ( i=0; idirax_step != 0 ) { printf("DX: DirAx not idle\n"); return; } - - dirax_sendline(line, ctx); + dirax_sendline(line, ctx); } -static void dirax_send_next(ControlContext *ctx) { - +static void dirax_send_next(ControlContext *ctx) +{ switch ( ctx->dirax_step ) { - + case 1 : { dirax_sendline("\\echo off\n", ctx); ctx->dirax_step++; break; } - + case 2 : { dirax_sendline("read dtr.drx\n", ctx); ctx->dirax_step++; break; } - + case 3 : { dirax_sendline("dmax 10\n", ctx); ctx->dirax_step++; break; } - + case 4 : { dirax_sendline("indexfit 6\n", ctx); ctx->dirax_step++; break; } - + case 5 : { dirax_sendline("levelfit 300\n", ctx); ctx->dirax_step++; break; } - + case 6 : { dirax_sendline("go\n", ctx); ctx->dirax_step++; break; } - + case 7 : { dirax_sendline("cell\n", ctx); ctx->dirax_step++; break; } - + default: { ctx->dirax_step = 0; printf("DX: Prompt. DirAx is idle\n"); } - + } - } -static gboolean dirax_readable(GIOChannel *dirax, GIOCondition condition, ControlContext *ctx) { - +static gboolean dirax_readable(GIOChannel *dirax, GIOCondition condition, + ControlContext *ctx) +{ int rval; - - rval = read(ctx->dirax_pty, ctx->dirax_rbuffer+ctx->dirax_rbufpos, ctx->dirax_rbuflen-ctx->dirax_rbufpos); - + + rval = read(ctx->dirax_pty, ctx->dirax_rbuffer+ctx->dirax_rbufpos, + ctx->dirax_rbuflen-ctx->dirax_rbufpos); + if ( (rval == -1) || (rval == 0) ) { - + printf("DX: Lost connection to DirAx\n"); waitpid(ctx->dirax_pid, NULL, 0); g_io_channel_shutdown(ctx->dirax, FALSE, NULL); ctx->dirax = NULL; displaywindow_update_dirax(ctx, ctx->dw); return FALSE; - + } else { - + int no_string = 0; - + ctx->dirax_rbufpos += rval; assert(ctx->dirax_rbufpos <= ctx->dirax_rbuflen); - + while ( (!no_string) && (ctx->dirax_rbufpos > 0) ) { - + int i; int block_ready = 0; DirAxInputType type = DIRAX_INPUT_NONE; - + /* See if there's a full line in the buffer yet */ - for ( i=0; idirax_rbufpos-1; i++ ) { /* Means the last value looked at is rbufpos-2 */ - + for ( i=0; idirax_rbufpos-1; i++ ) { + /* Means the last value looked at is rbufpos-2 */ + /* Is there a prompt in the buffer? */ if ( i+7 <= ctx->dirax_rbufpos ) { - if ( (strncmp(ctx->dirax_rbuffer+i, "Dirax> ", 7) == 0) - || (strncmp(ctx->dirax_rbuffer+i, "PROMPT:", 7) == 0) ) { + if ( (strncmp(ctx->dirax_rbuffer+i, + "Dirax> ", 7) == 0) + || (strncmp(ctx->dirax_rbuffer+i, + "PROMPT:", 7) == 0) ) { block_ready = 1; type = DIRAX_INPUT_PROMPT; break; } } - - if ( (ctx->dirax_rbuffer[i] == '\r') && (ctx->dirax_rbuffer[i+1] == '\n') ) { + + if ( (ctx->dirax_rbuffer[i] == '\r') + && (ctx->dirax_rbuffer[i+1] == '\n') ) { block_ready = 1; type = DIRAX_INPUT_LINE; break; } - + } - + if ( block_ready ) { - + unsigned int new_rbuflen; unsigned int endbit_length; - + switch ( type ) { - + case DIRAX_INPUT_LINE : { - + char *block_buffer = NULL; - + block_buffer = malloc(i+1); - memcpy(block_buffer, ctx->dirax_rbuffer, i); + memcpy(block_buffer, + ctx->dirax_rbuffer, i); block_buffer[i] = '\0'; - + if ( block_buffer[0] == '\r' ) { - memmove(block_buffer, block_buffer+1, i); + memmove(block_buffer, + block_buffer+1, i); } - - dirax_parseline(block_buffer, ctx); + + dirax_parseline(block_buffer, + ctx); free(block_buffer); endbit_length = i+2; - + break; - + } - + case DIRAX_INPUT_PROMPT : { - + dirax_send_next(ctx); endbit_length = i+7; break; - + } - + default : { - printf("DX: Unrecognised input mode (this never happens!)\n"); + printf("DX: Unrecognised input " + "mode (this never happens!)\n"); abort(); } - + } - - /* Now the block's been parsed, it should be forgotten about */ - memmove(ctx->dirax_rbuffer, ctx->dirax_rbuffer + endbit_length, ctx->dirax_rbuflen - endbit_length); - + + /* Now the block's been parsed, it should + * be forgotten about */ + memmove(ctx->dirax_rbuffer, ctx->dirax_rbuffer + + endbit_length, ctx->dirax_rbuflen + - endbit_length); + /* Subtract the number of bytes removed */ - ctx->dirax_rbufpos = ctx->dirax_rbufpos - endbit_length; - new_rbuflen = ctx->dirax_rbuflen - endbit_length; + ctx->dirax_rbufpos = ctx->dirax_rbufpos + - endbit_length; + new_rbuflen = ctx->dirax_rbuflen-endbit_length; if ( new_rbuflen == 0 ) { new_rbuflen = 256; } - ctx->dirax_rbuffer = realloc(ctx->dirax_rbuffer, new_rbuflen); + ctx->dirax_rbuffer = realloc(ctx->dirax_rbuffer, + new_rbuflen); ctx->dirax_rbuflen = new_rbuflen; - + } else { - - if ( ctx->dirax_rbufpos == ctx->dirax_rbuflen ) { - + + if ( ctx->dirax_rbufpos==ctx->dirax_rbuflen ) { + /* More buffer space is needed */ - ctx->dirax_rbuffer = realloc(ctx->dirax_rbuffer, ctx->dirax_rbuflen + 256); - ctx->dirax_rbuflen = ctx->dirax_rbuflen + 256; - /* The new space gets used at the next read, shortly... */ - + ctx->dirax_rbuffer = realloc( + ctx->dirax_rbuffer, + ctx->dirax_rbuflen + 256); + ctx->dirax_rbuflen = ctx->dirax_rbuflen + + 256; + /* The new space gets used at the next + * read, shortly... */ + } no_string = 1; - + } - + } - + } - + return TRUE; - + } -void dirax_stop(ControlContext *ctx) { +void dirax_stop(ControlContext *ctx) +{ dirax_sendline_if_idle("end\n", ctx); } -void dirax_rerun(ControlContext *ctx) { +void dirax_rerun(ControlContext *ctx) +{ dirax_sendline_if_idle("go\n", ctx); ctx->dirax_step = 7; } -void dirax_invoke(ControlContext *ctx) { - +void dirax_invoke(ControlContext *ctx) +{ FILE *fh; Reflection *ref; unsigned int opts; int saved_stderr; - + if ( ctx->dirax ) { dirax_rerun(ctx); return; } - + printf("DX: Starting DirAx...\n"); - + fh = fopen("dtr.drx", "w"); if ( !fh ) { printf("DX: Couldn't open temporary file dtr.drx\n"); return; } - fprintf(fh, "%f\n", 0.5); /* Lie about the wavelength. DirAx can't handle the truth. */ + fprintf(fh, "%f\n", 0.5); /* Lie about the wavelength. DirAx can't + * handle the truth. */ ref = ctx->reflectionlist->reflections; while ( ref ) { - fprintf(fh, "%10f %10f %10f %8f\n", ref->x/1e10, ref->y/1e10, ref->z/1e10, ref->intensity); + fprintf(fh, "%10f %10f %10f %8f\n", ref->x/1e10, ref->y/1e10, + ref->z/1e10, ref->intensity); ref = ref->next; } fclose(fh); - + saved_stderr = dup(STDERR_FILENO); ctx->dirax_pid = forkpty(&ctx->dirax_pty, NULL, NULL, NULL); if ( ctx->dirax_pid == -1 ) { @@ -359,28 +392,28 @@ void dirax_invoke(ControlContext *ctx) { return; } if ( ctx->dirax_pid == 0 ) { - + /* Child process: invoke DirAx */ struct termios t; - + /* Turn echo off */ tcgetattr(STDIN_FILENO, &t); t.c_lflag &= ~(ECHO | ECHOE | ECHOK | ECHONL); tcsetattr(STDIN_FILENO, TCSANOW, &t); - + /* Reconnect stderr */ dup2(saved_stderr, STDERR_FILENO); - + execlp("dirax", "", (char *)NULL); printf("(from the Other Side) Failed to invoke DirAx.\n"); _exit(0); - + } - + ctx->dirax_rbuffer = malloc(256); ctx->dirax_rbuflen = 256; ctx->dirax_rbufpos = 0; - + /* Set non-blocking */ opts = fcntl(ctx->dirax_pty, F_GETFL); fcntl(ctx->dirax_pty, F_SETFL, opts | O_NONBLOCK); @@ -388,55 +421,60 @@ void dirax_invoke(ControlContext *ctx) { ctx->dirax_step = 1; /* This starts the "initialisation" procedure */ ctx->dirax_read_cell = 0; ctx->dirax = g_io_channel_unix_new(ctx->dirax_pty); - g_io_add_watch(ctx->dirax, G_IO_IN | G_IO_HUP, (GIOFunc)dirax_readable, ctx); - + g_io_add_watch(ctx->dirax, G_IO_IN | G_IO_HUP, (GIOFunc)dirax_readable, + ctx); + displaywindow_update_dirax(ctx, ctx->dw); - + return; - -} -/* Despite being part of the same module, this has very little to do with invoking DirAx */ -ReflectionList *dirax_load(const char *filename) { +} +/* Despite being part of the same module, this has very little to do with + * invoking DirAx */ +ReflectionList *dirax_load(const char *filename) +{ FILE *fh; char line[256]; ReflectionList *list; int lambda_set = 0; - + fh = fopen(filename, "r"); if ( !fh ) { printf("Couldn't open file '%s'\n", filename); return 0; } - + list = reflectionlist_new(); - + while ( !feof(fh) ) { - + size_t ptr; float lambda, theta, chib, phib; - + fgets(line, 255, fh); ptr = skipspace(line); if ( line[ptr] == '!' ) continue; if ( line[ptr] == '\n' ) continue; if ( line[ptr] == '\r' ) continue; - if ( sscanf(line+ptr, "%f %f %f\n", &theta, &phib, &chib) == 3 ) { - + if ( sscanf(line+ptr, "%f %f %f\n", + &theta, &phib, &chib) == 3 ) { + double s, x, y, z; float blah, intensity; - - /* Try to find an intensity value. Use dummy value if it fails */ - if ( sscanf(line+ptr, "%f %f %f %f\n", &blah, &blah, &blah, &intensity) != 4 ) { + + /* Try to find an intensity value. Use dummy value if + * it fails */ + if ( sscanf(line+ptr, "%f %f %f %f\n", &blah, &blah, + &blah, &intensity) != 4 ) { intensity = 1.0; } - + if ( !lambda_set ) { printf("DX: Wavelength not specified\n"); continue; } - + chib = deg2rad(chib); phib = deg2rad(phib); theta = deg2rad(theta); @@ -445,42 +483,43 @@ ReflectionList *dirax_load(const char *filename) { y = +s*cos(chib)*cos(phib); z = +s*sin(chib); reflection_add(list, x, y, z, 1.0, REFLECTION_NORMAL); - + continue; - + } if ( sscanf(line+ptr, "%f\n", &lambda) == 1 ) { if ( lambda_set ) { - printf("DX: Warning: Found something which looks like a second wavelength\n"); + printf("DX: Warning: Found something which " + "looks like a second wavelength\n"); } lambda /= 1e10; /* Convert from A to m */ lambda_set = 1; } - + } - + fclose(fh); - + return list; - } -int dirax_is_drxfile(const char *filename) { +int dirax_is_drxfile(const char *filename) +{ FILE *fh; float lambda; char line[256]; - + fh = fopen(filename, "r"); if ( !fh ) { printf("Couldn't open file '%s'\n", filename); return 0; } - + while ( !feof(fh) ) { - + size_t ptr; - + fgets(line, 255, fh); ptr = skipspace(line); if ( line[ptr] == '!' ) continue; @@ -493,12 +532,11 @@ int dirax_is_drxfile(const char *filename) { } else { return 0; } - + } - + fclose(fh); - + return 0; - + } - -- cgit v1.2.3 From 899919470fe997388b195b005b31dd29d7d2c06d Mon Sep 17 00:00:00 2001 From: Thomas White Date: Thu, 19 Mar 2009 12:06:53 +0000 Subject: Linewrap and conventions in src/dirax.c --- src/dirax.c | 15 +++++++++++---- 1 file changed, 11 insertions(+), 4 deletions(-) diff --git a/src/dirax.c b/src/dirax.c index 439d03b..2924380 100644 --- a/src/dirax.c +++ b/src/dirax.c @@ -35,12 +35,14 @@ #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; @@ -121,6 +123,7 @@ 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; @@ -136,6 +139,8 @@ static void dirax_sendline(const char *line, ControlContext *ctx) printf("DX: Sent '%s'\n", copy); /* No newline here */ free(copy); } + + /* Send a "user" command to DirAx, failing if DirAx is not idle */ static void dirax_sendline_if_idle(const char *line, ControlContext *ctx) { @@ -147,6 +152,7 @@ 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 ) { @@ -201,6 +207,7 @@ static void dirax_send_next(ControlContext *ctx) } } + static gboolean dirax_readable(GIOChannel *dirax, GIOCondition condition, ControlContext *ctx) { @@ -345,11 +352,13 @@ static gboolean dirax_readable(GIOChannel *dirax, GIOCondition condition, } + void dirax_stop(ControlContext *ctx) { dirax_sendline_if_idle("end\n", ctx); } + void dirax_rerun(ControlContext *ctx) { dirax_sendline_if_idle("go\n", ctx); @@ -425,11 +434,9 @@ void dirax_invoke(ControlContext *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 */ ReflectionList *dirax_load(const char *filename) @@ -503,6 +510,7 @@ ReflectionList *dirax_load(const char *filename) return list; } + int dirax_is_drxfile(const char *filename) { @@ -538,5 +546,4 @@ int dirax_is_drxfile(const char *filename) fclose(fh); return 0; - } -- cgit v1.2.3 From dc76f055cefcfc70f97f4c5520e65eec7e662fcf Mon Sep 17 00:00:00 2001 From: Thomas White Date: Thu, 19 Mar 2009 12:07:53 +0000 Subject: Send a random selection of reflections to DirAx if more than 1,000 are present --- src/dirax.c | 67 ++++++++++++++++++++++++++++++++++++++++++++++++++++++++----- 1 file changed, 62 insertions(+), 5 deletions(-) diff --git a/src/dirax.c b/src/dirax.c index 2924380..7c79412 100644 --- a/src/dirax.c +++ b/src/dirax.c @@ -365,12 +365,59 @@ void dirax_rerun(ControlContext *ctx) ctx->dirax_step = 7; } + +static void dirax_send_random_selection(ReflectionList *r, int n, FILE *fh) +{ + char *used; + int i; + + used = malloc(n*sizeof(char)); + + for ( i=0; ireflections; + for ( j=0; jnext; + } + fprintf(fh, "%10f %10f %10f %8f\n", + ref->x/1e10, ref->y/1e10, + ref->z/1e10, ref->intensity); + used[ra] = 'U'; + + } +} + + void dirax_invoke(ControlContext *ctx) { FILE *fh; Reflection *ref; unsigned int opts; int saved_stderr; + int n; if ( ctx->dirax ) { dirax_rerun(ctx); @@ -386,11 +433,21 @@ void dirax_invoke(ControlContext *ctx) } fprintf(fh, "%f\n", 0.5); /* Lie about the wavelength. DirAx can't * handle the truth. */ - ref = ctx->reflectionlist->reflections; - while ( ref ) { - fprintf(fh, "%10f %10f %10f %8f\n", ref->x/1e10, ref->y/1e10, - ref->z/1e10, ref->intensity); - ref = ref->next; + + n = ctx->reflectionlist->n_reflections; + printf("DX: There are %i reflections - ", n); + if ( n > 1000 ) { + printf("sending a random selection to DirAx\n"); + dirax_send_random_selection(ctx->reflectionlist, n, fh); + } else { + printf("sending them all to DirAx\n"); + ref = ctx->reflectionlist->reflections; + while ( ref ) { + fprintf(fh, "%10f %10f %10f %8f\n", + ref->x/1e10, ref->y/1e10, + ref->z/1e10, ref->intensity); + ref = ref->next; + } } fclose(fh); -- cgit v1.2.3 From 2639321f2aa8e489f60b7cc60b4e4c89f990826e Mon Sep 17 00:00:00 2001 From: Thomas White Date: Fri, 20 Mar 2009 12:36:41 +0000 Subject: Limit the resolution of reflections passed to DirAx --- src/dirax.c | 28 +++++++++++++++++----------- 1 file changed, 17 insertions(+), 11 deletions(-) diff --git a/src/dirax.c b/src/dirax.c index 7c79412..2c82833 100644 --- a/src/dirax.c +++ b/src/dirax.c @@ -377,23 +377,20 @@ static void dirax_send_random_selection(ReflectionList *r, int n, FILE *fh) used[i] = '-'; } - for ( i=0; i<1000; i++ ) { + i = 0; + while ( i < 1000 ) { Reflection *ref; - int done; int j; long long int ra; - done = 0; - while ( !done ) { - ra = ((long long int)random() * (long long int)n); - ra /= RAND_MAX; - if ( used[ra] == 'U' ) { - printf("%lli - already used\n", ra); - } else { - done = 1; - } + ra = ((long long int)random() * (long long int)n); + ra /= RAND_MAX; + if ( used[ra] == 'U' ) { + printf("%lli - already used\n", ra); + continue; } + printf("Selected reflection %lli\n", ra); /* Dig out the correct reflection. A little faffy @@ -402,11 +399,20 @@ static void dirax_send_random_selection(ReflectionList *r, int n, FILE *fh) for ( j=0; jnext; } + + /* Limit resolution of reflections */ + if ( ( (ref->x/1e9)*(ref->x/1e9) + + (ref->y/1e9)*(ref->y/1e9) + + (ref->z/1e9)*(ref->z/1e9) ) > (20e9)*(20e9) ) + continue; + fprintf(fh, "%10f %10f %10f %8f\n", ref->x/1e10, ref->y/1e10, ref->z/1e10, ref->intensity); used[ra] = 'U'; + i++; + } } -- cgit v1.2.3 From db8113d6f55baf4e4ad9c2ecfa423418e73b37c3 Mon Sep 17 00:00:00 2001 From: Thomas White Date: Fri, 20 Mar 2009 12:36:52 +0000 Subject: levelfit=200 --- src/dirax.c | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/dirax.c b/src/dirax.c index 2c82833..d27cdac 100644 --- a/src/dirax.c +++ b/src/dirax.c @@ -182,7 +182,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; } -- cgit v1.2.3 From 69d1342278729353799d0fee646e0f44517ebed8 Mon Sep 17 00:00:00 2001 From: Thomas White Date: Sat, 28 Mar 2009 18:04:06 +0000 Subject: Loading and saving of cells, and general tidy-up --- data/displaywindow.ui | 3 + src/basis.c | 290 +++++++++++++++++++++++++++++++++----------------- src/basis.h | 13 ++- src/dirax.c | 172 ++++++++---------------------- src/displaywindow.c | 19 +++- src/intensities.c | 154 +++++++++++++-------------- 6 files changed, 339 insertions(+), 312 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 @@ + + + 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 + * (c) 2007-2009 Thomas White * * 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 + * (c) 2007-2009 Thomas White * * 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 d27cdac..19a6e26 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 ( (idirax_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; idirax_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,81 +338,25 @@ 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; - int i; - - used = malloc(n*sizeof(char)); - - for ( i=0; ireflections; - for ( j=0; jnext; - } - - /* Limit resolution of reflections */ - if ( ( (ref->x/1e9)*(ref->x/1e9) - + (ref->y/1e9)*(ref->y/1e9) - + (ref->z/1e9)*(ref->z/1e9) ) > (20e9)*(20e9) ) - continue; - - fprintf(fh, "%10f %10f %10f %8f\n", - ref->x/1e10, ref->y/1e10, - ref->z/1e10, ref->intensity); - used[ra] = 'U'; - - i++; - - } -} - - void dirax_invoke(ControlContext *ctx) { FILE *fh; Reflection *ref; unsigned int opts; int saved_stderr; - int n; if ( ctx->dirax ) { dirax_rerun(ctx); @@ -437,23 +370,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); @@ -494,11 +416,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 */ @@ -527,14 +450,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; @@ -571,12 +494,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 + * (c) 2007-2009 Thomas White * * 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; iimages->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; jrflist->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); - } - -- cgit v1.2.3 From 655a56776090b18f8bf66ed74c0fb23226513979 Mon Sep 17 00:00:00 2001 From: Thomas White Date: Mon, 30 Mar 2009 12:12:23 +0100 Subject: indexfit=2 --- src/dirax.c | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/dirax.c b/src/dirax.c index d27cdac..c4afbce 100644 --- a/src/dirax.c +++ b/src/dirax.c @@ -176,7 +176,7 @@ static void dirax_send_next(ControlContext *ctx) } case 4 : { - dirax_sendline("indexfit 6\n", ctx); + dirax_sendline("indexfit 2\n", ctx); ctx->dirax_step++; break; } -- cgit v1.2.3 From 3ffa0ce286bac2c3c4d9e7d41bd2ca8972d6288a Mon Sep 17 00:00:00 2001 From: Thomas White Date: Mon, 30 Mar 2009 12:14:17 +0100 Subject: Tidy up and remove debug --- src/dirax.c | 10 +++---- src/itrans-stat.c | 6 +--- src/reproject.c | 83 +++++++++++++++++++++++++++---------------------------- 3 files changed, 46 insertions(+), 53 deletions(-) diff --git a/src/dirax.c b/src/dirax.c index c4afbce..79b20ea 100644 --- a/src/dirax.c +++ b/src/dirax.c @@ -387,11 +387,9 @@ static void dirax_send_random_selection(ReflectionList *r, int n, FILE *fh) ra = ((long long int)random() * (long long int)n); ra /= RAND_MAX; if ( used[ra] == 'U' ) { - printf("%lli - already used\n", ra); continue; } - printf("Selected reflection %lli\n", ra); /* Dig out the correct reflection. A little faffy * because of the linked list */ @@ -401,10 +399,10 @@ static void dirax_send_random_selection(ReflectionList *r, int n, FILE *fh) } /* Limit resolution of reflections */ - if ( ( (ref->x/1e9)*(ref->x/1e9) - + (ref->y/1e9)*(ref->y/1e9) - + (ref->z/1e9)*(ref->z/1e9) ) > (20e9)*(20e9) ) - continue; +// if ( ( (ref->x/1e9)*(ref->x/1e9) +// + (ref->y/1e9)*(ref->y/1e9) +// + (ref->z/1e9)*(ref->z/1e9) ) > (20e9)*(20e9) ) +// continue; fprintf(fh, "%10f %10f %10f %8f\n", ref->x/1e10, ref->y/1e10, diff --git a/src/itrans-stat.c b/src/itrans-stat.c index ee41817..0e03ead 100644 --- a/src/itrans-stat.c +++ b/src/itrans-stat.c @@ -441,7 +441,7 @@ static gsl_matrix *itrans_peaksearch_stat_floodfill(gsl_matrix *m, gsl_matrix_set(p, 2, n, val); n++; if ( n == size ) { - printf("expanding %i->%i\n", size, size*2); + p = itrans_peaksearch_stat_matrix_expand(p, size, size*2); size *= 2; } @@ -450,14 +450,10 @@ static gsl_matrix *itrans_peaksearch_stat_floodfill(gsl_matrix *m, } } } - //printf("ff: ending loop, found %d\n",n); *count = n; if ( n > 0 ) { - //printf("pcheck s1=%d s2=%d\n",p->size1,p->size2); - printf("expandingg %i->%i\n", size, n); p = itrans_peaksearch_stat_matrix_expand(p, size, n); - //printf("pcheck s1=%d s2=%d\n",p->size1,p->size2); } return p; diff --git a/src/reproject.c b/src/reproject.c index e1ce13e..c41aa3f 100644 --- a/src/reproject.c +++ b/src/reproject.c @@ -23,24 +23,24 @@ void reproject_partner_features(ImageFeatureList *rflist, ImageRecord *image) { int i; - + for ( i=0; in_features; i++ ) { - + //if ( !reflection_is_easy(rflist->features[i].reflection) ) continue; - + double d = 0.0; ImageFeature *partner; int idx; - + partner = image_feature_closest(image->features, rflist->features[i].x, rflist->features[i].y, &d, &idx); - + if ( (d <= 20.0) && partner ) { rflist->features[i].partner = partner; rflist->features[i].partner_d = d; } else { rflist->features[i].partner = NULL; } - + } } @@ -55,12 +55,12 @@ ImageFeatureList *reproject_get_reflections(ImageRecord *image, ReflectionList * 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; - + flist = image_feature_list_new(); - + tilt = image->tilt; omega = image->omega; - + /* Calculate the (normalised) incident electron wavevector * n is rotated "into" the reconstruction, so only one omega step. */ nxt = 0.0; nyt = 0.0; nzt = 1.0; @@ -71,7 +71,7 @@ ImageFeatureList *reproject_get_reflections(ImageRecord *image, ReflectionList * 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; @@ -79,19 +79,19 @@ ImageFeatureList *reproject_get_reflections(ImageRecord *image, ReflectionList * uxt = ux; uyt = uy; uzt = uz; ux = uxt*cos(-omega) + uyt*-sin(omega); uy = -uxt*sin(omega) + uyt*cos(omega); uz = uzt; //reflection_add(ctx->reflectionlist, ux*50e9, uy*50e9, uz*50e9, 1, REFLECTION_VECTOR_MARKER_2); - + reflection = reflectionlist->reflections; do { - + double xl, yl, zl; double a, b, c; double A1, A2, s1, s2, s; - + /* Get the coordinates of the reciprocal lattice point */ xl = reflection->x; yl = reflection->y; zl = reflection->z; - + /* Next, solve the relrod equation to calculate the excitation error */ a = 1.0; b = -2.0*(xl*nx + yl*ny + zl*nz); @@ -101,37 +101,37 @@ ImageFeatureList *reproject_get_reflections(ImageRecord *image, ReflectionList * s1 = 1.0/image->lambda - A1; s2 = 1.0/image->lambda - A2; if ( fabs(s1) < fabs(s2) ) s = s1; else s = s2; - + /* Skip this reflection if s is large */ if ( fabs(s) <= smax ) { - + double xi, yi, zi; double gx, gy, gz; double cx, cy, cz; double theta; double x, y; double rx, ry, rz; - + /* Determine the intersection point */ xi = xl + s*nx; yi = yl + s*ny; zi = zl + s*nz; - + /* Calculate Bragg angle */ gx = xi - kx; gy = yi - ky; gz = zi - kz; /* This is the vector from the centre of the sphere to the intersection */ theta = angle_between(-kx, -ky, -kz, gx, gy, gz); - + if ( theta > 0.0 ) { - + double 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); @@ -154,15 +154,15 @@ ImageFeatureList *reproject_get_reflections(ImageRecord *image, ReflectionList * psi = 3*M_PI_2 + psi; /* Case 4 */ // if ( (i==20) && !done_debug ) printf("case 4\n"); } - + // if ( (i==20) && !done_debug ) printf("final psi=%f clockwise from 'up'\n", rad2deg(psi)); psi = M_PI_2 - psi; /* Anticlockwise from "+x" instead of clockwise from "up" */ // if ( (i==20) && !done_debug ) printf("final psi=%f anticlockwise from +x\n", rad2deg(psi)); - + psi += omega; // if ( (i==20) && !done_debug ) printf("final psi=%f anticlockwise from +tilt axis\n", rad2deg(psi)); // if ( (i==20) && !done_debug ) done_debug = 1; - + /* Calculate image coordinates from polar representation */ if ( image->fmode == FORMULATION_CLEN ) { x = image->camera_len*tan(theta)*cos(psi); @@ -178,36 +178,36 @@ ImageFeatureList *reproject_get_reflections(ImageRecord *image, ReflectionList * fprintf(stderr, "Unrecognised formulation mode in reproject_get_reflections()\n"); return NULL; } - + x += image->x_centre; y += image->y_centre; - + /* Sanity check */ if ( (x>=0) && (xwidth) && (y>=0) && (yheight) ) { - + /* Record the reflection */ image_add_feature_reflection(flist, x, y, image, reflection->intensity, reflection); /* Intensity should be multiplied by relrod spike function except * reprojected reflections aren't used quantitatively for anything */ - + //printf("Reflection at %f, %f\n", x, y); - + } /* else it's outside the picture somewhere */ - + } /* else this is the central beam so don't worry about it */ - + } - + reflection = reflection->next; - + } while ( reflection ); - + /* Partner features only if the image has a feature list. This allows the test * program to use this function to generate simulated data. */ if ( image->features != NULL ) { reproject_partner_features(flist, image); } - + return flist; } @@ -216,21 +216,21 @@ ImageFeatureList *reproject_get_reflections(ImageRecord *image, ReflectionList * void reproject_cell_to_lattice(ControlContext *ctx) { int i; - + if ( ctx->cell_lattice ) { reflection_list_from_new_cell(ctx->cell_lattice, ctx->cell); } else { ctx->cell_lattice = reflection_list_from_cell(ctx->cell); } - + displaywindow_enable_cell_functions(ctx->dw, TRUE); - + /* Invalidate all the reprojected feature lists */ for ( i=0; iimages->n_images; i++ ) { image_feature_list_free(ctx->images->images[i].rflist); ctx->images->images[i].rflist = NULL; } - + } /* Notify that ctx->cell has changed (also need to call displaywindow_update) */ @@ -240,4 +240,3 @@ void reproject_lattice_changed(ControlContext *ctx) { displaywindow_update_imagestack(ctx->dw); } - -- cgit v1.2.3 From b97b7f02373229765038a181dfb03c15b4fc5b87 Mon Sep 17 00:00:00 2001 From: Thomas White Date: Mon, 30 Mar 2009 12:22:19 +0100 Subject: Restore chunk which got clobbered by 69d1342278729353799d0fee646e0f44517ebed8 --- src/dirax.c | 21 ++++++++++++++++----- 1 file changed, 16 insertions(+), 5 deletions(-) diff --git a/src/dirax.c b/src/dirax.c index 6e97d21..5063c60 100644 --- a/src/dirax.c +++ b/src/dirax.c @@ -405,6 +405,7 @@ void dirax_invoke(ControlContext *ctx) Reflection *ref; unsigned int opts; int saved_stderr; + int n; if ( ctx->dirax ) { dirax_rerun(ctx); @@ -419,11 +420,21 @@ void dirax_invoke(ControlContext *ctx) return; } 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; + + n = ctx->reflectionlist->n_reflections; + printf("DX: There are %i reflections - ", n); + if ( n > 1000 ) { + printf("sending a random selection to DirAx\n"); + dirax_send_random_selection(ctx->reflectionlist, n, fh); + } else { + printf("sending them all to DirAx\n"); + ref = ctx->reflectionlist->reflections; + while ( ref ) { + fprintf(fh, "%10f %10f %10f %8f\n", + ref->x/1e10, ref->y/1e10, + ref->z/1e10, ref->intensity); + ref = ref->next; + } } fclose(fh); -- cgit v1.2.3