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

#define MAXDIM 1010
#define ANIME_STEP 1

void set_xt (int nx, int nt, double x[], double t[]);
void initial_condition(int nx, int nt, double delta_t, double alpha,
		       double x[], double u[][MAXDIM]);
double phi(double x);
double psi(double x);
void  write_cal_result(char *plot_data_file_name,
		       int nx, int ntime,
		       double x[],double u[][MAXDIM]);
void plot_2d (FILE * gp, char *data_file, int max_nt);

/*=====================================================================*/
/*   main                                                              */
/*=====================================================================*/
int main(void){
  double u[MAXDIM][MAXDIM];
  double x[MAXDIM], t[MAXDIM];
  double xmin, xmax, tmin, tmax;
  double delta_x, delta_t, alpha;
  int i, j;
  int nx, nt;
  char *plot_data_file_name;
  FILE *gp;

  plot_data_file_name = "./data/result";

  nt = 1000;
  nx = 1000;

  xmin=0.0;
  xmax=1.0;

  tmin=0.0;
  tmax=1.0;

  delta_x = (xmax-xmin)/nx;
  delta_t = (tmax-tmin)/nt;
  alpha = (delta_t/delta_x)*(delta_t/delta_x);

  set_xt(nx, nt, x, t);
  initial_condition(nx, nt, delta_t, alpha, x, u);

  for (i = 2; i <= nt; i++){
    for (j = 1; j < nx ; j++){
      u[j][i]=2.0*u[j][i-1]-u[j][i-2]
	+alpha*(u[j+1][i-1]-2.0*u[j][i-1]+u[j-1][i-1]);
    }
  }

  for(i=0; i<=nt; i+=ANIME_STEP){
    write_cal_result(plot_data_file_name, nx, i, x, u);
  }

  gp = popen("gnuplot -persist", "w");
  plot_2d(gp, plot_data_file_name,nt);
  pclose (gp);

  return 0;
}

/*=====================================================================*/
/*   set x-coordinate and time                                         */
/*=====================================================================*/
void set_xt(int nx, int nt, double x[], double t[]){

  int i;   	

  for (i = 0; i <= nx; i++){
	  x[i] = (double) i / nx;
  }

  for (i = 0; i <= nt; i++){
	  t[i] = (double) i / nt;
  }

}


/*=====================================================================*/
/*   set initial condition                                             */
/*=====================================================================*/
void initial_condition(int nx, int nt, double delta_t, double alpha,
		       double x[], double u[][MAXDIM]){

  int i;

  /* --------- boundary condition at x=0 and x=1  ------- */

  for(i = 0; i <= nt ; i++){
    u[0][i]=0.0;
    u[nx][i]=0.0;
  }

  /* --------- condition at t=0 ------------ */
  for (i = 1; i < nx ; i++){
    u[i][0]=phi(x[i]);
  }


  /* --------- condition at t=delta t ------------ */

  for (i = 1; i < nx ; i++){
    u[i][1]=u[i][0]+psi(x[i])*delta_t+
      alpha/2.0*(u[i+1][0]-2*u[i][0]+u[i-1][0]);
  }

}

/*=========================================================*/
/*  function phi(x)                                        */
/*=========================================================*/
double phi(double x){

  double a1, sigma1, x1, wave1;
  double a2, sigma2, x2, wave2;

  a1=0.5;
  sigma1 = 0.02;
  x1=0.8;

  wave1=a1*exp(-(x-x1)*(x-x1)/(sigma1*sigma1));

  a2=20.0;
  sigma2 = 0.02;
  x2=0.4;
  wave2=a2*(x-x2)*exp(-(x-x2)*(x-x2)/(sigma2*sigma2));


  return(wave1+wave2);

}


/*=========================================================*/
/*  function psi(x)                                        */
/*=========================================================*/
double psi(double x){

  double a1, sigma1, x1, wave1;
  double a2, sigma2, x2, wave2;

  a1=0.5;
  sigma1 = 0.02;
  x1=0.8;

  wave1=a1*exp(-(x-x1)*(x-x1)/(sigma1*sigma1))*
    (-2.0*(x-x1)/(sigma1*sigma1));


  a2=20.0;
  sigma2 = 0.02;
  x2=0.4;

  wave2=a2*exp(-(x-x2)*(x-x2)/(sigma2*sigma2))+
    a2*(x-x2)*exp(-(x-x2)*(x-x2)/(sigma2*sigma2))*
    (-2.0*(x-x2)/(sigma2*sigma2));

  return(wave1-wave2);

}


/*=========================================================*/
/*  write calculation resuts                               */
/*=========================================================*/
void  write_cal_result(char *plot_data_file_name,
		       int nx, int ntime,
		       double x[],double u[][MAXDIM]){

  int i;
  char num[20], file[50];
  FILE *fp;

  sprintf(num,"_%d.txt",ntime);
  sprintf(file, "%s%s",plot_data_file_name, num);

  printf("file name = %s\n", file);

  fp = fopen (file, "w");

  for (i = 0; i <= nx; i+=ANIME_STEP){
      fprintf (fp, "%f\t%f\t\n", x[i], u[i][ntime]);
  }

  fclose (fp);

}

/*==========================================================*/
/*   2D plot                                                */
/*==========================================================*/
void plot_2d(FILE * gp, char *data_file, int max_nt)
{

  int i,k;
  char num[20], file[50];

  fprintf (gp, "reset\n");
  fprintf (gp, "set terminal x11\n");
  fprintf (gp, "set nokey\n");
  fprintf (gp, "set xrange[0:1]\n");
  fprintf (gp, "set yrange[-1:1]\n");

  for(k=0; k<1; k++){
    for(i=0; i<=max_nt ; i+=ANIME_STEP){

      sprintf(num,"_%d.txt",i);
      sprintf(file, "%s%s",data_file, num);

      fprintf (gp, "plot \"%s\" using 1:2 with lines\n", file);
    }
  }


}