diff options
author | taw27 <taw27@bf6ca9ba-c028-0410-8290-897cf20841d1> | 2007-02-05 21:12:57 +0000 |
---|---|---|
committer | taw27 <taw27@bf6ca9ba-c028-0410-8290-897cf20841d1> | 2007-02-05 21:12:57 +0000 |
commit | 05b5d261682b9136fb46476a64eab6980b0dba64 (patch) | |
tree | d7faa450b69cf2104ffff7fc89a95914284d23af /src/reflections.c |
Initial import
git-svn-id: svn://cook.msm.cam.ac.uk:745/diff-tomo/dtr@1 bf6ca9ba-c028-0410-8290-897cf20841d1
Diffstat (limited to 'src/reflections.c')
-rw-r--r-- | src/reflections.c | 186 |
1 files changed, 186 insertions, 0 deletions
diff --git a/src/reflections.c b/src/reflections.c new file mode 100644 index 0000000..dbe2a7d --- /dev/null +++ b/src/reflections.c @@ -0,0 +1,186 @@ +/* + * reflections.c + * + * Data structures in 3D space + * + * (c) 2007 Thomas White <taw27@cam.ac.uk> + * dtr - Diffraction Tomography Reconstruction + * + */ + +#ifdef HAVE_CONFIG_H +#include <config.h> +#endif + +#include <stdlib.h> +#include <math.h> +#include <stdio.h> + +#include "reflections.h" + +static void reflection_addfirst(ReflectionContext *reflectionctx) { + /* Create first items on lists - saves faffing later. Corresponds to a central marker. + Reflections are only stored if they have non-zero value. */ + reflectionctx->reflections = malloc(sizeof(Reflection)); + reflectionctx->reflections->next = NULL; + reflectionctx->reflections->x = 0; + reflectionctx->reflections->y = 0; + reflectionctx->reflections->z = 0; + reflectionctx->reflections->type = REFLECTION_CENTRAL; + reflectionctx->last_reflection = reflectionctx->reflections; +} + +/* Size of the volume (arbitrary units - m-1 in DTR), number of reflections across. The centre is defined to be at nx/2, ny/2, nz/2. + This defines the three-dimensional grid of reflections */ +ReflectionContext *reflection_init() { + + ReflectionContext *reflectionctx = malloc(sizeof(ReflectionContext)); + reflection_addfirst(reflectionctx); + + return reflectionctx; + +} + +void reflection_clear(ReflectionContext *reflectionctx) { + + Reflection *reflection = reflectionctx->reflections; + do { + Reflection *next = reflection->next; + free(reflection); + reflection = next; + } while ( reflection ); + + reflection_addfirst(reflectionctx); + +} + +void reflection_add(ReflectionContext *reflectionctx, double x, double y, double z, double intensity, ReflectionType type) { + + Reflection *new_reflection; + + new_reflection = malloc(sizeof(Reflection)); + new_reflection->next = NULL; + new_reflection->x = x; + new_reflection->y = y; + new_reflection->z = z; + new_reflection->intensity = intensity; + new_reflection->type = type; + + reflectionctx->last_reflection->next = new_reflection; + reflectionctx->last_reflection = new_reflection; + +} + +void reflection_add_index(ReflectionContext *reflectionctx, signed int h, signed int k, signed int l, double intensity, ReflectionType type) { + + Reflection *new_reflection; + + new_reflection = malloc(sizeof(Reflection)); + new_reflection->next = NULL; + new_reflection->h = h; + new_reflection->k = k; + new_reflection->l = l; + new_reflection->intensity = intensity; + new_reflection->type = type; + + reflectionctx->last_reflection->next = new_reflection; + reflectionctx->last_reflection = new_reflection; + +} + +/* x and y in metres (in real space!) */ +void reflection_add_from_dp(ControlContext *ctx, double x, double y, double tilt_degrees, double intensity) { + + double nx, ny, nz; + + /* Real space */ + double L; + double d; + + /* Shared */ + double theta; + double psi; + + /* Reciprocal space */ + double r; + double tilt; + double omega; + double x_temp, y_temp, z_temp; + + tilt = 2*M_PI*(tilt_degrees/360); /* Convert to Radians */ + omega = 2*M_PI*(ctx->omega/360); /* Likewise */ + + d = sqrt((x*x) + (y*y)); + L = ctx->camera_length; + theta = atan2(d, L); + psi = atan2(y, x); + + r = 1/ctx->lambda; + x_temp = r*sin(theta)*cos(psi); + y_temp = -r*sin(theta)*sin(psi); /* Minus sign to define axes as y going upwards */ + z_temp = r- r*cos(theta); + + /* Apply the rotations... + First: rotate image clockwise until tilt axis is aligned horizontally. */ + nx = x_temp*cos(omega) + y_temp*-sin(omega); + ny = x_temp*sin(omega) + y_temp*cos(omega); + nz = z_temp; + + /* Now, tilt about the x-axis ANTICLOCKWISE around +x, i.e. the "wrong" way. + This is because the crystal is rotated in the experiment, not the Ewald sphere. */ + x_temp = nx; y_temp = ny; z_temp = nz; + nx = x_temp; + ny = cos(tilt)*y_temp - sin(tilt)*z_temp; + nz = -sin(tilt)*y_temp - cos(tilt)*z_temp; + + reflection_add(ctx->reflectionctx, nx, ny, nz, intensity, REFLECTION_NORMAL); + +} + +/* Alternative formulation for when the input data is already in reciprocal space units + x and y in metres^-1 (in reciprocal space */ +void reflection_add_from_reciprocal(ControlContext *ctx, double x, double y, double tilt_degrees, double intensity) { + + double nx, ny, nz; + + /* Input (reciprocal space) */ + double dr; + + /* Shared */ + double theta; + double psi; + double r; + + /* Reciprocal space */ + double tilt; + double omega; + double x_temp, y_temp, z_temp; + + tilt = M_PI*(tilt_degrees/180); /* Convert to Radians */ + omega = M_PI*(ctx->omega/180); /* Likewise */ + + dr = sqrt((x*x) + (y*y)); + r = 1/ctx->lambda; + theta = atan2(dr, r); + psi = atan2(y, x); + + x_temp = r*sin(theta)*cos(psi); + y_temp = -r*sin(theta)*sin(psi); /* Minus sign to define axes as y going upwards */ + z_temp = r - r*cos(theta); + + /* Apply the rotations... + First: rotate image clockwise until tilt axis is aligned horizontally. */ + nx = x_temp*cos(omega) + y_temp*-sin(omega); + ny = x_temp*sin(omega) + y_temp*cos(omega); + nz = z_temp; + + /* Now, tilt about the x-axis ANTICLOCKWISE around +x, i.e. the "wrong" way. + This is because the crystal is rotated in the experiment, not the Ewald sphere. */ + x_temp = nx; y_temp = ny; z_temp = nz; + nx = x_temp; + ny = cos(tilt)*y_temp - sin(tilt)*z_temp; + nz = -sin(tilt)*y_temp - cos(tilt)*z_temp; + + reflection_add(ctx->reflectionctx, nx, ny, nz, intensity, REFLECTION_NORMAL); + +} |