/*

ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc

c  twodsol.c:  Solves the sine-Gordon equation for a 2D soliton        c

c								       c

c  taken from: "Projects in Computational Physics" by Landau and Paez  c 

c	       copyrighted by John Wiley and Sons, New York            c 

c              code copyrighted by RH Landau                           c

c  supported by: US National Science Foundation, Northwest Alliance    c

c                for Computational Science and Engineering (NACSE),    c

c                US Department of Energy 	                       c

c								       c

c    Christiansen & Lomdahl, Physics 2D (1981) 482-494                 c

c    U_xx + U_yy-U_tt=j(x,y)sin (U)                                    c

c        -x_0 <x< x_0,    -y_0 <y< y_0, t >=0,                         c

c    i. c. U(x,y,0)=4 arctan (exp(3-sqrt(x^2+y^2)),                    c

c    i. c. d U(x,y,0)/dt =0,  j(x,y)=1,  x_0=y_0=7                     c

cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc

*/

#include<stdio.h>

#include<math.h>

 

#define D 1200



main()

{ 



double u[D][D][3];

int nint;

void initial(double u[][D][3]);

void solution(double u[][D][3], int nint);



/*    input a positive integer which is proportional to the time

      you want to see the position of the wave packet.		*/

/*printf(" Enter a positive integer from 1(initial time)\n");*/

/*printf("to 1800 to get wave packet position at that time:\n");*/

scanf("%d", &nint);



/* initializes the constant values and the wave */

initial(u);



/* time dependent equation is solved */

solution(u,nint);



}



void initial(double u[][D][3])

{

double dx,dy,dt,xx,yy,dts,time, tmp;

int i,j;

dx=14./200.;

dy=dx;

dt=dx/sqrt(2.);

dts=(dt/dx)*(dt/dx);

yy=-42.;

time=0.;

for(i=0;i<=D-1;i++){

	xx=-42.;

	for(j=0;j<=D-1;j++){

		tmp=3.-sqrt(xx*xx+yy*yy);

		u[i][j][0]=4.*atan(tmp);

		xx=xx+dx;

		}

	yy=yy+dy;

	}

}

	

void solution(double u[][D][3], int nint)	

{

double dx,dy,dt,time,a2,zz,dts,a1, tmp;

int l,m,mm,k,j,i;

dx=14./200.;

dy=dx;

dt=dx/sqrt(2.);

time=0.;

time=time+dt;

dts=(dt/dx)*(dt/dx);

tmp=0.;

for(m=1;m<=D-2;m++){

	for(l=1;l<=D-2;l++){

		a2=u[m+1][l][0]+u[m-1][l][0]+u[m][l+1][0]+u[m][l-1][0];

		tmp=.25*a2;

		u[m][l][1]=.5*(dts*a2-dt*dt*sin(tmp));

		

		}

	}

	

/* The borders in the second iteration */

for(mm=1;mm<=D-2;mm++){

	u[mm][0][1]=u[mm][1][1];

	u[mm][D-1][1]=u[mm][D-2][1];

	u[0][mm][1]=u[1][mm][1];

	u[D-1][mm][1]=u[D-2][mm][1];

	}

	

/* The still undefined terms */

u[0][0][1]=u[1][0][1];

u[D-1][0][1]=u[D-2][0][1];

u[0][D-1][1]=u[1][D-1][1];

u[D-1][D-1][1]=u[D-2][D-1][1];

/*Third and following iterations k=1,2... */

tmp=0.;	

for(k=0;k<=nint;k++){

	for(m=1;m<=D-2;m++){

		for(l=1;l<=D-2;l++){

			a1=u[m+1][l][1]+u[m-1][l][1]+u[m][l+1][1]+u[m][l-1][1];

			tmp=.25*a1;

			u[m][l][2]=-u[m][l][0]+dts*a1-dt*dt*sin(tmp);

			u[m][0][2]=u[m][1][2];

			u[m][D-1][2]=u[m][D-2][2];

			}

		}

		for(mm=1;mm<=D-2;mm++){

			u[mm][0][2]=u[mm][1][2];

			u[mm][D-1][2]=u[mm][D-2][2];

			u[0][mm][2]=u[1][mm][2];

			u[D-1][mm][2]=u[D-2][mm][2];

			}

		u[0][0][2]=u[1][0][2];

		u[D-1][0][2]=u[D-2][0][2];

		u[0][D-1][2]=u[1][D-1][2];

		u[D-1][D-1][2]=u[D-2][D-1][2];

		

	/* New iterations are now old, reuse, recycle, recover */

		for(l=0;l<=D-1;l++){

			for(m=0;m<=D-1;m++){

				u[l][m][0]=u[l][m][1];

				u[l][m][1]=u[l][m][2];

				}

			}

		if(k==nint) {

			for(i=0;i<=D-1;i=i+5){

				for(j=0;j<=D-1;j=j+5){

					printf("%e\n",sin(u[i][j][2]/2.)

					);

					}

				printf("\n");

				}

			}

		time=time+dt;

		}

}