/************************************************************************************************
*												*
*	Ginz2D.c										*
* 												*
*												*
*	Solves the 2D Ginzburg-Landau-equation using a pseudospectral scheme in spatial domain	*
*	and a exponential time-differencing scheme in temporal domain.				*
*	The program displays its results using dislin.						*
*												*
*												*
*	Some parameter sets:									*
*		(alpha, beta)	=	( 1.0,  2.0)	//plain waves				*
*				=	( 2.0, -1.0)	//phase turbulence			*
*				=	( 2.0, -2.0)	//defect turbulence			*
*				=	( 0.5, -1.5)	//intermittency				*
*				=	( 0.0, -4.0)	//intermittency				*
*				=	( 0.0,  1.5)	//coherent structures, spirals		*
*												*
*		These parameters are actually borrowed from the 1D simulation, but it shows 	*
*		that the parameters produce similar structures in 2D.				*
*												*
*												*
*	compile with										*
*		gcc -std=c99 -O3 -Wall -o Ginz2D Ginz2D.c -I(FFTW3_PATH)/include		*
*			-L(FFTW_PATH)/lib -lfftw3 -I(DISLIN_PATH) -L(DISLIN_PATH) -ldislin -lm	*
*												*
*												*
*	author: Johannes Luelff	(johannes.luelff@uni-muenster.de)				*
*												*
************************************************************************************************/


#include <stdio.h>
#include <stdlib.h>
//we want to use the C99 complex data type, so we need complex.h
#include <complex.h>
#include <math.h>
#include "fftw3.h"
#include "dislin.h"


//these two defines control whether we...

//... want to use dislin (we surely want to, cause otherwise
//    this program wouldnt produce any output)
#define USE_DISLIN
//... want dislin to print its result into a PNG file or onto the screen
#undef USE_PNG


//constants etc.
#define PI 3.14159265358979323846264338327950288419716939937510582097494459230781
const int N = 256;		//Grid resolution
const double L = 200;		//Physical length
const double alpha = 0.0;	//coefficients of GL-equation
const double beta  = 1.5;	//             -"-

const int MAX_ITER = 1000;	//maximal number of iterations;
const int EVERY_DISPLAY = 5;	//display every 'Every_DISPLAY'-th timestep

//time(step)
double t = 0;
double dt = 0.05;

//arrays etc.
complex *A;			//the 'main' field
complex *Aedt_N1;		//field for nonlinearity...
complex *Aedt_N2;		//... and backup field for nonlin. from previous timestep
complex *Atmp;
double *kabsvec;		//vector that stores |k|**2

float *A_dislin_left, *A_dislin_right;	//used for display with dislin

fftw_plan planS2F, planF2S;		//plan for spatial to fourier, fourier to spatial


//multiplies array with scalar
void smultArray(complex *arr, double s) {
	for(int i=0; i<N*N; i++)
		arr[i] *= s;
}

//fft from spatial to fourier domain
void fftS2F(complex *arr) {
	fftw_execute_dft(planS2F, arr, arr);
	smultArray(arr, 1/sqrt(N*N));
}

//fft from fourier to spatial domain
void fftF2S(complex *arr) {
	fftw_execute_dft(planF2S, arr, arr);
	smultArray(arr, 1/sqrt(N*N));
}

//assigns one array to the other
void assignArrays(complex *out, const complex *in) {
	for(int i=0; i<N*N; i++)
		out[i] = in[i];
}

//returns |z|**2
double cnorm(const complex z) {
	return creal(z)*creal(z) + cimag(z)*cimag(z);
}


//calculates the nonlinear part of the GL-equation in spatial domain
void calcNonlin(const complex *in, complex *out) {
	assignArrays(out, in);
	fftF2S(out);

	for(int i=0; i<N*N; i++)
		out[i] = -(1.0+I*beta) * cnorm(out[i]) * out[i];

	fftS2F(out);
}

//calculates the next time step using a 2nd order exponential time differencing scheme
void doTimeStep_EDT2() {
	//switch nonlinearities N1 / N2
	complex *ctmp = Aedt_N1;
	Aedt_N1 = Aedt_N2;
	Aedt_N2 = ctmp;

	//calc new nonlinearity
	calcNonlin(A, Aedt_N1);

	//put all together according to EDT2
	//this may look scary, but its just the formula for EDT2 in a quite straight-forward manner.
	//q is the factor of the linear part of the GL-eq.
	//since we use a pseudospectral scheme, the 2nd derivative is also just a factor
	for(int i=0; i<N*N; i++) {
		complex q = (-(1+I*alpha)*kabsvec[i]+1);
		A[i] = A[i] * cexp(q*dt)
			+ Aedt_N1[i] * ((1+dt*q)*cexp(dt*q) - 1 - 2*dt*q) / (dt*q*q)
			+ Aedt_N2[i] * (-cexp(dt*q)+1+dt*q) / (dt*q*q);
	}

	t += dt;
}


//displays real part and absolute value of the field using dislin
//dont take all the numbers regarding dislin in this routine as 'the only truth';
//dealing with dislin is always kind of tricky and try-and-error-like
void display(complex *in) {
	//to spatial domain
	assignArrays(Atmp, in);
	fftF2S(Atmp);

	//prepare output array for dislin (left panel)
	double min_left = HUGE_VAL;
	double max_left = -HUGE_VAL;
	for(int i=0; i<N*N; i++) {
			A_dislin_left[i] = creal(Atmp[i]);

			if(A_dislin_left[i]<min_left)
				min_left = A_dislin_left[i];
			if(A_dislin_left[i]>max_left)
				max_left = A_dislin_left[i];
	}
	if(fabs(min_left-max_left)<1e-6) {
		min_left=-0.1;
		max_left= 0.1;
	}

	//prepare output array for dislin (right panel)
	double min_right = HUGE_VAL;
	double max_right = -HUGE_VAL;
	for(int i=0; i<N*N; i++) {
			A_dislin_right[i] = cabs(Atmp[i]);

			if(A_dislin_right[i]<min_right)
				min_right = A_dislin_right[i];
			if(A_dislin_right[i]>max_right)
				max_right = A_dislin_right[i];
	}
	if(fabs(min_right-max_right)<1e-6) {
		min_right=-0.1;
		max_right= 0.1;
	}

	//plot with dislin
#ifdef USE_PNG
	//init dislin
	metafl("PNG");
	setfil("out.png");
	scrmod("REVERS");
	winsiz(1280,640);
	page(2970*2,2970);
	disini();
	ax3len(2970*0.75, 2970*0.75, 2970*0.75);
#endif //USE_PNG

	erase();
	axspos(2970*0.1,2970*0.85);
	graf3(0,L,0,L/2, 0,L,0,L/2, min_left, max_left, min_left, (max_left-min_left)/2);
	crvmat(A_dislin_left, N, N, ((512/N)<=0)?1:(512/N), ((512/N)<=0)?1:(512/N));
	titlin("Real(A)", 3);
	title();
	endgrf();

	axspos(2970*1.1,2970*0.85);
	graf3(0,L,0,L/2, 0,L,0,L/2, min_right, max_right, min_right, (max_right-min_right)/2);
	crvmat(A_dislin_right, N, N, ((512/N)<=0)?1:(512/N), ((512/N)<=0)?1:(512/N));
	titlin("Abs(A)", 3);
	title();
	endgrf();

#ifdef USE_PNG
	disfin();
#endif //USE_PNG
}


int main(int argc, char *argv[]) {

	//get memory
	A = fftw_malloc(sizeof(complex)*N*N);
	Atmp = fftw_malloc(sizeof(complex)*N*N);

	Aedt_N1 = fftw_malloc(sizeof(complex)*N*N);
	Aedt_N2 = fftw_malloc(sizeof(complex)*N*N);

	kabsvec = fftw_malloc(sizeof(double)*N*N);

	A_dislin_left = fftw_malloc(sizeof(float)*N*N);
	A_dislin_right = fftw_malloc(sizeof(float)*N*N);


	//make plans for fftw
	planS2F = fftw_plan_dft_2d(N, N, (fftw_complex*)A, (fftw_complex*)A, FFTW_FORWARD,  FFTW_ESTIMATE);
	planF2S = fftw_plan_dft_2d(N, N, (fftw_complex*)A, (fftw_complex*)A, FFTW_BACKWARD, FFTW_ESTIMATE);

	//init kvec
	for(int i=0; i<N; i++)
		for(int j=0; j<N; j++) {
			double kx = 2*PI * ((i<N/2+1) ? (i) : (i-N)) / L;
			double ky = 2*PI * ((j<N/2+1) ? (j) : (j-N)) / L;
			kabsvec[i*N+j] = kx*kx + ky*ky;
		}

	//init field with white noise of amplitude 0.01
	//and go to fourier space, cause we will do the timestepping in fourier space
	for(int i=0; i<N*N; i++)
		A[i] = 0.02*(0.5 - rand()/(double)RAND_MAX + I*(0.5 - rand()/(double)RAND_MAX));
	fftS2F(A);


#ifdef USE_DISLIN
#ifndef USE_PNG
	//init dislin
	metafl("XWIN");
	scrmod("REVERS");
	winsiz(1280,640);
	page(2970*2,2970);
	disini();
	ax3len(2970*0.75, 2970*0.75, 2970*0.75);
#endif //USE_PNG
#endif //USE_DISLIN

	//main timestepping loop
	for(int iter=0; iter<=MAX_ITER; iter++) {

		doTimeStep_EDT2();

#ifdef USE_DISLIN
		if(iter%EVERY_DISPLAY==0) {
			display(A);
			printf("t = %e , dt = %e\n", t, dt);
		}
#endif //USE_DISLIN

	}

#ifdef USE_DISLIN
#ifndef USE_PNG
	//finalize dislin
	disfin();
#endif //USE_PNG
#endif //USE_DISLIN

	//clean up stuff
	free(A);
	free(Atmp);
	free(Aedt_N1);
	free(Aedt_N2);
	free(kabsvec);
	free(A_dislin_left);
	free(A_dislin_right);

	fftw_destroy_plan(planF2S);
	fftw_destroy_plan(planS2F);

	//yay were done
	printf("READY\n");
}