diff options
| -rw-r--r-- | src/render_hkl.c | 119 |
1 files changed, 108 insertions, 11 deletions
diff --git a/src/render_hkl.c b/src/render_hkl.c index 142438e7..c08f0776 100644 --- a/src/render_hkl.c +++ b/src/render_hkl.c @@ -41,6 +41,9 @@ #include <cairo.h> #include <cairo-pdf.h> #endif +#include <gsl/gsl_matrix.h> +#include <gsl/gsl_linalg.h> +#include <gsl/gsl_blas.h> #include "utils.h" #include "symmetry.h" @@ -165,12 +168,97 @@ static void draw_circles(double xh, double xk, double xl, Reflection *refl; RefListIterator *iter; SymOpMask *m; - double nx, ny; + gsl_matrix *basis; + gsl_vector *ind; + gsl_permutation *p; + int signum; + double adx, ady, adz; + double bdx, bdy, bdz; + double cdx, cdy, cdz; + gsl_vector *za; + double asx, asy, asz; + double bsx, bsy, bsz; + double csx, csy, csz; + gsl_matrix *A; + double za_len; m = new_symopmask(sym); - nx = xh*xh + xk*xk + xl*xl; - ny = yh*yh + yk*yk + yl*yl; + /* Get the zone axis direction in cartesian coordinates */ + za = gsl_vector_alloc(3); + if ( za == NULL ) { + ERROR("Couldn't allocate za\n"); + return; + } + if ( cell_get_cartesian(cell, &adx, &ady, &adz, + &bdx, &bdy, &bdz, + &cdx, &cdy, &cdz) ) { + ERROR("Couldn't get cartesian parameters\n"); + return; + } + gsl_vector_set(za, 0, adx*zh + bdx*zk + cdx*zl); + gsl_vector_set(za, 1, ady*zh + bdy*zk + cdy*zl); + gsl_vector_set(za, 2, adz*zh + bdz*zk + cdz*zl); + + /* Normalise it (to unit length), then set its length to the + * interplanar spacing in reciprocal space, which is 1/the length of + * the direct space ZA vector */ + za_len = gsl_blas_dnrm2(za); + gsl_blas_dscal(1.0/za_len, za); + gsl_blas_dscal(1.0/za_len, za); + + /* Express it in terms of the basis vectors of the reciprocal lattice */ + if ( cell_get_reciprocal(cell, &asx, &asy, &asz, + &bsx, &bsy, &bsz, + &csx, &csy, &csz) ) { + ERROR("Couldn't get reciprocal parameters\n"); + return; + } + + A = gsl_matrix_alloc(3, 3); + if ( A == NULL ) { + ERROR("Couldn't allocate A\n"); + return; + } + gsl_matrix_set(A, 0, 0, asx); + gsl_matrix_set(A, 1, 0, asy); + gsl_matrix_set(A, 2, 0, asz); + gsl_matrix_set(A, 0, 1, bsx); + gsl_matrix_set(A, 1, 1, bsy); + gsl_matrix_set(A, 2, 1, bsz); + gsl_matrix_set(A, 0, 2, csx); + gsl_matrix_set(A, 1, 2, csy); + gsl_matrix_set(A, 2, 2, csz); + + p = gsl_permutation_alloc(3); + + gsl_linalg_LU_decomp(A, p, &signum); + gsl_linalg_LU_svx(A, p, za); + STATUS("Zone axis along %5.2e %5.2e %5.2e in the reciprocal lattice\n", + gsl_vector_get(za, 0), + gsl_vector_get(za, 1), + gsl_vector_get(za, 2)); + + gsl_matrix_free(A); + + basis = gsl_matrix_alloc(3, 3); + if ( basis == NULL ) return; + + gsl_matrix_set(basis, 0, 0, xh); + gsl_matrix_set(basis, 1, 0, xk); + gsl_matrix_set(basis, 2, 0, xl); + gsl_matrix_set(basis, 0, 1, yh); + gsl_matrix_set(basis, 1, 1, yk); + gsl_matrix_set(basis, 2, 1, yl); + gsl_matrix_set(basis, 0, 2, gsl_vector_get(za, 0)); + gsl_matrix_set(basis, 1, 2, gsl_vector_get(za, 1)); + gsl_matrix_set(basis, 2, 2, gsl_vector_get(za, 2)); + gsl_linalg_LU_decomp(basis, p, &signum); + + gsl_vector_free(za); + + ind = gsl_vector_alloc(3); + if ( ind == NULL ) return; /* Iterate over all reflections */ for ( refl = first_refl(list, &iter); @@ -195,10 +283,14 @@ static void draw_circles(double xh, double xk, double xl, get_equiv(sym, m, i, ha, ka, la, &h, &k, &l); /* Is the reflection in the zone? */ - if ( h*zh + k*zk + l*zl != zone ) continue; + if ( h*zh + k*zk + l*zl != zone) continue; - xi = (h*xh + k*xk + l*xl) / nx; - yi = (h*yh + k*yk + l*yl) / ny; + gsl_vector_set(ind, 0, h); + gsl_vector_set(ind, 1, k); + gsl_vector_set(ind, 2, l); + gsl_linalg_LU_svx(basis, p, ind); + xi = gsl_vector_get(ind, 0); + yi = gsl_vector_get(ind, 1); switch ( wght) { @@ -243,6 +335,10 @@ static void draw_circles(double xh, double xk, double xl, } + gsl_matrix_free(basis); + gsl_vector_free(ind); + gsl_permutation_free(p); + free_symopmask(m); } @@ -326,10 +422,10 @@ static void render_za(UnitCell *cell, RefList *list, double rmin, rmax; /* Vector product to determine the zone axis. */ - zh = xk*yl - xl*yk; - zk = - xh*yl + xl*yh; - zl = xh*yk - xk*yh; - STATUS("Zone axis is %i %i %i\n", zh, zk, zl); + zh = yk*xl - yl*xk; + zk = - yh*xl + yl*xh; + zl = yh*xk - yk*xh; + STATUS("Zone axis is [%i %i %i]\n", zh, zk, zl); /* Size of output and centre definition */ wh = 1024; @@ -428,7 +524,7 @@ static void render_za(UnitCell *cell, RefList *list, /* Draw indexing lines */ cairo_set_line_cap(dctx, CAIRO_LINE_CAP_ROUND); - cairo_set_line_width(dctx, 4.0); + cairo_set_line_width(dctx, 2.0); cairo_move_to(dctx, (double)cx, (double)cy); u = rmax*sin(theta); v = rmax*cos(theta); @@ -447,6 +543,7 @@ static void render_za(UnitCell *cell, RefList *list, snprintf(tmp, 255, "%i%i%i", abs(yh), abs(yk), abs(yl)); cairo_text_extents(dctx, tmp, &size); + cairo_set_line_width(dctx, 2.0); cairo_move_to(dctx, (double)cx, (double)cy); cairo_line_to(dctx, cx, cy+rmax*scale); cairo_set_source_rgb(dctx, 0.0, 1.0, 0.0); |