/* Calculates time evolution of initial data (u, u_t)=(A*r^5*Gaussian, B*r^5*Gaussian) for cubic focusing Klein-Gordon equation 
u_tt-u_rr=-u+u^3/r^2 in radial 3d using centered time and centered space discretization.
Saves (A,B) in blowup.dat, dispersion.dat, or indecisive.dat depending on outcome */

/*****************************************************
** Parameters ****************************************
******************************************************/

/* minimum amplitude A */
#define MINA 0.0

/* maximum amplitude A */
#define MAXA 4.5

/* stepsize A */
#define INCA 0.007

/* minimum amplitude B */
#define MINB -5.0

/* maximum amplitude B */
#define MAXB 5.0

/* stepsize B */
#define INCB 0.007


/* Grid mesh dr in multiples of DRBASE, see below */
#define DRMULT 30


/* Number of spatial gridpoints */
#define KD 6000

/* Courant constant dt/dr */
#define R 0.9



/* Energy of Q */

#define Qenergy 1.5038


/*********************************************************/

#include <stdio.h>
#include <stdlib.h>
#include <math.h>


/* define resolution of soliton file */
#define DRBASE 1.0e-4

/* define grid mesh */
#define DR ((double)DRMULT*DRBASE)

/* global variable for Q */
double *Q;
 
/****************** initial data ******************/

/* initial position */

double f(long k, double A)
{
	double r=k*DR; /* compute radial coordinate from gridpoint */
	
	
	return A*r*sin(10.0*r)*exp(-r); 
}

/* initial velocity */

double g(long k, double B)
{
	
	double r=k*DR; /* compute radial coordinate from gridpoint */
	
	
	return B*r*exp(-r); 
}


/***********************************************/





/* import routine for soliton (or other functions) 
 data are assumed to be given with dr=DRBASE
 caller has to provide pointer to array */
int import(char *file, double *data, long size)
{
	register long k, j; /* loop variables */
	double dummy; /* auxiliary variable to store skipped entries */
	FILE *fp; /* file pointer */
	
	/* try to open file for reading */
	if((fp=fopen(file, "r"))==NULL)
		return 0; /* no success, return error */
	
	/* loop for reading data */
	for(k=0; k<size; k++)
	{
		/* check if end of file has been reached */
		if(feof(fp))
			data[k]=0.0; /* if yes, set 0.0 */
		else
		{
			/* if no, read data */
			fscanf(fp, "%lf", &data[k]);
			
			/* skip data if DRMULT>1 */
			for(j=1; j<DRMULT; j++)
			{
				/* read dummy data */
				fscanf(fp, "%lf", &dummy);
			}
		}
	}
	
	/* close file */
	fclose(fp);
	
	/* exit with success */
	return 1;
}

/* discretization of 1d Laplacian */
double D(double u[], long k, long K)
{
	return (u[k+1]-2.0*u[k]+u[k-1])/(DR*DR);
}


/* definition of the nonlinearity on the right-hand side of the equation */
double F(double u[], long k)
{
	return -u[k]+u[k]*u[k]*u[k]/(k*DR*k*DR);
}


/* calculates the free energy 1/2 integral u_t^2+(u_r-u/r)^2+u^2
 by using trapezoidal rule and second order accurate approximations to derivatives 
 needs one time step in the past and one in the future */
double kineticenergy(double u2[], double u1[], double u0[], long K, double dt)
{
	double en=0.0; /* stores 2 times the energy */
	register long k; /* loop variable */
	
	/* loop over all gridpoints
	 contribution from k=0 vanishes identically, i.e., start at k=1 and stop at k=K-2
	 (trapezoidal rule) */
	en=0.0;
	for(k=1; k<K-1; k++)
	{
		/* trapezoidal rule for u_t^2, use center difference to approximate u_t */
		en+=DR/(2.0*dt*dt)*((u2[k]-u0[k])/2.0*(u2[k]-u0[k])/2.0
							+(u2[k+1]-u0[k+1])/2.0*(u2[k+1]-u0[k+1])/2.0);
		
		/* trapezoidal rule for (u_r-u/r)^2 */
		if(k==1)
		/* special case k=1, exploit vanishing of field at origin for center difference */
			en+=1.0/(2.0*DR)*((u1[2]/2.0-u1[1])*(u1[2]/2.0-u1[1])
							  +((u1[3]-u1[1])/2.0-u1[2]/2.0)*((u1[3]-u1[1])/2.0-u1[2]/2.0));
		else if(k==K-2)
		/* special case last gridpoint, use simple backward difference for u_r, shouldn't do
		 too much harm 'cause nothing interesting happens far away from the origin */
			en+=1.0/(2.0*DR)*(((u1[K-1]-u1[K-3])/2.0-u1[K-2]/((double)(K-2)))
							  *((u1[K-1]-u1[K-3])/2.0-u1[K-2]/((double)(K-2)))
							  +(u1[K-1]-u1[K-2]-u1[K-1]/((double)(K-1)))
							  *(u1[K-1]-u1[K-2]-u1[K-1]/((double)(K-1))));
		else
		/* general case, use central difference for u_r */
			en+=1.0/(2.0*DR)*(((u1[k+1]-u1[k-1])/2.0-u1[k]/((double)k))
							  *((u1[k+1]-u1[k-1])/2.0-u1[k]/((double)k))
							  +((u1[k+2]-u1[k])/2.0-u1[k+1]/((double)(k+1)))
							  *((u1[k+2]-u1[k])/2.0-u1[k+1]/((double)(k+1))));
		
		/* trapezoidal rule for u^2 */
		en+=DR/2.0*(u1[k]*u1[k]+u1[k+1]*u1[k+1]);
		
	}
	
	/* return energy */
	return en/2.0;
}



/* calculates the nonlinear  energy  integral -1/4 u^4/r^2 
 by using trapezoidal rule and second order accurate approximations to derivatives 
 needs one time step in the past and one in the future */
double potentialenergy(double u2[], double u1[], double u0[], long K, double dt)
{
	double en=0.0; /* stores 2 times the energy */
	register long k; /* loop variable */
	
	/* loop over all gridpoints
	 contribution from k=0 vanishes identically, i.e., start at k=1 and stop at k=K-2
	 (trapezoidal rule) */
	en=0.0;
	for(k=1; k<K-1; k++)
	{
				
		/* trapezoidal rule for -1/2 u^4/r^2 */
		en+=1.0/(4.0*DR)*(u1[k]*u1[k]*u1[k]*u1[k]/((double)k*k)
						   +u1[k+1]*u1[k+1]*u1[k+1]*u1[k+1]/((double)(k+1)*(k+1)));
	}
	
	/* return energy */
	return en/2.0;
}



int main(void)
{
	long K; /* number of space points actually calculated, equals N+Kd */
	double dt; /* time mesh  */
	double *u0, *u1, *u2; /* Field at three different time steps */
	double A, B; /* amplitudes */
	double J, K10; /* functionals */
	register long k; /* loop variables */
	FILE *fKMplus, *fKMminus; /* file pointers */
	
	/* Set K */
	K=KD+3;
	
	/* set dt */
	dt=R*DR;
	
	/* allocate memory for arrays */
	if((u0=calloc(K+1, sizeof(double)))==NULL ||
		   (u1=calloc(K+1, sizeof(double)))==NULL ||
		   (u2=calloc(K+1, sizeof(double)))==NULL)
	{
		fputs("Could not allocate memory!\n", stderr);
		exit(EXIT_FAILURE);
	}	
	
	/* allocate memory for soliton */
	if((Q=calloc(K+1, sizeof(double)))==NULL)
	{
		fputs("Could not allocate memory for Q!\n", stderr);
		return EXIT_FAILURE;
	}
	
	/* import soliton */
	puts("Importing soliton...");
	if(!import("solitonQ.dat", Q, K+1))
		fputs("Error importing soliton, assuming Q=0!\n", stderr);
	else
		puts("Ok.");
	
	/* open files */
	if((fKMplus=fopen("JKplus.dat", "w"))==NULL)
	{
		fputs("Could not open file \"KMplus.dat\"!\n", stderr);
		exit(EXIT_FAILURE);
	}
	if((fKMminus=fopen("JKminus.dat", "w"))==NULL)
	{
		fputs("Could not open file \"KMminus.dat\"!\n", stderr);
		exit(EXIT_FAILURE);
	}
	
	
	/* main loop */
	for(B=MINB; B<=MAXB; B+=INCB)
		for(A=MINA; A<=MAXA; A+=INCA)
		{
			
			/* copy initial data to u0 */
			for(k=1; k<K; k++)
				u0[k]=f(k,A);
			
			/* calculate first time step by Taylor expansion */
			for(k=1; k<K; k++)
				u1[k]=u0[k]+g(k,B)*dt+dt*dt/2.0*(D(u0, k, K)+F(u0, k));
			
			/* calculate next time step */
			for(k=1; k<K; k++)
				u2[k]=2.0*u1[k]-u0[k]+dt*dt*(D(u1, k, K)
											 +F(u1, k));	
			
			/* calculate J and K functions */
			J=kineticenergy(u2, u1, u0, K, dt)-potentialenergy(u2, u1, u0, K, dt);
			K10=2.0*kineticenergy(u2, u1, u0, K, dt)-4.0*potentialenergy(u2, u1, u0, K, dt);
			
			/* KMminus */
			if(J<Qenergy) 
			{ 
				if(K10<0)
				{
					printf("Amplitudes A=%f, B=%f: KMminus.\n", A, B);
					fflush(stdout);
					/* write to file */
					fprintf(fKMminus, "%f\t%f\n", A, B);
				}
				/* KMplus */
				else 
				{
					printf("Amplitudes A=%f, B=%f: KMplus.\n", A, B);
					fflush(stdout);
					/* write to file */
					fprintf(fKMplus, "%f\t%f\n", A, B);
				}
			}
			
		}
	printf("Done.\n\nUse\n");
	printf("plot \"KMplus.dat\", \"KMminus.dat\"\n");
	printf("in gnuplot to plot the output!\n\n");
	
	/* close files */
	fclose(fKMplus);
	fclose(fKMminus);
	
	
	/* free memory */
	free(u0); 
	free(u1);
	free(u2);
	
	return EXIT_SUCCESS;
}