1 files changed, 126 insertions, 0 deletions
diff --git a/src/multislice.c b/src/multislice.c
new file mode 100644
index 0000000..b314098
--- /dev/null
+++ b/src/multislice.c
@@ -0,0 +1,126 @@
+/*
+ * multislice.c
+ *
+ * Multislice Dynamical Simulations
+ *
+ * (c) 2006-2007 Thomas White <taw27@cam.ac.uk>
+ *
+ *  synth2d - Two-Dimensional Crystallographic Fourier Synthesis
+ *
+ */
+
+#ifdef HAVE_CONFIG_H
+#include <config.h>
+#endif
+
+#include <stdlib.h>
+#include <math.h>
+#include <stdio.h>
+#include <fftw3.h>
+#include <stdlib.h>
+#include <string.h>
+
+#include "reflist.h"
+#include "model.h"
+#include "data.h"
+
+/* Relativistic Electron Interaction Constant.  Voltage in volts. */
+static double multislice_interaction(double voltage) {
+	
+	double a, b;
+	
+	a = (1 + voltage*1.9569341e-6) * 0.25615739;
+	b = sqrt(voltage) * sqrt(1 + voltage*0.97846707e-6);
+	
+	return a / b;
+	
+}
+
+/* Unit cell volume */
+static double multislice_volume(double a, double b, double c, double alpha, double beta, double gamma) {
+
+	
+
+}
+
+static fftw_complex *multislice_phasegrating(AtomicModel *model_orig, size_t width, size_t height) {
+
+	AtomicModel	*model;
+	fftw_complex	*phase_grating_real;
+	fftw_complex	*phase_grating_reciprocal;
+	fftw_plan	plan;
+	ReflectionList	*V;
+	ReflectionList	*template;
+	
+	template = reflist_new();
+	
+	model = model_copy(model_orig);
+	model->thickness = 0.0;
+	V = model_calculate_f(template, model, 69);
+	model_free(model);
+	reflist_free(template);
+	
+	plan = fftw_plan_dft_2d(width, height, phase_grating_real, phase_grating_reciprocal, FFTW_FORWARD, FFTW_MEASURE);
+	fftw_execute(plan);
+	fftw_destroy_plan(plan);
+	
+	return phase_grating_reciprocal;
+	
+}
+
+static fftw_complex *multislice_propogator(AtomicModel *model_orig, size_t width, size_t height, double slice) {
+
+	
+
+	return NULL;
+
+}
+
+static void multislice_multiply(fftw_complex *a, fftw_complex *b, size_t size) {
+	size_t i;
+	for ( i=0; i<size; i++ ) {
+		a[i][0] *= b[i][0];
+		a[i][1] *= b[i][1];
+	}
+}
+
+static void multislice_convolve(fftw_complex *a, fftw_complex *b, size_t size) {
+
+	
+
+}
+
+/* Calculate dynamical diffraction amplitudes and phases */
+ReflectionList *multislice_calculate_f_dyn(AtomicModel *model, ReflectionList *template) {
+
+	size_t		width, height;
+	fftw_complex	*wavefunction;
+	fftw_complex	*phasegrating;
+	fftw_complex	*propogator;
+	double		t;
+	double		c;
+	
+	c = data_c();
+	width = 128;  height = 128;
+	
+	printf("MS: Calculating phase grating function...\n"); fflush(stdout);
+	phasegrating = multislice_phasegrating(model, width, height);
+	
+	printf("MS: Calculating propogration function...\n"); fflush(stdout);
+	propogator = multislice_propogator(model, width, height, c);
+	
+	/* Initial value of wavefunction */
+	wavefunction = fftw_malloc(width*height*sizeof(fftw_complex));
+	memcpy(wavefunction, phasegrating, width*height*sizeof(fftw_complex));
+	
+	/* Iterate */
+	for ( t=0.0; t<model->thickness; t+=c ) {
+		multislice_multiply(wavefunction, propogator, width*height);
+		multislice_convolve(wavefunction, phasegrating, width*height);
+		printf("MS: Current thickness: %f\n", t); fflush(stdout);
+	}
+	
+	return NULL;
+
+}
+