#include<complex.h>
#include<fftw3.h>
#include<vector>
#include<iostream>
#include<fstream>
#include<cmath>
#include<cstdlib>

using namespace std;

const int N = 512;
const double DX = 0.06;

void init_vector(double _Complex* v, int n, double dx, double lambda) {
	double k = 2.0*M_PI/lambda;
	for(int i=0; i<n; i++) {
		double x = i*dx;
		v[i] = cos(k*x) + I*0.0;
	}
}

void init_gaussian(double _Complex* v, int n, double dx, double A, double x0, double w) {
	for(int i=0; i<n; i++) {
		double x = i*dx;
		v[i] = A*exp(-pow(x-x0,2)/pow(w,2));
	}
}

void write_vector(ofstream& f, double _Complex* v, int n, double dx) {
	for(int i=0; i<n; i++)
		f << i*dx << " " << creal(v[i]) << " " << cimag(v[i]) << endl;
	f << endl;
}

void derivar_kspace(double _Complex* v, double _Complex* vout, int n, double dx) {
}

int main() {
	double _Complex *phi = new double _Complex[N];
	double _Complex *phi_k = new double _Complex[N];
	ofstream fout("fft.dat");

	fftw_plan transformar = fftw_plan_dft_1d(N, phi, phi_k, FFTW_FORWARD, FFTW_MEASURE);
	fftw_plan antitransformar = fftw_plan_dft_1d(N, phi_k, phi, FFTW_BACKWARD, FFTW_MEASURE);

// 	init_vector(phi, N, DX, 2.0);
	init_gaussian(phi, N, DX, 1.0, DX*N/2, 1.0);

	write_vector(fout, phi, N, DX);
	
	fftw_execute(transformar);
	for(int i=0;i<N;i++) phi_k[i] = phi_k[i]*1.0/N;
	fftw_execute(antitransformar);

	write_vector(fout, phi, N, DX);
	
	fout.close();
	fftw_destroy_plan(antitransformar);
	fftw_destroy_plan(transformar);
	delete(phi);
	delete(phi_k);
	return 0;
}