#include <iostream>
#include <iomanip>
#include <fstream>
#include <map>
#include <set>
#include <sstream>
 
#include <boost/shared_ptr.hpp>

#include "include/linalg.hpp"
#include "include/rbf.hpp"
#include "include/bvp.hpp"
#include "include/error.hpp"
#include "include/utils.hpp"
#include "include/diff_op.hpp"
#include "include/interpolator.hpp"
#include "include/func.hpp"
#include "include/vtkplot.hpp"

using namespace kwantix;

//Define the RBF type here.
typedef kwantix::conical RBF_TYPE;

const double pi = 3.141592653589793238462643383279502;

class wave_op : public linear_diff_op2{
public:
  wave_op(double dt_in): dt2(dt_in*dt_in){};
  double at(const realfunc& f, const point& p) const
  {
    return dt2*L(f,p) - f(p);
  }
private:
  Laplacian L;
  double dt2;
};

class zero_func : public realfunc{
public:
  double at(const point&) const{return 0;};
};

void save_iteration(const interpolator<RBF_TYPE>& u,
                    size_t n);

int main()
{  
  gsl_set_error_handler(&kwantix::errorHandler);

  try
  {
    using namespace std;
    using boost::shared_ptr;
   
    map<point, double> u_init = kwantix::read_pd_map("data/wave_init.map");
      
    shared_ptr<domain> Omega(new domain("data/circ_intr.matrix",
                                        "data/circ_bdry.matrix",
                                        "data/circ_nrml.matrix"));    
    zero_func g;

    //timestep
    double dt=1e-2;
    
    shared_ptr<wave_op> W(new wave_op(dt) );
    shared_ptr<dirichlet_op> D(new dirichlet_op);
    shared_ptr<linear_BVP2> the_bvp(new linear_BVP2(Omega, W, D,
                                                    u_init, g));

    kwantix::conical::set_n(5);     
    interpolator<RBF_TYPE> u(the_bvp);
    u.precompute_ev();
    interpolator<RBF_TYPE> u0 = u, u1 = u;
    
    u = 5*u;

    vtkplot plotter(u,-0.1, 0.1);
    plotter.set_view_direction(0,1,1);
    plotter.begin_interaction();

    //Main loop
    size_t maxiter = 500;
    for(size_t n = 1; n <= maxiter; n++)
    {
      plotter.plot();
      u0 = u1;
      u1 = u;
      cout << "Iteration: " << n << endl;
      u.set_f( (u0 - 2*u1).at() );
      plotter.update_values(u.at());
    }

  }
  catch(error exc)
  {
    kwantix::show_exception(exc);
    return 1;
  }

}