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

// Initial conditions: the velocity is zero

#define sqr(x) ((x)*(x))


//#define L 160.0
#define L 40.0
//#define N 1600
#define N 400
//#define T 800.0
#define T 400.0
//#define M 80000
#define M 40000
#define m0 800
#define a ((L)/(N))
#define h ((T)/(M))
#define v 0.25
#define C sqr((h)*(v)/(a)) // should be <1
#define pi acos(-1.0L)
#define z0 (0.5*(L))
#define dz ((z0)/5.0)
#define var sqr(0.02*(L))
#define dv 1.0



double phi_new[N+1], phi_curr[N+1], phi_old[N+1];

void initc(){// initial conditions
  int n;
  double z;
  double gauss_G, gauss_D; 
  for(n=1;n<N;n++){
    z = n*a;
    //    phi_curr[n] = sin(2*pi*z/L); // one mode
    gauss_G = exp(-sqr(z-z0-dz)/(2.0*var*dv))/sqrt(2.0*pi*var*dv); //gaussian_G
    gauss_D = exp(-sqr(z-z0+dz)/(2.0*var))/sqrt(2.0*pi*var); //gaussian_D

    phi_curr[n] = 0.5*(gauss_G + gauss_D);
  }
}

void bc(){// boundary conditions
  phi_old[0] = 0.0; phi_old[N] = 0.0;
  phi_curr[0] = 0.0; phi_curr[N] = 0.0;
  phi_new[0] = 0.0; phi_new[N] = 0.0;
}

void prempas(){// t=0->t=1
  int n;

  for(n=1;n<N;n++){
    phi_new[n] = phi_curr[n] + (C/2.0)*(phi_curr[n+1]-2.0*phi_curr[n]+phi_curr[n-1]);
  }

}

void update(){
  int n;

  for(n=1;n<N;n++){
    phi_new[n] = 2.0*phi_curr[n]+C*(phi_curr[n+1]-2.0*phi_curr[n]+phi_curr[n-1])-phi_old[n];
  }

}

void replace(){
  int n;

  for(n=1;n<N;n++){
    phi_old[n] = phi_curr[n];
  }

  for(n=1;n<N;n++){
    phi_curr[n] = phi_new[n];
  }
}

void sortie(double t){
  int n;
    for(n=0;n<=N;n++){
      printf("%g\t%g\t%g\n", t, n*a, phi_curr[n]);
    }
}

main(){
  int m;
  double t;

  //  printf("%g\n", C);

    initc();
  bc();

  prempas(); replace();

  for(m=2;m<M;m++){
    t = m*h;
    update(); 

    if((m%m0)==0){
      sortie(t);}

    replace();

  }


}