aboutsummaryrefslogtreecommitdiff
path: root/src/basis.c
diff options
context:
space:
mode:
Diffstat (limited to 'src/basis.c')
-rw-r--r--src/basis.c290
1 files changed, 192 insertions, 98 deletions
diff --git a/src/basis.c b/src/basis.c
index 6e7aed8..37fedba 100644
--- a/src/basis.c
+++ b/src/basis.c
@@ -3,7 +3,7 @@
*
* Handle basis structures
*
- * (c) 2007 Thomas White <taw27@cam.ac.uk>
+ * (c) 2007-2009 Thomas White <taw27@cam.ac.uk>
*
* dtr - Diffraction Tomography Reconstruction
*
@@ -21,135 +21,229 @@
#include "reflections.h"
#include "basis.h"
#include "utils.h"
-
-double basis_efom(ReflectionList *reflectionlist, Basis *basis) {
-
- int n_indexed, n_counted;
- Reflection *cur;
-
- cur = reflectionlist->reflections;
- n_indexed = 0;
- n_counted = 0;
- while ( cur ) {
-
- if ( cur->type == REFLECTION_NORMAL ) {
-
- /* Can this basis "approximately" account for this reflection? */
- double det;
- double a11, a12, a13, a21, a22, a23, a31, a32, a33;
- double h, k, l;
-
- /* Set up the coordinate transform from hkl to xyz */
- a11 = basis->a.x; a12 = basis->a.y; a13 = basis->a.z;
- a21 = basis->b.x; a22 = basis->b.y; a23 = basis->b.z;
- a31 = basis->c.x; a32 = basis->c.y; a33 = basis->c.z;
-
- /* Invert the matrix to get hkl from xyz */
- det = a11*(a22*a33 - a23*a32) - a12*(a21*a33 - a23*a31) + a13*(a21*a32 - a22*a31);
- h = ((a22*a33-a23*a32)*cur->x + (a23*a31-a21*a33)*cur->y + (a21*a32-a22*a31)*cur->z) / det;
- k = ((a13*a32-a12*a33)*cur->x + (a11*a33-a13*a31)*cur->y + (a12*a31-a11*a32)*cur->z) / det;
- l = ((a12*a23-a13*a22)*cur->x + (a13*a21-a11*a23)*cur->y + (a11*a22-a12*a21)*cur->z) / det;
-
- /* Calculate the deviations in terms of |a|, |b| and |c| */
- h = fabs(h); k = fabs(k); l = fabs(l);
- h -= floor(h); k -= floor(k); l -= floor(l);
- if ( h == 1.0 ) h = 0.0;
- if ( k == 1.0 ) k = 0.0;
- if ( l == 1.0 ) l = 0.0;
-
- /* Define "approximately" here. Circle in basis space becomes an ellipsoid in reciprocal space */
- if ( h*h + k*k + l*l <= 0.1*0.1*0.1 ) n_indexed++;
-
- n_counted++;
-
- }
-
- cur = cur->next;
-
- }
-
- return (double)n_indexed / n_counted;
-
+#include "displaywindow.h"
+
+static void basis_print(Basis *cell)
+{
+ printf("%12.8f %12.8f %12.8f\n",
+ cell->a.x/1e9, cell->a.y/1e9, cell->a.z/1e9);
+ printf("%12.8f %12.8f %12.8f\n",
+ cell->b.x/1e9, cell->b.y/1e9, cell->b.z/1e9);
+ printf("%12.8f %12.8f %12.8f\n",
+ cell->c.x/1e9, cell->c.y/1e9, cell->c.z/1e9);
}
-Basis basis_add(Basis u, Basis v) {
-
- Basis ans;
-
- ans.a.x = u.a.x + v.a.x; ans.a.y = u.a.y + v.a.y; ans.a.z = u.a.z + v.a.z;
- ans.b.x = u.b.x + v.b.x; ans.b.y = u.b.y + v.b.y; ans.b.z = u.b.z + v.b.z;
- ans.c.x = u.c.x + v.c.x; ans.c.y = u.c.y + v.c.y; ans.c.z = u.c.z + v.c.z;
-
- return ans;
-
+static void cell_print(UnitCell *cell)
+{
+ printf("%12.8f %12.8f %12.8f nm\n",
+ cell->a*1e9, cell->b*1e9, cell->c*1e9);
+ printf("%12.8f %12.8f %12.8f deg\n",
+ rad2deg(cell->alpha), rad2deg(cell->beta), rad2deg(cell->gamma));
}
-static void basis_print(Basis *cell) {
-
- printf("%12.8f %12.8f %12.8f\n", cell->a.x/1e9, cell->a.y/1e9, cell->a.z/1e9);
- printf("%12.8f %12.8f %12.8f\n", cell->b.x/1e9, cell->b.y/1e9, cell->b.z/1e9);
- printf("%12.8f %12.8f %12.8f\n", cell->c.x/1e9, cell->c.y/1e9, cell->c.z/1e9);
-
-}
-
-static void cell_print(UnitCell *cell) {
-
- printf("%12.8f %12.8f %12.8f nm\n", cell->a*1e9, cell->b*1e9, cell->c*1e9);
- printf("%12.8f %12.8f %12.8f deg\n", rad2deg(cell->alpha), rad2deg(cell->beta), rad2deg(cell->gamma));
-
-}
-
-UnitCell basis_get_cell(Basis *basis) {
-
+UnitCell basis_get_cell(Basis *basis)
+{
UnitCell cell;
gsl_matrix *m;
gsl_matrix *inv;
gsl_permutation *perm;
double ax, ay, az, bx, by, bz, cx, cy, cz;
int s;
-
-// basis->a.x = 0.5; basis->a.y = 0.0; basis->a.z = 0.0;
-// basis->b.x = 0.0; basis->b.y = 0.2; basis->b.z = 0.0;
-// basis->c.x = 0.0; basis->c.y = 0.0; basis->c.z = 1.0;
+
printf("Reciprocal-space cell (nm^-1):\n");
basis_print(basis);
-
+
m = gsl_matrix_alloc(3, 3);
- gsl_matrix_set(m, 0, 0, basis->a.x); gsl_matrix_set(m, 0, 1, basis->b.x); gsl_matrix_set(m, 0, 2, basis->c.x);
- gsl_matrix_set(m, 1, 0, basis->a.y); gsl_matrix_set(m, 1, 1, basis->b.y); gsl_matrix_set(m, 1, 2, basis->c.y);
- gsl_matrix_set(m, 2, 0, basis->a.z); gsl_matrix_set(m, 2, 1, basis->b.z); gsl_matrix_set(m, 2, 2, basis->c.z);
-
+ gsl_matrix_set(m, 0, 0, basis->a.x);
+ gsl_matrix_set(m, 0, 1, basis->b.x);
+ gsl_matrix_set(m, 0, 2, basis->c.x);
+ gsl_matrix_set(m, 1, 0, basis->a.y);
+ gsl_matrix_set(m, 1, 1, basis->b.y);
+ gsl_matrix_set(m, 1, 2, basis->c.y);
+ gsl_matrix_set(m, 2, 0, basis->a.z);
+ gsl_matrix_set(m, 2, 1, basis->b.z);
+ gsl_matrix_set(m, 2, 2, basis->c.z);
+
perm = gsl_permutation_alloc(m->size1);
inv = gsl_matrix_alloc(m->size1, m->size2);
gsl_linalg_LU_decomp(m, perm, &s);
gsl_linalg_LU_invert(m, perm, inv);
gsl_permutation_free(perm);
gsl_matrix_free(m);
-
+
gsl_matrix_transpose(inv);
-
- ax = gsl_matrix_get(inv, 0, 0); bx = gsl_matrix_get(inv, 0, 1); cx = gsl_matrix_get(inv, 0, 2);
- ay = gsl_matrix_get(inv, 1, 0); by = gsl_matrix_get(inv, 1, 1); cy = gsl_matrix_get(inv, 1, 2);
- az = gsl_matrix_get(inv, 2, 0); bz = gsl_matrix_get(inv, 2, 1); cz = gsl_matrix_get(inv, 2, 2);
-
+
+ ax = gsl_matrix_get(inv, 0, 0);
+ bx = gsl_matrix_get(inv, 0, 1);
+ cx = gsl_matrix_get(inv, 0, 2);
+ ay = gsl_matrix_get(inv, 1, 0);
+ by = gsl_matrix_get(inv, 1, 1);
+ cy = gsl_matrix_get(inv, 1, 2);
+ az = gsl_matrix_get(inv, 2, 0);
+ bz = gsl_matrix_get(inv, 2, 1);
+ cz = gsl_matrix_get(inv, 2, 2);
+
printf("Real-space cell (nm):\n");
printf("%12.8f %12.8f %12.8f\n", ax*1e9, ay*1e9, az*1e9);
printf("%12.8f %12.8f %12.8f\n", bx*1e9, by*1e9, bz*1e9);
printf("%12.8f %12.8f %12.8f\n", cx*1e9, cy*1e9, cz*1e9);
-
+
cell.a = sqrt(ax*ax + ay*ay + az*az);
cell.b = sqrt(bx*bx + by*by + bz*bz);
cell.c = sqrt(cx*cx + cy*cy + cz*cz);
cell.alpha = acos((bx*cx + by*cy + bz*cz)/(cell.b * cell.c));
cell.beta = acos((ax*cx + ay*cy + az*cz)/(cell.a * cell.c));
cell.gamma = acos((bx*ax + by*ay + bz*az)/(cell.b * cell.a));
-
+
gsl_matrix_free(inv);
-
- printf("Cell parameters:\n");
+
+ printf("Cell parameters:\n");
cell_print(&cell);
-
+
return cell;
-
}
+
+static int basis_do_save(Basis *cell, const char *filename)
+{
+ FILE *fh;
+ UnitCell rcell;
+
+ fh = fopen(filename, "w");
+
+ fprintf(fh, "# DTR unit cell description\n");
+
+ /* A "human-readable" form */
+ rcell = basis_get_cell(cell);
+ fprintf(fh, "# a %12.8f nm\n", rcell.a*1e9);
+ fprintf(fh, "# b %12.8f nm\n", rcell.b*1e9);
+ fprintf(fh, "# c %12.8f nm\n", rcell.c*1e9);
+ fprintf(fh, "# alpha %12.8f deg\n", rad2deg(rcell.alpha));
+ fprintf(fh, "# beta %12.8f deg\n", rad2deg(rcell.beta));
+ fprintf(fh, "# gamma %12.8f deg\n", rad2deg(rcell.gamma));
+
+ /* The useful form */
+ fprintf(fh, "a %12.8f %12.8f %12.8f\n",
+ cell->a.x, cell->a.y, cell->a.z);
+ fprintf(fh, "b %12.8f %12.8f %12.8f\n",
+ cell->b.x, cell->b.y, cell->b.z);
+ fprintf(fh, "c %12.8f %12.8f %12.8f\n",
+ cell->c.x, cell->c.y, cell->c.z);
+
+ fclose(fh);
+
+ return 0;
+}
+
+static gint basis_save_response(GtkWidget *widget, gint response,
+ ControlContext *ctx)
+{
+ if ( response == GTK_RESPONSE_ACCEPT ) {
+ char *filename;
+ filename = gtk_file_chooser_get_filename(
+ GTK_FILE_CHOOSER(widget));
+ if ( basis_do_save(ctx->cell, filename) ) {
+ displaywindow_error("Failed to save unit cell.",
+ ctx->dw);
+ }
+ g_free(filename);
+ }
+
+ gtk_widget_destroy(widget);
+
+ return 0;
+}
+
+void basis_save(ControlContext *ctx)
+{
+ GtkWidget *save;
+
+ save = gtk_file_chooser_dialog_new("Save Unit Cell to File",
+ GTK_WINDOW(ctx->dw->window),
+ GTK_FILE_CHOOSER_ACTION_SAVE,
+ GTK_STOCK_CANCEL, GTK_RESPONSE_CANCEL,
+ GTK_STOCK_SAVE, GTK_RESPONSE_ACCEPT,
+ NULL);
+ g_signal_connect(G_OBJECT(save), "response",
+ G_CALLBACK(basis_save_response), ctx);
+ gtk_widget_show_all(save);
+}
+
+static int basis_do_load(Basis *cell, const char *filename)
+{
+ FILE *fh;
+ float x, y, z;
+ int got_a = 0;
+ int got_b = 0;
+ int got_c = 0;
+
+ fh = fopen(filename, "r");
+
+ while ( !feof(fh) ) {
+
+ char line[256];
+
+ if ( fgets(line, 255, fh) != NULL ) {
+
+ if ( sscanf(line, "a %f %f %f\n", &x, &y, &z) == 3 ) {
+ cell->a.x = x; cell->a.y = y; cell->a.z = z;
+ got_a = 1;
+ }
+ if ( sscanf(line, "b %f %f %f\n", &x, &y, &z) == 3 ) {
+ cell->b.x = x; cell->b.y = y; cell->b.z = z;
+ got_b = 1;
+ }
+ if ( sscanf(line, "c %f %f %f\n", &x, &y, &z) == 3 ) {
+ cell->c.x = x; cell->c.y = y; cell->c.z = z;
+ got_c = 1;
+ }
+
+ }
+
+ }
+
+ fclose(fh);
+
+ return !(got_a && got_b && got_c);
+}
+
+static gint basis_load_response(GtkWidget *widget, gint response,
+ ControlContext *ctx)
+{
+ if ( response == GTK_RESPONSE_ACCEPT ) {
+ char *filename;
+ filename = gtk_file_chooser_get_filename(
+ GTK_FILE_CHOOSER(widget));
+ if ( ctx->cell ) {
+ free(ctx->cell);
+ }
+ ctx->cell = malloc(sizeof(Basis));
+
+ if ( basis_do_load(ctx->cell, filename) ) {
+ displaywindow_error("Failed to load unit cell.",
+ ctx->dw);
+ } else {
+ displaywindow_update(ctx->dw);
+ }
+ g_free(filename);
+ }
+
+ gtk_widget_destroy(widget);
+
+ return 0;
+}
+
+void basis_load(ControlContext *ctx)
+{
+ GtkWidget *load;
+
+ load = gtk_file_chooser_dialog_new("Load Unit Cell from File",
+ GTK_WINDOW(ctx->dw->window),
+ GTK_FILE_CHOOSER_ACTION_OPEN,
+ GTK_STOCK_CANCEL, GTK_RESPONSE_CANCEL,
+ GTK_STOCK_SAVE, GTK_RESPONSE_ACCEPT,
+ NULL);
+ g_signal_connect(G_OBJECT(load), "response",
+ G_CALLBACK(basis_load_response), ctx);
+ gtk_widget_show_all(load);
+}
generated by cgit v1.2.3 (git 2.25.1) at 2024-05-17 18:36:07 +0000