/* Numerical integration of the transformed cubic focusing Klein-Gordon equation 
u_tt-u_rr=-u+u^3/r^2 in radial 3d using centered time and centered space discretization */

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

/* Amplitude of initial position */
#define A 1.0

/* Amplitude of initial velocity */
#define B 0.0

/* Grid mesh dr */
#define DR 0.0002

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

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

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

/* Height of output window */
#define MAXY 3.0

/* delay factor for movie, increase if movie runs too fast, 0.0 means no delay at all */
#define DELAY 0.0

/* time steps to write out for movie, 1=write every step, 2=write every second step, etc.
 increase if movie runs too slow */
#define SKIPT 38

/* gridpoints to write out for movie, 1=write every point, 2=write every second point, etc.
increase if movie doesn't run smoothly */
#define SKIPR 20

#define MAX 10000

#define DISP .4

#define PI 3.1415926


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

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

/****************** initial data *****************/



/* initial position */
double f(double r)
{

    if (r>.5 && r<1) 
	{
		return 2*A*sin(2*PI*(r-.75));
	}
	
	else if (r<=0.5) 
	{
		return -16*A*r*r*r*(1+(0.5-r)*(0.5-r));
	}
	
	else if (r>=1)
	{
		return 2*A*(1+(r-1)*(r-1))*exp(-r*r+1);
	}


}

/* initial velocity */
double g(double r)
{

	
	if (r>.25 && r<.5) 
	{
		return B*sin(2*PI*(2*r-.75))/2;
	}
	
	else if (r<=0.25) 
	{
		return -32*B*r*r*r*(1+(0.25-r)*(0.25-r));
	}
	
	else if (r>=.5)
	{
		return B*(1+(r-.5)*(r-.5))*exp(-r*r+.25)/2;
	}
	
}

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

/* 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 energy 1/2 integral u_t^2+(u_r-u/r)^2+u^2-1/2 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 energy(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]);
		
		/* 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;
}
	
	
/* 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;
}

	
int main(void)
{
	long K; /* number of space points actually calculated, equals N+Kd */
	double dt, groesse; /* time mesh  */
	double *u0, *u1, *u2; /* Field at three different time steps */
	double *ut; /* auxiliary pointer for cycling */
	
	register long k, n; /* loop variables */
	long nsave; /* counter for saved time steps */
	FILE *fout, *fanim, *fenergy;
	
	/* Set K */
	K=KD+N;
	
	/* 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 output file */
	if((fout=fopen("data.dat", "w"))==NULL)
	{
		fputs("Could not open output file data.dat!\n", stderr);
		exit(EXIT_FAILURE);
	}
	
	/* open anim file */
	if((fanim=fopen("anim", "w"))==NULL)
	{
		fputs("Could not open animation file anim!\n", stderr);
		exit(EXIT_FAILURE);
	}	
	
	/* open energy file */
	if((fenergy=fopen("energy.dat", "w"))==NULL)
	{
		fputs("Could not open energy file energy.dat!\n", stderr);
		exit(EXIT_FAILURE);
	}	
	
	/* set dt */
	dt=R*DR;
	
	/* copy initial data to u0 */
	for(k=1; k<K; k++)
		u0[k]=f(k*DR);
	
	/* calculate first time step by Taylor expansion */
	for(k=1; k<K; k++)
		u1[k]=u0[k]+g(k*DR)*dt+dt*dt/2.0*(D(u0, k, K)+F(u0, k));
	
	/* calculate further time steps */
	nsave=0; /* set counter for saved time steps */
	for(n=1; n<=N; n++)
	{
		printf("Calculating time step number %ld\r", n);
		fflush(stdout); /* force output of text */
		
		/* 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) 
		{
		fputs("\n Finite time blowup!\n", stderr);
		printf("\nDone!\nUse \"gnuplot anim\" to watch a movie!\n\n");
		return EXIT_SUCCESS;
		}
		
		else if (groesse < DISP)  
		{
		fputs("\n Dispersion!\n", stderr);
		printf("\nDone!\nUse \"gnuplot anim\" to watch a movie!\n\n");
		return EXIT_SUCCESS;
		}
		
		/* write data to file */
		if(n%SKIPT==0)
		{
			for(k=0; k<KD; k++)
			{
				if(k%SKIPR==0)
					fprintf(fout, "%f\t%f\n", k*DR, u2[k]);
			}
			fprintf(fout, "\n\n");
			
			/* write anim command */
			fprintf(fanim, "plot [0:%f] [-%f:%f] \"data.dat\" index %ld\n",
		       		KD*DR, MAXY, MAXY, nsave);
			fprintf(fanim, "pause %f\n", DELAY);
			
			/* write energy for last time step */
			fprintf(fenergy, "%f\t%f\n", (n-1)*dt, energy(u2, u1, u0, K, dt));
			
			/* increase counter */
			nsave++; 
		}
		
		/* cycle arrays */
		ut=u2; /* 2 -> t */
		u2=u0; /* 0 -> 2 */
		u0=u1; /* 1 -> 0 */
		u1=ut; /* t -> 1 */
	}

	/* close files */
	fclose(fout);
	fclose(fanim);
	fclose(fenergy);
	
	/* free memory */
	free(u0); 
	free(u1);
	free(u2);
	
	printf("\nDone!\nUse \"gnuplot anim\" to watch a movie!\n\n");
	
	return EXIT_SUCCESS;
}