aboutsummaryrefslogtreecommitdiff
diff options
context:
space:
mode:
Diffstat
-rw-r--r--.gitignore1
-rw-r--r--src/Makefile.am8
-rw-r--r--src/cell.c148
-rw-r--r--src/cell.h55
-rw-r--r--src/image.c174
-rw-r--r--src/image.h121
-rw-r--r--src/main.c2
-rw-r--r--src/relrod.c180
-rw-r--r--src/relrod.h24
-rw-r--r--src/sim-main.c82
-rw-r--r--src/utils.c124
-rw-r--r--src/utils.h42
12 files changed, 959 insertions, 2 deletions
diff --git a/.gitignore b/.gitignore
index 3c13d553..692c0da4 100644
--- a/.gitignore
+++ b/.gitignore
@@ -10,3 +10,4 @@ stamp-h1
src/*.o
src/.deps
src/template_index
+src/simulate_patterns
diff --git a/src/Makefile.am b/src/Makefile.am
index ac173bb1..b480cedd 100644
--- a/src/Makefile.am
+++ b/src/Makefile.am
@@ -1,4 +1,8 @@
-bin_PROGRAMS = template_index
-template_index_SOURCES = main.c
+bin_PROGRAMS = template_index simulate_patterns
+
+template_index_SOURCES = main.c relrod.c utils.c image.c cell.c
template_index_LDADD = @LIBS@ -lm
AM_CFLAGS = -Wall -g
+
+simulate_patterns_SOURCES = sim-main.c relrod.c utils.c image.c cell.c
+simulate_patterns_LDADD = @LIBS@ -lm
diff --git a/src/cell.c b/src/cell.c
new file mode 100644
index 00000000..d2495b0f
--- /dev/null
+++ b/src/cell.c
@@ -0,0 +1,148 @@
+/*
+ * cell.c
+ *
+ * Unit Cell Calculations
+ *
+ * (c) 2007-2009 Thomas White <thomas.white@desy.de>
+ *
+ * template_index - Indexing diffraction patterns by template matching
+ *
+ */
+
+#ifdef HAVE_CONFIG_H
+#include <config.h>
+#endif
+
+#include <math.h>
+#include <stdlib.h>
+#include <stdio.h>
+
+#include "cell.h"
+#include "utils.h"
+
+
+/* Update the cartesian representation from the crystallographic one */
+static void cell_update_cartesian(UnitCell *cell)
+{
+ double tmp, V, cosalphastar, cstar;
+
+ if ( !cell ) return;
+
+ /* a in terms of x, y and z
+ * +a (cryst) is defined to lie along +x (cart) */
+ cell->ax = cell->a;
+ cell->ay = 0.0;
+ cell->az = 0.0;
+
+ /* b in terms of x, y and z
+ * b (cryst) is defined to lie in the xy (cart) plane */
+ cell->bx = cell->b*cos(cell->gamma);
+ cell->by = cell->b*sin(cell->gamma);
+ cell->bz = 0.0;
+
+ tmp = cos(cell->alpha)*cos(cell->alpha)
+ - cos(cell->beta)*cos(cell->beta)
+ - cos(cell->gamma)*cos(cell->gamma)
+ + 2.0*cos(cell->alpha)*cos(cell->beta)*cos(cell->gamma);
+ V = cell->a * cell->b * cell->c * sqrt(1.0 - tmp);
+
+ cosalphastar = cos(cell->beta)*cos(cell->gamma) - cos(cell->alpha);
+ cosalphastar /= sin(cell->beta)*sin(cell->gamma);
+
+ cstar = (cell->a * cell->b * sin(cell->gamma))/V;
+
+ /* c in terms of x, y and z */
+ cell->cx = cell->c*cos(cell->beta);
+ cell->cy = -cell->c*sin(cell->beta)*cosalphastar;
+ cell->cz = 1.0/cstar;
+}
+
+
+/* Update the crystallographic representation from the cartesian one */
+static void cell_update_crystallographic(UnitCell *cell)
+{
+ if ( !cell ) return;
+
+ cell->a = modulus(cell->ax, cell->ay, cell->az);
+ cell->b = modulus(cell->bx, cell->by, cell->bz);
+ cell->c = modulus(cell->cx, cell->cy, cell->cz);
+
+ cell->alpha = angle_between(cell->bx, cell->by, cell->bz,
+ cell->cx, cell->cy, cell->cz);
+ cell->beta = angle_between(cell->ax, cell->ay, cell->az,
+ cell->cx, cell->cy, cell->cz);
+ 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("alpha = %f deg\n", rad2deg(cell->alpha));
+ printf(" beta = %f deg\n", rad2deg(cell->beta));
+ printf("gamma = %f deg\n", rad2deg(cell->gamma));
+}
+
+
+UnitCell *cell_new()
+{
+ UnitCell *cell;
+
+ cell = malloc(sizeof(UnitCell));
+ if ( !cell ) return NULL;
+
+ cell->a = 1.0;
+ cell->b = 1.0;
+ cell->c = 1.0;
+ cell->alpha = M_PI_2;
+ cell->beta = M_PI_2;
+ cell->gamma = M_PI_2;
+
+ cell_update_cartesian(cell);
+
+ return cell;
+}
+
+
+void cell_set_parameters(UnitCell *cell, double a, double b, double c,
+ double alpha, double beta, double gamma)
+{
+ if ( !cell ) return;
+
+ cell->a = a;
+ cell->b = b;
+ cell->c = c;
+ cell->alpha = alpha;
+ cell->beta = beta;
+ cell->gamma = gamma;
+
+ cell_update_cartesian(cell);
+}
+
+
+void cell_set_cartesian(UnitCell *cell,
+ double ax, double ay, double az,
+ double bx, double by, double bz,
+ double cx, double cy, double cz)
+{
+ if ( !cell ) return;
+
+ cell->ax = ax; cell->ay = ay; cell->az = az;
+ cell->bx = bx; cell->by = by; cell->bz = bz;
+ cell->cx = cx; cell->cy = cy; cell->cz = cz;
+
+ cell_update_crystallographic(cell);
+}
+
+
+UnitCell *cell_new_from_parameters(double a, double b, double c,
+ double alpha, double beta, double gamma)
+{
+ UnitCell *cell;
+
+ cell = cell_new();
+ if ( !cell ) return NULL;
+
+ cell_set_parameters(cell, a, b, c, alpha, beta, gamma);
+
+ return cell;
+}
diff --git a/src/cell.h b/src/cell.h
new file mode 100644
index 00000000..e99c88df
--- /dev/null
+++ b/src/cell.h
@@ -0,0 +1,55 @@
+/*
+ * cell.h
+ *
+ * Unit Cell Calculations
+ *
+ * (c) 2007-2009 Thomas White <thomas.white@desy.de>
+ *
+ * template_index - Indexing diffraction patterns by template matching
+ *
+ */
+
+#ifndef CELL_H
+#define CELL_H
+
+#ifdef HAVE_CONFIG_H
+#include <config.h>
+#endif
+
+typedef struct {
+
+ /* Crystallographic representation */
+ double a; /* nm */
+ double b; /* nm */
+ double c; /* nm */
+ double alpha; /* Radians */
+ double beta; /* Radians */
+ double gamma; /* Radians */
+
+ /* Cartesian representation */
+ double ax; double bx; double cx;
+ double ay; double by; double cy;
+ double az; double bz; double cz;
+
+ /* Cartesian representation of reciprocal axes */
+ double axs; double bxs; double cxs;
+ double ays; double bys; double cys;
+ double azs; double bzs; double czs;
+
+} UnitCell;
+
+extern UnitCell *cell_new(void);
+
+/* Lengths in nm, angles in radians */
+extern UnitCell *cell_new_from_parameters(double a, double b, double c,
+ double alpha, double beta, double gamma);
+
+extern void cell_set_cartesian(UnitCell *cell,
+ double ax, double ay, double az,
+ double bx, double by, double bz,
+ double cx, double cy, double cz);
+
+extern void cell_set_parameters(UnitCell *cell, double a, double b, double c,
+ double alpha, double beta, double gamma);
+
+#endif /* CELL_H */
diff --git a/src/image.c b/src/image.c
new file mode 100644
index 00000000..b3825b38
--- /dev/null
+++ b/src/image.c
@@ -0,0 +1,174 @@
+/*
+ * image.c
+ *
+ * Handle images and image features
+ *
+ * (c) 2006-2009 Thomas White <thomas.white@desy.de>
+ *
+ * template_index - Indexing diffraction patterns by template matching
+ *
+ */
+
+
+#define _GNU_SOURCE 1
+#include <stdlib.h>
+#include <assert.h>
+#include <math.h>
+
+#include "image.h"
+#include "utils.h"
+
+
+int image_add(ImageList *list, struct image *image)
+{
+ if ( list->images ) {
+ list->images = realloc(list->images,
+ (list->n_images+1)*sizeof(struct image));
+ } else {
+ assert(list->n_images == 0);
+ list->images = malloc(sizeof(struct image));
+ }
+
+ /* Copy the metadata */
+ list->images[list->n_images] = *image;
+
+ /* Fill out some extra fields */
+ list->images[list->n_images].features = NULL;
+ list->images[list->n_images].rflist = NULL;
+
+ list->n_images++;
+
+ return list->n_images - 1;
+}
+
+
+ImageList *image_list_new()
+{
+ ImageList *list;
+
+ list = malloc(sizeof(ImageList));
+
+ list->n_images = 0;
+ list->images = NULL;
+
+ return list;
+}
+
+
+void image_add_feature_reflection(ImageFeatureList *flist, double x, double y,
+ struct image *parent, double intensity)
+{
+ if ( flist->features ) {
+ flist->features = realloc(flist->features,
+ (flist->n_features+1)*sizeof(ImageFeature));
+ } else {
+ assert(flist->n_features == 0);
+ flist->features = malloc(sizeof(ImageFeature));
+ }
+
+ flist->features[flist->n_features].x = x;
+ flist->features[flist->n_features].y = y;
+ flist->features[flist->n_features].intensity = intensity;
+ flist->features[flist->n_features].parent = parent;
+ flist->features[flist->n_features].partner = NULL;
+ flist->features[flist->n_features].partner_d = 0.0;
+
+ flist->n_features++;
+
+}
+
+
+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;
+
+ flist = malloc(sizeof(ImageFeatureList));
+
+ flist->n_features = 0;
+ flist->features = NULL;
+
+ return flist;
+}
+
+
+void image_feature_list_free(ImageFeatureList *flist)
+{
+ if ( !flist ) return;
+
+ if ( flist->features ) free(flist->features);
+ free(flist);
+}
+
+
+ImageFeature *image_feature_closest(ImageFeatureList *flist, double x, double y,
+ double *d, int *idx)
+{
+ int i;
+ double dmin = +HUGE_VAL;
+ int closest = 0;
+
+ for ( i=0; i<flist->n_features; i++ ) {
+
+ double d;
+
+ d = distance(flist->features[i].x, flist->features[i].y, x, y);
+
+ if ( d < dmin ) {
+ dmin = d;
+ closest = i;
+ }
+
+ }
+
+ if ( dmin < +HUGE_VAL ) {
+ *d = dmin;
+ *idx = closest;
+ return &flist->features[closest];
+ }
+
+ *d = +INFINITY;
+ return NULL;
+}
+
+
+ImageFeature *image_feature_second_closest(ImageFeatureList *flist,
+ double x, double y, double *d,
+ int *idx)
+{
+ int i;
+ double dmin = +HUGE_VAL;
+ int closest = 0;
+ double dfirst;
+ int idxfirst;
+
+ image_feature_closest(flist, x, y, &dfirst, &idxfirst);
+
+ for ( i=0; i<flist->n_features; i++ ) {
+
+ double d;
+
+ d = distance(flist->features[i].x, flist->features[i].y, x, y);
+
+ if ( (d < dmin) && (i != idxfirst) ) {
+ dmin = d;
+ closest = i;
+ }
+
+ }
+
+ if ( dmin < +HUGE_VAL ) {
+ *d = dmin;
+ *idx = closest;
+ return &flist->features[closest];
+ }
+
+ return NULL;
+
+}
diff --git a/src/image.h b/src/image.h
new file mode 100644
index 00000000..f53f7d92
--- /dev/null
+++ b/src/image.h
@@ -0,0 +1,121 @@
+/*
+ * image.h
+ *
+ * Handle images and image features
+ *
+ * (c) 2006-2009 Thomas White <thomas.white@desy.de>
+ *
+ * template_index - Indexing diffraction patterns by template matching
+ *
+ */
+
+
+#ifdef HAVE_CONFIG_H
+#include <config.h>
+#endif
+
+#ifndef IMAGE_H
+#define IMAGE_H
+
+#include <stdint.h>
+
+
+typedef enum {
+ FORMULATION_CLEN,
+ FORMULATION_PIXELSIZE
+} FormulationMode;
+
+
+typedef struct imagefeature_struct {
+
+ struct image *parent;
+ double x;
+ double y;
+ double intensity;
+
+ /* Partner for this feature (in another feature list) or NULL */
+ struct imagefeature_struct *partner;
+
+ /* Distance between this feature and its partner, if any. */
+ double partner_d;
+
+ /* The reflection this was projected from, if any */
+ struct reflection_struct *reflection;
+
+} ImageFeature;
+
+
+typedef struct {
+
+ ImageFeature *features;
+ int n_features;
+
+} ImageFeatureList;
+
+
+struct image {
+
+ uint16_t *data;
+
+ /* Radians. Defines where the pattern lies in reciprocal space */
+ double tilt;
+
+ /* Radians. Defines where the pattern lies in reciprocal space */
+ double omega;
+
+ /* Image parameters can be given as camera length or pixel size.
+ * If FORMULATION_CLEN, then camera_len and resolution must be given.
+ * If FORMULATION_PIXELSIZE, then pixel_size only is needed.*/
+ FormulationMode fmode;
+ double pixel_size;
+ double camera_len;
+ double resolution; /* pixels per metre */
+
+ /* Wavelength must always be given */
+ double lambda;
+
+ int width;
+ int height;
+ double x_centre;
+ double y_centre;
+
+ ImageFeatureList *features; /* "Experimental" features */
+ ImageFeatureList *rflist; /* "Predicted" features */
+
+};
+
+
+typedef struct imagelist_struct {
+
+ int n_images;
+ struct image *images;
+
+} ImageList;
+
+
+extern ImageList *image_list_new(void);
+
+extern int image_add(ImageList *list, struct image *image);
+
+extern ImageFeatureList *image_feature_list_new(void);
+
+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 ImageFeature *image_feature_closest(ImageFeatureList *flist,
+ double x, double y, double *d,
+ int *idx);
+
+extern ImageFeature *image_feature_second_closest(ImageFeatureList *flist,
+ double x, double y, double *d,
+ int *idx);
+
+
+#endif /* IMAGE_H */
diff --git a/src/main.c b/src/main.c
index ffb269e4..79dff45a 100644
--- a/src/main.c
+++ b/src/main.c
@@ -60,6 +60,8 @@ int main(int argc, char *argv[])
return 1;
}
+ printf("Generating templates...\n");
+
printf("Input files (%i):\n", nin);
for ( i=0; i<nin; i++ ) {
printf("%6i: %s\n", i+1, in_files[i]);
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;
+}
diff --git a/src/relrod.h b/src/relrod.h
new file mode 100644
index 00000000..1dac6057
--- /dev/null
+++ b/src/relrod.h
@@ -0,0 +1,24 @@
+/*
+ * relrod.h
+ *
+ * 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
+ *
+ */
+
+#ifdef HAVE_CONFIG_H
+#include <config.h>
+#endif
+
+#ifndef RELROD_H
+#define RELROD_H
+
+#include "image.h"
+#include "cell.h"
+
+extern void get_reflections(struct image *image, UnitCell *cell);
+
+#endif /* RELROD_H */
diff --git a/src/sim-main.c b/src/sim-main.c
new file mode 100644
index 00000000..38ad45bc
--- /dev/null
+++ b/src/sim-main.c
@@ -0,0 +1,82 @@
+/*
+ * sim-main.c
+ *
+ * Simulate test data
+ *
+ * (c) 2006-2009 Thomas White <thomas.white@desy.de>
+ *
+ * template_index - Indexing diffraction patterns by template matching
+ *
+ */
+
+
+#ifdef HAVE_CONFIG_H
+#include <config.h>
+#endif
+
+#include <stdarg.h>
+#include <stdlib.h>
+#include <stdio.h>
+#include <string.h>
+#include <unistd.h>
+
+#include "image.h"
+#include "relrod.h"
+#include "cell.h"
+#include "utils.h"
+
+
+static void main_show_help(const char *s)
+{
+ printf("Syntax: %s <file1.h5> <file2.h5> [...]\n\n", s);
+ printf("Index diffraction patterns\n\n");
+ printf(" -h Display this help message\n");
+}
+
+
+int main(int argc, char *argv[])
+{
+ int c;
+ ImageList *list;
+ UnitCell *cell;
+ struct image image;
+
+ while ((c = getopt(argc, argv, "h")) != -1) {
+
+ switch ( c ) {
+
+ case 'h' : {
+ main_show_help(argv[0]);
+ return 0;
+ }
+
+ }
+
+ }
+
+ printf("Generating test data...\n");
+ list = image_list_new();
+ image.width = 512;
+ image.height = 512;
+ 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.data = malloc(512*512*2);
+ image_add(list, &image);
+
+ cell = cell_new_from_parameters(1.0,
+ 1.0,
+ 1.0,
+ deg2rad(90.0),
+ deg2rad(90.0),
+ deg2rad(90.0));
+
+ get_reflections(&image, cell);
+
+ return 0;
+}
diff --git a/src/utils.c b/src/utils.c
new file mode 100644
index 00000000..0c7d7ad4
--- /dev/null
+++ b/src/utils.c
@@ -0,0 +1,124 @@
+/*
+ * utils.c
+ *
+ * Utility stuff
+ *
+ * (c) 2006-2009 Thomas White <thomas.white@desy.de>
+ *
+ * template_index - Indexing diffraction patterns by template matching
+ *
+ */
+
+#include <math.h>
+#include <string.h>
+
+#include "utils.h"
+
+
+/* Return the MOST POSITIVE of two numbers */
+unsigned int biggest(signed int a, signed int b)
+{
+ if ( a>b ) {
+ return a;
+ }
+ return b;
+}
+
+
+/* Return the LEAST POSITIVE of two numbers */
+unsigned int smallest(signed int a, signed int b)
+{
+ if ( a<b ) {
+ return a;
+ }
+ return b;
+}
+
+
+double distance(double x1, double y1, double x2, double y2)
+{
+ return sqrt((x2-x1)*(x2-x1) + (y2-y1)*(y2-y1));
+}
+
+
+double modulus(double x, double y, double z)
+{
+ return sqrt(x*x + y*y + z*z);
+}
+
+
+double modulus_squared(double x, double y, double z) {
+ return x*x + y*y + z*z;
+}
+
+
+double distance3d(double x1, double y1, double z1,
+ double x2, double y2, double z2)
+{
+ return modulus(x1-x2, y1-y2, z1-z2);
+}
+
+
+/* Angle between two vectors. Answer in radians */
+double angle_between(double x1, double y1, double z1,
+ double x2, double y2, double z2)
+{
+ double mod1 = modulus(x1, y1, z1);
+ double mod2 = modulus(x2, y2, z2);
+ return acos( (x1*x2 + y1*y2 + z1*z2) / (mod1*mod2) );
+}
+
+
+/* As above, answer in degrees */
+double angle_between_d(double x1, double y1, double z1,
+ double x2, double y2, double z2)
+{
+ return rad2deg(angle_between(x1, y1, z1, x2, y2, z2));
+}
+
+
+/* Wavelength of an electron (in m) given accelerating potential (in V) */
+double lambda(double V)
+{
+ double m = 9.110E-31;
+ double h = 6.625E-34;
+ double e = 1.60E-19;
+ double c = 2.998E8;
+
+ return h / sqrt(2*m*e*V*(1+((e*V) / (2*m*c*c))));
+}
+
+
+size_t skipspace(const char *s)
+{
+ size_t i;
+
+ for ( i=0; i<strlen(s); i++ ) {
+ if ( (s[i] != ' ') && (s[i] != '\t') ) return i;
+ }
+
+ return strlen(s);
+}
+
+
+void chomp(char *s)
+{
+ size_t i;
+
+ if ( !s ) return;
+
+ for ( i=0; i<strlen(s); i++ ) {
+ if ( (s[i] == '\n') || (s[i] == '\r') ) {
+ s[i] = '\0';
+ return;
+ }
+ }
+}
+
+
+int sign(double a)
+{
+ if ( a < 0 ) return -1;
+ if ( a > 0 ) return +1;
+ return 0;
+}
diff --git a/src/utils.h b/src/utils.h
new file mode 100644
index 00000000..ef37bc73
--- /dev/null
+++ b/src/utils.h
@@ -0,0 +1,42 @@
+/*
+ * utils.h
+ *
+ * Utility stuff
+ *
+ * (c) 2006-2009 Thomas White <thomas.white@desy.de>
+ *
+ * template_index - Indexing diffraction patterns by template matching
+ *
+ */
+
+#ifndef UTILS_H
+#define UTILS_H
+
+#ifdef HAVE_CONFIG_H
+#include <config.h>
+#endif
+
+#include <math.h>
+
+extern unsigned int biggest(signed int a, signed int b);
+extern unsigned int smallest(signed int a, signed int b);
+extern double distance(double x1, double y1, double x2, double y2);
+extern double modulus(double x, double y, double z);
+extern double modulus_squared(double x, double y, double z);
+extern double angle_between(double x1, double y1, double z1,
+ double x2, double y2, double z2);
+extern double angle_between_d(double x1, double y1, double z1,
+ double x2, double y2, double z2);
+extern double lambda(double voltage);
+extern double distance3d(double x1, double y1, double z1,
+ double x2, double y2, double z2);
+extern size_t skipspace(const char *s);
+extern void chomp(char *s);
+extern int sign(double a);
+
+#define rad2deg(a) ((a)*180/M_PI)
+#define deg2rad(a) ((a)*M_PI/180)
+
+#define is_odd(a) ((a)%2==1)
+
+#endif /* UTILS_H */
generated by cgit v1.2.3 (git 2.25.1) at 2024-07-02 12:47:46 +0000