/* 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 */

/*  THIS IS THE FORMER 4DAT */

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

/* minimum amplitude A */
#define MINA 0.0

/* maximum amplitude A */
#define MAXA 8.5

/* stepsize A */
#define INCA 0.006

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

/* maximum amplitude B */
#define MAXB 10.0

/* stepsize B */
#define INCB 0.006

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

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

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

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

/* threshold of L^infty norm for blow up */
#define MAX 5000.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 */


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

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


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



/* 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)
{
	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 */
	int result; /* result of time evolution */
	FILE  *fdisp, *findec; /* file pointers */
	pid_t pid; /* stores return value of fork() */
	int status; /* stores status of wait */
	int n;
	int blocks; /* number of chunks into which the B region is split according to MAXP */
	
	/* set blocks */
	
	blocks=ceil((MAXB-MINB)/(MAXP*INCB));
	
	
	/* Set K */
	K=KD+N;
	
	/* set dt */
	dt=R*DR;
	
	/* 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.");
	
	/* 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);
	}	
	
	/* open files */
	if((fdisp=fopen("dispersion.dat", "w"))==NULL)
	{
		fputs("Could not open file \"dispersion_spitz1.dat\"!\n", stderr);
		exit(EXIT_FAILURE);
	}
	if((findec=fopen("indecisive.dat", "w"))==NULL)
	{
		fputs("Could not open file \"indecisive.dat\"!\n", stderr);
		exit(EXIT_FAILURE);
	}
	
	/* set pid to 1 initially, parent process is identified by pid!=0 */
	pid=1;
	
	/* main loop */
	for(n=0; n<blocks; n+=1)
	{
		if(pid!=0)
		{
			for(B=MINB+n*MAXP*INCB; B<MINB+(n+1)*MAXP*INCB; B+=INCB)
			{
		
				/* create process */
				if(pid!=0) 
				{
					pid=fork();
			
					/* child process has pid=0 or pid=-1 if fork failed */
					if(pid==0 || pid==-1)
					{				
						for(A=MINA; A<=MAXA; A+=INCA)
						{
							if((A>4 || fabs(B+A)>5) && (A+B>-11))
							{
							/* calculate time evolution */
							result=evolution(A, B, u0, u1, u2, K, dt);
					
					
							/* dispersion */
							if(result==DISP) 
							{
								printf("Amplitudes A=%f, B=%f: dispersion.\n", A, B);
								fflush(stdout);
								/* write to file */
								fprintf(fdisp, "%f\t%f\n", A, B);
								fflush(fdisp);
							}
							/* indecisive */
							else if(result==INDEC)
							{
								printf("Amplitudes A=%f, B=%f: indecisive.\n", A, B);
								fflush(stdout);
								/* write to file */
								fprintf(findec, "%f\t%f\n", A, B);
								fflush(findec);
							}
						}
						}
					}
				}
			}
		}
		while(wait(&status)!=-1);	
		
	}
		
	
	
	/* pid contains the pid of the last fork, only parent has != 0 */
	if(pid!=0)
	{
		/* this is only executed in the parent process */
		
		/* wait for all children to finish */
		while(wait(&status)!=-1);
		
		/*fill in what's missing with dispersion */
		for(A=MINA; A<=MAXA; A+=INCA)
		{
			for(B=MINB; B<=MAXB; B+=INCB)
			{  
				if( A<=4 && fabs(B+A)<=5 )
				{
					fprintf(fdisp, "%f\t%f\n", A, B);
				fflush(fdisp);
				}
			}
		}
		
		printf("Done.\n\nUse\n");
		printf("plot   \"dispersion.dat\", \"indecisive.dat\"\n");
		printf("in gnuplot to plot the output!\n\n");
	}
	
	/* close files and free memory has to be done by every child and the the parent */
	/* close files */
	fclose(findec);
	fclose(fdisp);
	
	
	/* free memory */
	free(Q);
	free(u0); 
	free(u1);
	free(u2);
	
	/* exit is performed by all children and parent */
	exit(EXIT_SUCCESS);
}