Initial import

git-svn-id: svn://cook.msm.cam.ac.uk:745/diff-tomo/dtr@1 bf6ca9ba-c028-0410-8290-897cf20841d1

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);
+
+}