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