#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"

using namespace linalg;
using namespace bvp;

#define RBF_TYPE rbf::conical

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(&error_handling::errorHandler);

  try{
    using namespace std;
    using boost::shared_ptr;
   
    map<point, double> u_init = utils::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));

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

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

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

}

void 
save_iteration(const interpolator<RBF_TYPE>& u,
               size_t n)
{
  using namespace std;
  stringstream ss;
  string filename = "results/u";
  ss << filename << setw(5) << setfill('0') << n << ".map";
  ss >> filename;
  ofstream ofs(filename.c_str());
  ofs << u;
  ofs.close();
}