#include<math.h>
#include<iostream>
#include<fstream>
#include<string>

using namespace std;

class String
{
  private:
    String() {}

  public:
    double l;
    double dx;
    int nx;
    double dt;
    double c2;
    double v2;
    double tmax;
    double *u0;
    double *u1;
    double *u2;

    String(double _l, double k, double rho, double _dx): l(_l), dx(_dx)
      { u0 = u1 = u2 = NULL; v2 = k/rho; init(); }

    ~String() { delete u0; delete u1; delete u2; u0 = u1 = u2 = NULL; }

    void init()
      {
        // delete all arrays
        if(u0) delete u0;
        if(u1) delete u1;
        if(u2) delete u2;

        // create new arrays
        nx = int(l/dx)+1;
        u0 = new double[nx];
        u1 = new double[nx]; 
        u2 = new double[nx]; 

        // initial condition u(x,t=0)
        for(int ix = 0; ix < nx; ix++) {
          double x = dx*ix;
//          u0[ix] = x <= 0.8*l ? 1.25*x/l : 5*(1-x/l);
//          u0[ix] = 0.001*sin(2*3.141592654*x/2/l);

          if(x <= 0.1)
            u0[ix] = 0.;
          if(x > 0.1 && x <= 0.2) 
            u0[ix] = 10*x-1;
          if(x > 0.2 && x <= 0.3) 
            u0[ix] = -10*x+3;
          if(x > 0.3 && x <= 0.7) 
            u0[ix] = 0;
          if(x > 0.7 && x <= 0.8) 
            u0[ix] = 10*x-7;
          if(x > 0.8 && x <= 0.9) 
            u0[ix] = -10*x+9;
          if(x > 0.9) u0[ix] = 0.;

        }
      }

    void solve(double _dt, double _tmax)
      {
        dt = _dt;
        tmax = _tmax;
        c2 = dx*dx/dt/dt;

        int nt = int(tmax/dt)+1; 

        // output initial configuration
        ofstream fout;
        fout.open("s.dat",std::ios::out);
        for(int ix = 0; ix < nx; ix++){ 
          fout << ix*dx << " " << u0[ix] << " " << 0 << endl;
        }
        fout << endl;
//        fout.close();
        /////////////////////////////////////////////

        u1[0] = u1[nx-1] = 0;
        for(int ix = 1; ix < nx-1; ix++) {
          u1[ix] = u0[ix] + 0.5*v2/c2*(u0[ix+1]+u0[ix-1]-2*u0[ix]);
        }
        int file_index = 1;
        // Loop in time
        for(int it = 2; it < nt; it++) { 
          // calculate the new u(x,t) from u0 and u1
          u2[0] = u2[nx-1] = 0;
          for(int ix = 1; ix < nx-1; ix++) {
            u2[ix] = 2*u1[ix] - u0[ix] + v2/c2*(u1[ix+1]+u1[ix-1]-2*u1[ix]);
//cout << ix << " " << u0[ix] << endl;
          }
          ////////////////////////////////////////////////////////////
          // output
          if(it%20000 == 0){
            char filename[200];
            sprintf(filename,"s%d.dat",file_index++);
//            fout.open(filename,std::ios::out);
            for(int ix = 0; ix < nx; ix++) 
              fout << ix*dx << " " << u2[ix] << " " << it*dt << endl;
            fout << endl;
//            fout.close();
          }
          ////////////////////////////////////////////////////////////
          // move one time slice
          for(int ix = 0; ix < nx; ix++) { u0[ix] = u1[ix]; u1[ix] = u2[ix]; }
        }
        fout.close();
      }
};

int
main()
{
  double l, dx, dt, k, rho, tmax;

  cout << "Lenght (l):" << endl;
  cin >> l;
  cout << "Tension (k):" << endl;
  cin >> k;
  cout << "Linear density (rho):" << endl;
  cin >> rho;
  cout << "Dx: " << endl;
  cin >> dx;
  cout << "Dt: " << endl;
  cin >> dt;
  cout << "Tmax: " << endl;
  cin >> tmax;
  cout << "V2=" << k/rho << endl;
  cout << "C2=" << dx*dx/dt/dt << endl;

  String s(l, k, rho, dx);
  s.solve(dt, tmax);
  return 0;
}