/* Calculates time evolution of initial data (u, u_t)=(A*some function, B*some function) for cubic focusing Klein-Gordon equation 
 u_tt-u_rr=-u+u3/r2 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.0

/* stepsize A */
#define INCA 0.008

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

/* maximum amplitude B */    
#define MAXB 6.0

/* stepsize B */
#define INCB 0.008

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

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

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

/* Time steps to calculate */
#define N  3000

/* threshold of L^infty norm for blow up */
#define MAX 100000.0

/* threshold for dispersion */
#define MIN 0.8

/* threshold for dispersion in time, i.e., this is the minimum number of time steps that are calculated until the evolution can be classified as "dispersion" */
#define DISPT 0

/* Maximal number of processes */
#define MAXP 50

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

#define BLOWUP -1
#define INDEC 0
#define DISP 1

#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*exp(-r); 
}

/* initial velocity */
double g(long k, double B)
{
	double r=k*DR; /* compute radial coordinate from gridpoint */
	
	return B*r*exp(-r*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);
}

/* returns the maximum of |u|+|u_t| over specified gridpoints */
double maxuut(double unew[], double uold[], long K, double dt)
{
	register long k; /* loop variable */
	double m=0.0; /* stores maximum */
	
	/* loop over all gridpoints */
	for(k=0; k<K; k++)
	/* store maximum |u|+|u_t|, calculate time derivative by simple backward difference */
		if(fabs(unew[k])+fabs((unew[k]-uold[k])/dt)>m) 
			m=fabs(unew[k])+fabs((unew[k]-uold[k])/dt);
	
	return m;
}

/* calculates time evolution, returns BLOWUP, DISP or INDEC */
int evolution(double A, double B, double u0[], double u1[], double u2[], long K, double dt)
{	
	register long k, n; /* loop variables */
	double *ut; /* auxiliary pointer for cycling */
	double groesse; 
	
	/* 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 further time steps */
	for(n=1; n<=N; n++)
	{	
		/* calculate 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));	
		
		/* compute size */ 
		groesse=maxuut(u2, u1, KD, dt); 
		
		/* check for blow up */
		if(groesse > MAX) return BLOWUP;
		
		/* check for dispersion */
		if(n>DISPT && groesse < MIN) return DISP;
		
		/* cycle arrays */
		ut=u2; /* 2 -> t */
		u2=u0; /* 0 -> 2 */
		u0=u1; /* 1 -> 0 */
		u1=ut; /* t -> 1 */
	}
	
	/* if we ever reach this point then the time evolution was indecisive */
	return INDEC;
}

int main(void)
{
	double A, B; /* amplitudes */
	FILE  *filler; /* file pointers */
	 
	
	/* open files */
	if((filler=fopen("filler.dat", "w"))==NULL)
	{
		fputs("Could not open file \"filler.dat\"!\n", stderr);
		exit(EXIT_FAILURE);
	}
		

	
		/*fill in what's missing with dispersion*/
		for(A=MINA; A<=MAXA; A+=INCA)
		{
			for(B=MINB; B<=MAXB; B+=INCB)
			{  
				if(fabs(A)<=2 && fabs(B)<=2)
					fprintf(filler, "%f\t%f\n", A, B);
				
			}
		}
		
		fflush(filler);
		
	
	
	/* close files and free memory has to be done by every child and the the parent */
	/* close files */
	fclose(filler);
	
	
	
	/* exit is performed by all children and parent */
	exit(EXIT_SUCCESS);
}