/************************************************************************************************
*												*
*	Ginz1D.c										*
* 												*
*												*
*	Solves the 1D Ginzburg-Landau-equation using a pseudospectral scheme in spatial domain	*
*	and a exponential time-differencing scheme in temporal domain.				*
*	The program displays its results as a space-time plot 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			*
*												*
*												*
*	compile with										*
*		gcc -std=c99 -O3 -Wall -o Ginz1D Ginz1D.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 = 512;		//Grid resolution
const double L = 200;		//Physical length
const double alpha = 0.0;	//coefficients of GL-equation
const double beta  = -4.0;	//             -"-

const int MAX_ITER = 10000;	//maximal number of iterations;
const int EVERY_OUT = 5;	//use every 'EVERY_OUT'-th timestep for display
const int EVERY_DISLIN = 10;	//display, once 'EVERY_DISLIN' timesteps are gathered

//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
complex *A_display;			//2D field for the space-time plot

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; 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));
}

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

//assigns one array to the other
void assignArrays(complex *out, const complex *in) {
	for(int i=0; i<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; 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; i++) {
		complex q = (-(1.0+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(const complex *in) {
	//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(in[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(in[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);
	Atmp = fftw_malloc(sizeof(complex)*N);

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

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

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

	A_display = fftw_malloc(sizeof(complex)*N*N);
	//its a space-time plot that fills constantly, so we init it to 0
	for(int i=0; i<N*N; i++)
		A_display[i] = 0;


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

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

	//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; 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
		/****************************************************************************************
		 *	now heres the tricky part regarding the space-time plot:			*
		 *											*
		 *	i wanted the space-time plot to 'scroll up' as the time passes.			*
		 *	so first of all i move every timestep (i.e. every line in x-dir) up by one step	*
		 *	(this is done in the 2 for-loops)						*
		 *	after that, the current timestep is copied into the last line to get the	*
		 *	'scrolling-effect'								*
		 *	(of course we have to switch from fourier to spatial domain first because we 	*
		 *	want to	see the solutions in the spatial domain)				*
		 ***************************************************************************************/
		if(iter%EVERY_OUT==0) {
			//save timestep
			for(int i=0; i<N; i++) {
				for(int j=N-1; j>=1; j--) {
					A_display[i*N+j] = A_display[i*N+j-1];
				}
			}
			assignArrays(Atmp, A);
			fftF2S(Atmp);			
			for(int i=0; i<N; i++) {
				A_display[i*N+0] = Atmp[i];
			}
			//because cause dislin is pretty slow, we just display
			//once 'EVERY_DISLIN' timesteps are gathered
			if((iter/EVERY_OUT)%EVERY_DISLIN==0)
				display(A_display);
			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);
	free(A_display);

	fftw_destroy_plan(planF2S);
	fftw_destroy_plan(planS2F);

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