#include<stdio.h>
#include<math.h>
#include<sys/time.h>

#define PI 3.141592653589793 //2384626433832795

#define real_type float

#define NX 1024 // Mesh size: X
#define NY 1024 // Mesh size: Y
#define NT 4000 // Number of time steps

// condition: max(p) * tau / sqrt(hx^2+hy^2) < 1

#define X0 0.0 // Border X left
#define X1 4.0 // Border X right

#define Y0 0.0 // Border Y left
#define Y1 4.0 // Border Y right

#define P1 0.1 // pressure value 1
#define P2 0.2 // pressure value 2

const real_type hx = (X1-X0)/(NX-1);
const real_type hy = (Y1-Y0)/(NY-1);
const real_type tau = 0.01;

real_type U1[NY][NX];
real_type U2[NY][NX];

real_type p(int i,int j) { if (j<NY/2) return P1*P1; else return P2*P2; }

double dtime() 
{ struct timeval t; 
  gettimeofday(&t,0);
  return t.tv_sec + t.tv_usec*1e-6;
}

void zero(real_type u[NY][NX])
{ int i,j;
  for (i=0;i<NY;i++)
    for (j=0;j<NX;j++)
      u[i][j]=0;
}

real_type f(int i,int j,real_type t)
{ real_type gamma = 4.0;
  real_type f0 = 1.0;
  real_type t0 = 1.5;
  real_type x1,x2,x3;
  if (i==NY/2 && j==1)
  { x1=2.0*PI*f0*(t-t0);
    x2=exp(-x1*x1/(gamma*gamma));
    x3=sin(x1);
    return x2*x3;
  }
  else return 0;
}

void count(real_type u1[NY][NX],real_type u2[NY][NX],int it)
{ real_type const rhx2=1.0/(2*hx*hx);
  real_type const rhy2=1.0/(2*hy*hy);
  real_type const tau2=tau*tau;
  real_type un;
  int i,j;
  for (i=1;i<NY-2;i++)
    for (j=1;j<NX-2;j++)
    { un=u1[i][j];
      u2[i][j]=2*un-u2[i][j]+tau2*(f(i,j,it*tau)+
             +rhx2*((u1[i+1][j]-un)*(p(i,j-1)+p(i,j))+(u1[i-1][j]-un)*(p(i-1,j-1)+p(i-1,j)))
             +rhy2*((u1[i][j+1]-un)*(p(i-1,j)+p(i,j))+(u1[i][j-1]-un)*(p(i-1,j-1)+p(i,j-1)))
                                  );
    }
}

void save(real_type u[NY][NX],int it)
{ char filename[20];
  int i,j;
  FILE *f;
  sprintf(filename,"%04d",it);
  f=fopen(filename,"wb");
  fwrite(u,sizeof(real_type),NX*NY,f);
  fclose(f);
}

void print(real_type u[NY][NX],int it)
{ char filename[20];
  int i,j;
  FILE *f;
  sprintf(filename,"%04d",it);
  f=fopen(filename,"w");
  for (i=0;i<NY;i++)
    for (j=0;j<NX;j++)
      fprintf(f,"%d\t%d\t%.10e\n",i,j,u[i][j]);
  fclose(f);
}

void print_y(real_type u[NY][NX],int i,int it)
{ char filename[20];
  int j;
  FILE *f;
  sprintf(filename,"%04dy",it);
  f=fopen(filename,"w");
  for (j=0;j<NX;j++)
    fprintf(f,"%e\t%.15e\n",j*hy,u[i][j]);
  fclose(f);
}

real_type max(real_type u[NY][NX])
{ real_type m,x;
  int i,j;
  m=fabs(u[0][0]);
  for (i=0;i<NY;i++)
    for (j=0;j<NX;j++)
    { x=fabs(u[i][j]);
      if (m<x) m=x;
    }
  return m;
}

int main()
{ int it,i,j;
  double t1,t2;

  zero(U1);
  zero(U2);

  t1=dtime();

  for (it=0;it<NT;it+=2)
  {
    count(U1,U2,it  ); //printf("%6d: %le\n",it+1,max(U2));
    if ((it+1)%100==0) { save(U2,it+1); print_y(U2,NX/2,it+1); fprintf(stderr,"%d\n",it+1); }
    count(U2,U1,it+1); //printf("%6d: %le\n",it+2,max(U1));
    if ((it+2)%100==0) { save(U1,it+2); print_y(U1,NX/2,it+2); fprintf(stderr,"%d\n",it+2); }
  }
  t2=dtime();
  printf("Time: %lf\n",t2-t1);
  return 0;
}