#include <stdio.h>
#include <unistd.h>

#define NT 800
#define NX 200
#define SLEEP_TIME 20000

void set_xt (int nx, int nt, double x[], double t[]);
void initial_condition(int nx, int nt, double alpha, double u[][NT+1]);
/*=====================================================================*/
/*   main                                                              */
/*=====================================================================*/
int main(){
  double u[NX+1][NT+1];
  double x[NX+1], t[NT+1];
  double xmin, xmax, tmin, tmax;
  double delta_x, delta_t, alpha;
  int i, j;
  FILE *gp;

  xmin=0.0;
  xmax=1.0;

  tmin=0.0;
  tmax=2.7;

  delta_x = (xmax-xmin)/NX;
  delta_t = (tmax-tmin)/NT;
  alpha = (delta_t/delta_x)*(delta_t/delta_x);

  printf("alpha = %f\n",alpha);

  if(alpha >= 1.0){
    printf("alpha must be less than 1. Now alpha = %f\n",alpha);
    return 1;
  }

  set_xt(NX, NT, x, t);
  initial_condition(NX, NT, alpha, u);

  gp = popen("gnuplot -persist","w");
  fprintf(gp, "set xrange [-0.1:1.1]\n");
  fprintf(gp, "set yrange [-0.6:0.6]\n");
  
  for(j=0; j<=1; j++){
    printf("t=%f\n",t[j]);
    fprintf(gp, "plot '-' with lines linetype 1\n");
    for(i=0; i<=NX; i++){
      fprintf(gp,"%f\t%f\n", x[i], u[i][j]);
    }
    fprintf(gp, "e\n");
    usleep(SLEEP_TIME);
  }


  for(j=2; j<=NT; j++){
    printf("%d\tt=%f\n",j,t[j]);
    for(i=1; i<NX ; i++){
      u[i][j]=2.0*u[i][j-1]-u[i][j-2]
	+alpha*(u[i+1][j-1]-2.0*u[i][j-1]+u[i-1][j-1]);
    }
    fprintf(gp, "plot '-' with lines linetype 1\n");
    for(i=0; i<=NX; i++){
      fprintf(gp,"%f\t%f\n", x[i], u[i][j]);
    }
    fprintf(gp,"e\n");
    usleep(SLEEP_TIME);
  }

  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 alpha, double u[][NT+1])
{
  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++){
    if(i <= nx/2){
      u[i][0]=(double)i/nx;
    }else{
      u[i][0]=1.0-(double)i/nx;
    }
  }

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

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

}