/*
 * 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;

}