/****************************************************************************************************************
*	Program to demonstrate the usage of the FFTW-library to calculate derivatives.				*
*														*
*	compile with												*
*		gcc -std=c99 -o testTrafo testTrafo.c -I(FFTW3_PATH)/include -L(FFTW_PATH)/lib -lfftw3 -lm	*
*														*
*	author: Johannes Luelff											*
*														*
****************************************************************************************************************/

#include <stdio.h>
#include <math.h>
#include "fftw3.h"

const double PI = 3.141592635;
const int N = 256;
const double L = 1.0;

// 0 <= x < L
double init(double x) {
	return sin(16*PI*x/L);
}


int main() {

	int i;
	double field_real_in[N], field_real_out[N];
	fftw_complex field_complex[N/2+1];
	
	fftw_plan planR2C = fftw_plan_dft_r2c_1d(N, field_real_in, 
		field_complex, FFTW_ESTIMATE);
	fftw_plan planC2R = fftw_plan_dft_c2r_1d(N, field_complex, 
		field_real_out, FFTW_ESTIMATE);

	for(i = 0; i<N; i++) 
		field_real_in[i] = init(L*i/(double)N);

	fftw_execute(planR2C);
	//now in complex space
	for(i=0; i<N/2+1; i++) {
		field_complex[i][0] /= sqrt(N);
		field_complex[i][1] /= sqrt(N);
	}
	
	//save FT
	FILE *f1 = fopen("complex.txt", "w");
	for(i=0; i<N/2+1; i++)
		fprintf(f1, "%f %f %f\n", 2*PI*i/(double)L, 
			field_complex[i][0], field_complex[i][1]);

	fftw_execute(planC2R);
	//back in real space
	for(i=0; i<N; i++)
		field_real_out[i] /= sqrt(N);

	FILE *f2 = fopen("real.txt", "w");
	for(i=0; i<N; i++) 
		fprintf(f2, "%f %f %f\n", L*i/(double)N, 
			field_real_in[i], field_real_out[i]);

	fclose(f1);
	fclose(f2);

	fftw_destroy_plan(planR2C);
	fftw_destroy_plan(planC2R);

	printf("READY\n");

	return 0;
}