Restructuring ready for new simulation method

6 files changed, 85 insertions, 232 deletions

diff --git a/src/Makefile.am b/src/Makefile.am
index d30736c7..3ab2d7f8 100644
--- a/src/Makefile.am
+++ b/src/Makefile.am
@@ -1,5 +1,4 @@
 bin_PROGRAMS = pattern_sim
 
-pattern_sim_SOURCES = main.c relrod.c utils.c image.c cell.c hdf5-file.c \
-                      templates.c
+pattern_sim_SOURCES = main.c diffraction.c utils.c image.c cell.c hdf5-file.c
 pattern_sim_LDADD = @LIBS@ -lm -lgsl -lgslcblas
diff --git a/src/diffraction.c b/src/diffraction.c
new file mode 100644
index 00000000..93c6a22c
--- /dev/null
+++ b/src/diffraction.c
@@ -0,0 +1,59 @@
+/*
+ * diffraction.c
+ *
+ * Calculate diffraction patterns by Fourier methods
+ *
+ * (c) 2007-2009 Thomas White <thomas.white@desy.de>
+ *
+ * pattern_sim - Simulate diffraction patterns from small crystals
+ *
+ */
+
+
+#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_diffraction(struct image *image, UnitCell *cell)
+{
+
+}
diff --git a/src/diffraction.h b/src/diffraction.h
new file mode 100644
index 00000000..c32e62c2
--- /dev/null
+++ b/src/diffraction.h
@@ -0,0 +1,24 @@
+/*
+ * diffraction.h
+ *
+ * Calculate diffraction patterns by Fourier methods
+ *
+ * (c) 2007-2009 Thomas White <thomas.white@desy.de>
+ *
+ * pattern_sim - Simulate diffraction patterns from small crystals
+ *
+ */
+
+#ifdef HAVE_CONFIG_H
+#include <config.h>
+#endif
+
+#ifndef DIFFRACTION_H
+#define DIFFRACTION_H
+
+#include "image.h"
+#include "cell.h"
+
+extern void get_diffraction(struct image *image, UnitCell *cell);
+
+#endif	/* DIFFRACTION_H */
diff --git a/src/main.c b/src/main.c
index 04f806f2..8ab11c85 100644
--- a/src/main.c
+++ b/src/main.c
@@ -85,24 +85,7 @@ int main(int argc, char *argv[])
 		memset(image.data, 0, 512*512*2);
 		image.tilt = deg2rad(t);
 
-		/* Calculate reflections */
-		get_reflections(&image, cell, 1.0/CRYSTAL_SIZE);
-
-		/* Construct the image */
-		nrefl = image_feature_count(image.rflist);
-		for ( i=0; i<nrefl; i++ ) {
-
-			struct imagefeature *f;
-			int x, y;
-
-			f = image_get_feature(image.rflist, i);
-
-			x = f->x;
-			y = f->y;  /* Discards digits after the decimal point */
-
-			image.data[y*image.width+x] = 1;
-
-		}
+		get_diffraction(&image, cell);
 
 		/* Write the output file */
 		snprintf(filename, 32, "simulated-%.0f.h5", t);
diff --git a/src/templates.c b/src/templates.c
deleted file mode 100644
index 1ad1c95a..00000000
--- a/src/templates.c
+++ /dev/null
@@ -1,183 +0,0 @@
-/*
- * templates.c
- *
- * Handle templates
- *
- * (c) 2006-2009 Thomas White <thomas.white@desy.de>
- *
- * pattern_sim - Simulate diffraction patterns from small crystals
- *
- */
-
-
-#define _GNU_SOURCE 1
-#include <stdlib.h>
-#include <assert.h>
-#include <math.h>
-#include <stdio.h>
-
-#include "templates.h"
-#include "cell.h"
-#include "image.h"
-#include "utils.h"
-#include "relrod.h"
-
-
-struct _templatelist
-{
-	int		n_templates;
-	struct template	*templates;
-};
-
-
-struct template_feature
-{
-	float		x;
-	float		y;
-	struct template_feature *next;
-};
-
-
-struct template
-{
-	float		omega;
-	float		tilt;
-	struct template_feature *features;
-};
-
-
-static int template_add(TemplateList *list, struct template *template)
-{
-	if ( list->templates ) {
-		list->templates = realloc(list->templates,
-		                 (list->n_templates+1)*sizeof(struct template));
-	} else {
-		assert(list->n_templates == 0);
-		list->templates = malloc(sizeof(struct template));
-	}
-
-	/* Copy the data */
-	list->templates[list->n_templates] = *template;
-
-	list->n_templates++;
-
-	return list->n_templates - 1;
-}
-
-
-static TemplateList *template_list_new()
-{
-	TemplateList *list;
-
-	list = malloc(sizeof(TemplateList));
-
-	list->n_templates = 0;
-	list->templates = NULL;
-
-	return list;
-}
-
-
-TemplateList *generate_templates(UnitCell *cell, struct image params)
-{
-	TemplateList *list;
-	double omega, tilt;
-
-	list = template_list_new();
-
-	omega = deg2rad(40.0);
-
-	//for ( omega=deg2rad(-180); omega<deg2rad(180); omega+=deg2rad(1) ) {
-
-		params.omega = omega;
-
-		for ( tilt=0; tilt<deg2rad(180); tilt+=deg2rad(1) ) {
-
-			struct template t;
-			struct template_feature *tfc;
-			int nrefl, i;
-
-			t.omega = omega;
-			t.tilt = tilt;
-			t.features = malloc(sizeof(struct template_feature));
-			t.features->next = NULL;
-			tfc = t.features;
-
-			params.tilt = tilt;
-			get_reflections(&params, cell, 0.01e9);
-
-			nrefl = image_feature_count(params.rflist);
-			for ( i=0; i<nrefl; i++ ) {
-
-				struct imagefeature *f;
-
-				f = image_get_feature(params.rflist, i);
-
-				tfc->x = f->x;
-				tfc->y = f->y;
-				tfc->next = malloc(sizeof(struct template_feature));
-				tfc = tfc->next;
-
-			}
-
-			template_add(list, &t);
-
-			image_feature_list_free(params.rflist);
-
-			printf("Generating templates... %+5.2f %+5.2f\r",
-			       rad2deg(omega), rad2deg(tilt));
-
-		}
-	//}
-
-	printf("Generating templates... done                      \n");
-
-	return list;
-}
-
-
-static int try_template(struct image *image, struct template template)
-{
-	int fit = 0;
-	struct template_feature *f;
-
-	f = template.features;
-	while ( f != NULL ) {
-
-		int x, y;
-
-		x = f->x;
-		y = f->y;	/* Discards digits after the decimal point */
-
-		fit += image->data[y*image->width+x];
-
-		f = f->next;
-
-	}
-
-	return fit;
-}
-
-
-int try_templates(struct image *image, TemplateList *list)
-{
-	int i;
-	int fit_max = 0;
-	int i_max = 0;
-
-	for ( i=0; i<list->n_templates; i++ ) {
-
-		int fit;
-
-		fit = try_template(image, list->templates[i]);
-		if ( fit > fit_max ) {
-			fit_max = fit;
-			i_max = i;
-		}
-
-	}
-	image->omega = list->templates[i_max].omega;
-	image->tilt = list->templates[i_max].tilt;
-
-	return 0;
-}
diff --git a/src/templates.h b/src/templates.h
deleted file mode 100644
index c9d053ad..00000000
--- a/src/templates.h
+++ /dev/null
@@ -1,29 +0,0 @@
-/*
- * templates.h
- *
- * Handle templates
- *
- * (c) 2006-2009 Thomas White <thomas.white@desy.de>
- *
- * pattern_sim - Simulate diffraction patterns from small crystals
- *
- */
-
-
-#ifdef HAVE_CONFIG_H
-#include <config.h>
-#endif
-
-#ifndef TEMPLATES_H
-#define TEMPLATES_H
-
-#include "cell.h"
-#include "image.h"
-
-/* An opaque type representing a list of templates */
-typedef struct _templatelist TemplateList;
-
-extern TemplateList *generate_templates(UnitCell *cell, struct image params);
-extern int try_templates(struct image *image, TemplateList *list);
-
-#endif	/* TEMPLAETS_H */