From 4e1e2ef4472e28a146e6c83d053f1fc4d2419f88 Mon Sep 17 00:00:00 2001 From: Thomas White Date: Wed, 14 Oct 2009 17:35:18 +0200 Subject: Calculate reflections using reciprocal cell --- src/Makefile.am | 4 +- src/cell.c | 67 +++++++++++++++++++++++++++++---- src/cell.h | 13 +++++-- src/image.c | 11 +----- src/image.h | 8 ++-- src/relrod.c | 114 ++++++++++++++++++++++++++++---------------------------- src/sim-main.c | 18 ++++----- 7 files changed, 140 insertions(+), 95 deletions(-) diff --git a/src/Makefile.am b/src/Makefile.am index b480cedd..46a5e9d9 100644 --- a/src/Makefile.am +++ b/src/Makefile.am @@ -1,8 +1,8 @@ bin_PROGRAMS = template_index simulate_patterns template_index_SOURCES = main.c relrod.c utils.c image.c cell.c -template_index_LDADD = @LIBS@ -lm +template_index_LDADD = @LIBS@ -lm -lgsl -lgslcblas AM_CFLAGS = -Wall -g simulate_patterns_SOURCES = sim-main.c relrod.c utils.c image.c cell.c -simulate_patterns_LDADD = @LIBS@ -lm +simulate_patterns_LDADD = @LIBS@ -lm -lgsl -lgslcblas diff --git a/src/cell.c b/src/cell.c index 6fa91cd3..81638016 100644 --- a/src/cell.c +++ b/src/cell.c @@ -16,6 +16,9 @@ #include #include #include +#include +#include +#include #include "cell.h" #include "utils.h" @@ -74,9 +77,9 @@ static void cell_update_crystallographic(UnitCell *cell) cell->gamma = angle_between(cell->ax, cell->ay, cell->az, cell->bx, cell->by, cell->bz); - printf("a=%f nm\n", cell->a); - printf("b=%f nm\n", cell->b); - printf("c=%f nm\n", cell->c); + printf("a=%f nm\n", cell->a/1e9); + printf("b=%f nm\n", cell->b/1e9); + printf("c=%f nm\n", cell->c/1e9); printf("alpha = %f deg\n", rad2deg(cell->alpha)); printf(" beta = %f deg\n", rad2deg(cell->beta)); printf("gamma = %f deg\n", rad2deg(cell->gamma)); @@ -104,7 +107,7 @@ UnitCell *cell_new() void cell_set_parameters(UnitCell *cell, double a, double b, double c, - double alpha, double beta, double gamma) + double alpha, double beta, double gamma) { if ( !cell ) return; @@ -135,7 +138,7 @@ void cell_set_cartesian(UnitCell *cell, UnitCell *cell_new_from_parameters(double a, double b, double c, - double alpha, double beta, double gamma) + double alpha, double beta, double gamma) { UnitCell *cell; @@ -149,9 +152,9 @@ UnitCell *cell_new_from_parameters(double a, double b, double c, void cell_get_cartesian(UnitCell *cell, - double *ax, double *ay, double *az, - double *bx, double *by, double *bz, - double *cx, double *cy, double *cz) + double *ax, double *ay, double *az, + double *bx, double *by, double *bz, + double *cx, double *cy, double *cz) { if ( !cell ) return; @@ -159,3 +162,51 @@ void cell_get_cartesian(UnitCell *cell, *bx = cell->bx; *by = cell->by; *bz = cell->bz; *cx = cell->cx; *cy = cell->cy; *cz = cell->cz; } + + +void cell_get_reciprocal(UnitCell *cell, + double *asx, double *asy, double *asz, + double *bsx, double *bsy, double *bsz, + double *csx, double *csy, double *csz) +{ + int s; + gsl_matrix *m; + gsl_matrix *inv; + gsl_permutation *perm; + + m = gsl_matrix_alloc(3, 3); + gsl_matrix_set(m, 0, 0, cell->ax); + gsl_matrix_set(m, 0, 1, cell->bx); + gsl_matrix_set(m, 0, 2, cell->cx); + gsl_matrix_set(m, 1, 0, cell->ay); + gsl_matrix_set(m, 1, 1, cell->by); + gsl_matrix_set(m, 1, 2, cell->cy); + gsl_matrix_set(m, 2, 0, cell->az); + gsl_matrix_set(m, 2, 1, cell->bz); + gsl_matrix_set(m, 2, 2, cell->cz); + + /* Invert */ + 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); + + /* Transpose */ + gsl_matrix_transpose(inv); + + *asx = gsl_matrix_get(inv, 0, 0); + *bsx = gsl_matrix_get(inv, 0, 1); + *csx = gsl_matrix_get(inv, 0, 2); + *asy = gsl_matrix_get(inv, 1, 0); + *bsy = gsl_matrix_get(inv, 1, 1); + *csy = gsl_matrix_get(inv, 1, 2); + *asz = gsl_matrix_get(inv, 2, 0); + *bsz = gsl_matrix_get(inv, 2, 1); + *csz = gsl_matrix_get(inv, 2, 2); + + printf("a* = %f %f %f\n", *asx/1e9, *asy/1e9, *asz/1e9); + printf("b* = %f %f %f\n", *bsx/1e9, *bsy/1e9, *bsz/1e9); + printf("c* = %f %f %f\n", *csx/1e9, *csy/1e9, *csz/1e9); +} diff --git a/src/cell.h b/src/cell.h index cd31aa66..196b9521 100644 --- a/src/cell.h +++ b/src/cell.h @@ -19,9 +19,9 @@ typedef struct { /* Crystallographic representation */ - double a; /* nm */ - double b; /* nm */ - double c; /* nm */ + double a; /* m */ + double b; /* m */ + double c; /* m */ double alpha; /* Radians */ double beta; /* Radians */ double gamma; /* Radians */ @@ -40,7 +40,7 @@ typedef struct { extern UnitCell *cell_new(void); -/* Lengths in nm, angles in radians */ +/* Lengths in m, angles in radians */ extern UnitCell *cell_new_from_parameters(double a, double b, double c, double alpha, double beta, double gamma); @@ -57,4 +57,9 @@ extern void cell_get_cartesian(UnitCell *cell, double *bx, double *by, double *bz, double *cx, double *cy, double *cz); +extern void cell_get_reciprocal(UnitCell *cell, + double *asx, double *asy, double *asz, + double *bsx, double *bsy, double *bsz, + double *csx, double *csy, double *csz); + #endif /* CELL_H */ diff --git a/src/image.c b/src/image.c index 8780c0dd..d3e65685 100644 --- a/src/image.c +++ b/src/image.c @@ -69,8 +69,8 @@ ImageList *image_list_new() } -void image_add_feature_reflection(ImageFeatureList *flist, double x, double y, - struct image *parent, double intensity) +void image_add_feature(ImageFeatureList *flist, double x, double y, + struct image *parent, double intensity) { if ( flist->features ) { flist->features = realloc(flist->features, @@ -93,13 +93,6 @@ void image_add_feature_reflection(ImageFeatureList *flist, double x, double y, } -void image_add_feature(ImageFeatureList *flist, double x, double y, - struct image *parent, double intensity) -{ - image_add_feature_reflection(flist, x, y, parent, intensity); -} - - ImageFeatureList *image_feature_list_new() { ImageFeatureList *flist; diff --git a/src/image.h b/src/image.h index 69b8cc79..fc8a1f2f 100644 --- a/src/image.h +++ b/src/image.h @@ -82,10 +82,13 @@ struct image { typedef struct _imagelist ImageList; +/* Image lists */ extern ImageList *image_list_new(void); extern int image_add(ImageList *list, struct image *image); + +/* Feature lists */ extern ImageFeatureList *image_feature_list_new(void); extern void image_feature_list_free(ImageFeatureList *flist); @@ -93,11 +96,6 @@ extern void image_feature_list_free(ImageFeatureList *flist); extern void image_add_feature(ImageFeatureList *flist, double x, double y, struct image *parent, double intensity); -extern void image_add_feature_reflection(ImageFeatureList *flist, - double x, double y, - struct image *parent, - double intensity); - extern struct imagefeature *image_feature_closest(ImageFeatureList *flist, double x, double y, double *d, int *idx); diff --git a/src/relrod.c b/src/relrod.c index df4c4df3..e344133a 100644 --- a/src/relrod.c +++ b/src/relrod.c @@ -62,14 +62,18 @@ void get_reflections(struct image *image, UnitCell *cell) double kx, ky, kz; /* Electron wavevector ("normal" times 1/lambda) */ double ux, uy, uz; /* "up" vector */ double rx, ry, rz; /* "right" vector */ - double ax, ay, az; /* "a" lattice parameter */ - double bx, by, bz; /* "b" lattice parameter */ - double cx, cy, cz; /* "c" lattice parameter */ + double asx, asy, asz; /* "a*" lattice parameter */ + double bsx, bsy, bsz; /* "b*" lattice parameter */ + double csx, csy, csz; /* "c*" lattice parameter */ signed int h, k, l; - flist = image_feature_list_new(); - cell_get_cartesian(cell, &ax, &ay, &az, &bx, &by, &bz, &cx, &cy, &cz); + /* Get the reciprocal unit cell */ + cell_get_reciprocal(cell, &asx, &asy, &asz, + &bsx, &bsy, &bsz, + &csx, &csy, &csz); + /* Prepare list and some parameters */ + flist = image_feature_list_new(); tilt = image->tilt; omega = image->omega; wavenumber = 1.0/image->lambda; @@ -96,9 +100,9 @@ void get_reflections(struct image *image, UnitCell *cell) double g_sq, gn; /* Get the coordinates of the reciprocal lattice point */ - xl = h*ax + k*bx + l*cx; - yl = h*ay + k*by + l*cy; - zl = h*az + k*bz + l*cz; + xl = h*asx + k*bsx + l*csx; + yl = h*asy + k*bsy + l*csy; + zl = h*asz + k*bsz + l*csz; g_sq = modulus_squared(xl, yl, zl); gn = xl*nx + yl*ny + zl*nz; @@ -119,6 +123,7 @@ void get_reflections(struct image *image, UnitCell *cell) double gx, gy, gz; double theta; double x, y; + double dx, dy, psi; /* Determine the intersection point */ xi = xl + s*nx; yi = yl + s*ny; zi = zl + s*nz; @@ -130,56 +135,49 @@ void get_reflections(struct image *image, UnitCell *cell) * the sphere to the intersection */ theta = angle_between(-kx, -ky, -kz, gx, gy, gz); - if ( theta > 0.0 ) { - - double dx, dy, psi; - - /* Calculate azimuth of point in image - * (anticlockwise from +x) */ - dx = xi*rx + yi*ry + zi*rz; - dy = xi*ux + yi*uy + zi*uz; - psi = atan2(dy, dx); - - /* Get image coordinates from polar - * representation */ - if ( image->fmode == FORMULATION_CLEN ) { - x = image->camera_len*tan(theta)*cos(psi); - y = image->camera_len*tan(theta)*sin(psi); - x *= image->resolution; - y *= image->resolution; - } else if ( image->fmode==FORMULATION_PIXELSIZE ) { - x = tan(theta)*cos(psi) / image->lambda; - y = tan(theta)*sin(psi) / image->lambda; - x /= image->pixel_size; - y /= image->pixel_size; - } else { - fprintf(stderr, - "Unrecognised formulation mode " - "in get_reflections\n"); - return; - } - - x += image->x_centre; - y += image->y_centre; - - /* Sanity check */ - if ( (x>=0) && (xwidth) - && (y>=0) && (yheight) ) { - - /* Record the reflection - * Intensity should be multiplied by - * relrod spike function except - * reprojected reflections aren't used - * quantitatively for anything - */ - image_add_feature_reflection(flist, x, - y, image, 1.0); - - } /* else it's outside the picture somewhere */ - - } /* else this is the central beam */ - - } + /* Calculate azimuth of point in image + * (anticlockwise from +x) */ + dx = xi*rx + yi*ry + zi*rz; + dy = xi*ux + yi*uy + zi*uz; + psi = atan2(dy, dx); + + /* Get image coordinates from polar + * representation */ + if ( image->fmode == FORMULATION_CLEN ) { + x = image->camera_len*tan(theta)*cos(psi); + y = image->camera_len*tan(theta)*sin(psi); + x *= image->resolution; + y *= image->resolution; + } else if ( image->fmode==FORMULATION_PIXELSIZE ) { + x = tan(theta)*cos(psi) / image->lambda; + y = tan(theta)*sin(psi) / image->lambda; + x /= image->pixel_size; + y /= image->pixel_size; + } else { + fprintf(stderr, + "Unrecognised formulation mode " + "in get_reflections\n"); + return; + } + + x += image->x_centre; + y += image->y_centre; + + /* Sanity check */ + if ( (x>=0) && (xwidth) + && (y>=0) && (yheight) ) { + + /* Record the reflection. + * Intensity should be multiplied by relrod + * spike function, except reprojected + * reflections aren't used quantitatively for + * anything. */ + printf("%i %i %i\n", h,k,l); + image_add_feature(flist, x, y, image, 1.0); + + } /* else it's outside the picture somewhere */ + + } /* else reflection is not excited in this orientation */ } } diff --git a/src/sim-main.c b/src/sim-main.c index 38ad45bc..23ba3e99 100644 --- a/src/sim-main.c +++ b/src/sim-main.c @@ -61,20 +61,20 @@ int main(int argc, char *argv[]) image.tilt = 0.0; image.omega = 0.0; image.fmode = FORMULATION_CLEN; - image.x_centre = 128; - image.y_centre = 128; - image.camera_len = 1.0; /* one metre */ - image.resolution = 5120; - image.lambda = 0.6e-9; + image.x_centre = 255.5; + image.y_centre = 255.5; + image.camera_len = 0.2; /* 20 cm */ + image.resolution = 5120; /* 512 pixels in 10 cm */ + image.lambda = 0.2e-9; /* LCLS wavelength */ image.data = malloc(512*512*2); image_add(list, &image); - cell = cell_new_from_parameters(1.0, - 1.0, - 1.0, + cell = cell_new_from_parameters(28.10e-9, + 28.10e-9, + 16.52e-9, deg2rad(90.0), deg2rad(90.0), - deg2rad(90.0)); + deg2rad(120.0)); get_reflections(&image, cell); -- cgit v1.2.3