Image, feature and unit cell infrastructure

Brought across from DTR and Synth3D

1 files changed, 180 insertions, 0 deletions

diff --git a/src/relrod.c b/src/relrod.c
new file mode 100644
index 00000000..674f7971
--- /dev/null
+++ b/src/relrod.c
@@ -0,0 +1,180 @@
+/*
+ * relrod.c
+ *
+ * Calculate reflection positions via line-sphere intersection test
+ *
+ * (c) 2007-2009 Thomas White <thomas.white@desy.de>
+ *
+ * template_index - Indexing diffraction patterns by template matching
+ *
+ */
+
+
+#include <stdlib.h>
+#include <math.h>
+#include <stdio.h>
+
+#include "image.h"
+#include "utils.h"
+#include "cell.h"
+
+
+static void mapping_rotate(double x, double y, double z,
+                           double *ddx, double *ddy, double *ddz,
+                           double omega, double tilt)
+{
+	double nx, ny, nz;
+	double x_temp, y_temp, z_temp;
+
+	/* First: rotate image clockwise until tilt axis is aligned
+	 * horizontally. */
+	nx = x*cos(omega) + y*sin(omega);
+	ny = -x*sin(omega) + y*cos(omega);
+	nz = z;
+
+	/* 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;
+
+	/* Finally, reverse the omega rotation to restore the location of the
+	 * image in 3D space */
+	x_temp = nx; y_temp = ny; z_temp = nz;
+	nx = x_temp*cos(-omega) + y_temp*sin(-omega);
+	ny = -x_temp*sin(-omega) + y_temp*cos(-omega);
+	nz = z_temp;
+
+	*ddx = nx;
+	*ddy = ny;
+	*ddz = nz;
+}
+
+
+void get_reflections(struct image *image, UnitCell *cell)
+{
+	ImageFeatureList *flist;
+	double smax = 0.2e9;
+	double tilt, omega, k;
+	double nx, ny, nz; /* "normal" vector */
+	double kx, ky, kz; /* Electron wavevector ("normal" times 1/lambda */
+	double ux, uy, uz; /* "up" vector */
+	double rx, ry, rz; /* "right" vector */
+
+	flist = image_feature_list_new();
+
+	tilt = image->tilt;
+	omega = image->omega;
+	k = 1.0/image->lambda;
+
+	/* Calculate the (normalised) incident electron wavevector */
+	mapping_rotate(0.0, 0.0, 1.0, &nx, &ny, &nz, omega, tilt);
+	kx = nx / image->lambda;
+	ky = ny / image->lambda;
+	kz = nz / image->lambda;	/* This is the centre of the Ewald sphere */
+
+	/* Determine where "up" is */
+	mapping_rotate(0.0, 1.0, 0.0, &ux, &uy, &uz, omega, tilt);
+
+	/* Determine where "right" is */
+	mapping_rotate(1.0, 0.0, 0.0, &rx, &ry, &rz, omega, tilt);
+
+	do {
+
+		double xl, yl, zl;
+		double a, b, c;
+		double s1, s2, s, t;
+		double g_sq, gn;
+
+		/* Get the coordinates of the reciprocal lattice point */
+	//	xl = reflection->x;
+	//	yl = reflection->y;
+	//	zl = reflection->z;
+		g_sq = modulus_squared(xl, yl, zl);
+		gn = xl*nx + yl*ny + zl*nz;
+
+		/* Next, solve the relrod equation to calculate
+		 * the excitation error */
+		a = 1.0;
+		b = 2.0*(gn - k);
+		c = -2.0*gn*k + g_sq;
+		t = -0.5*(b + sign(b)*sqrt(b*b - 4.0*a*c));
+		s1 = t/a;
+		s2 = c/t;
+		if ( fabs(s1) < fabs(s2) ) s = s1; else s = s2;
+
+		/* Skip this reflection if s is large */
+		if ( fabs(s) <= smax ) {
+
+			double xi, yi, zi;
+			double gx, gy, gz;
+			double theta;
+			double x, y;
+
+			/* Determine the intersection point */
+			xi = xl + s*nx;  yi = yl + s*ny;  zi = zl + s*nz;
+
+			/* Calculate Bragg angle */
+			gx = xi - kx;
+			gy = yi - ky;
+			gz = zi - kz;	/* This is the vector from the centre of
+					*  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 reproject_get_reflections\n");
+					return NULL;
+				}
+
+				x += image->x_centre;
+				y += image->y_centre;
+
+				/* Sanity check */
+				if ( (x>=0) && (x<image->width)
+				  && (y>=0) && (y<image->height) ) {
+
+					/* 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 */
+
+		}
+
+	} while ( 1 );
+
+	image->rflist = flist;
+}