#include <iostream>
#include <fstream>
#include <map>
#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;
using namespace rbf;

class wave_op : public linear_diff_op2{
public:
  wave_op(double h_in) : h(h_in) {} ;
  double at(const realfunc &f, const point &p) const;
  double at(const radial_basis_function &RBF, const point &p) const;
private:
  Laplacian L;
  double h;
};

double wave_op::at(const realfunc &f, const point &p) const{
  return f(p) - h*h*L(f,p);
}

double wave_op::at(const radial_basis_function &RBF, const point &p) const{
  return RBF(p) - h*h*L(RBF,p);
}


using boost::shared_ptr;

std::map<point,double> calc_f(const std::map<point,double>& u0,
			      const std::map<point,double>& u1,
			      const matrix& intr);

template<typename RBF>
void save_step(const interpolator<RBF>& u, 
	       const matrix& intr,
	       const matrix& bdry,
	       size_t n);

template<typename RBF>
map<point, double> give_u1(const interpolator<RBF>& u,
			   const matrix& intr);

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

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

  try{
    using namespace std;

    using linalg::vector;
    shared_ptr<domain> Omega(new domain("data/circ_intr.matrix",
					"data/circ_bdry.matrix",
					"data/circ_nrml.matrix")
			     );
    //Interior conditions, init them
    map<point,double> u_init = utils::read_pd_map("data/wave_init.map");

    //Boundary conditions, Dirichlet, set to zero (clamped membrane)
    zero_func g;
    
    //timestep
    double h = 0.01;
    
    shared_ptr<wave_op> W(new wave_op(h) );
    shared_ptr<dirichlet_op> D(new dirichlet_op);
    shared_ptr<linear_BVP2> the_bvp(new linear_BVP2(Omega,W,D,u_init,g) );

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

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

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

}

std::map<point,double> calc_f(const std::map<point,double>& u0,
			      const std::map<point,double>& u1,
			      const matrix& intr)
{
  std::map<point, double> f;
  slice s(1,2);
  for(size_t i = 1; i <= intr.rows(); i++){
    f[intr(i,s)] = 2*u1.at(intr(i,s)) - u0.at(intr(i,s));
  }
  return f;
}

template <typename RBF>
void save_step(const interpolator<RBF>& u,
	       const matrix& intr,
	       const matrix& bdry,
	       size_t n)
{
  using namespace std;

  string filename = "results/u";
  if(n < 10000)
    filename += "0";
  if(n < 1000)
    filename += "0";
  if(n < 100)
    filename += "0";
  if(n < 10)
    filename += "0";
  stringstream ss;
  string n_string;
  ss << n;
  ss >> n_string;
  //filename += n_string + ".vector";
  filename += n_string + ".matrix";
  ofstream ofs(filename.c_str());

//   slice s(1,2);
//   for(size_t i = 1; i <= intr.rows(); i++)
//     ofs << u.at(intr(i,s)) << endl;
//   for(size_t i = 1; i <= bdry.rows(); i++)
//     ofs << u.at(bdry(i,s)) << endl;
  intr.rows();
  bdry.rows();
  size_t rows = 20, cols = 30;
  matrix result(rows,cols);
  double r = 0;
  for(size_t i = 1; i <= rows; i++){
    double th = 0;
    for(size_t j = 1; j <= cols; j++){
      point p(2);
      p(1) = r*cos(th);
      p(2) = r*sin(th);
      result(i,j) = u.at(p);
      th += 2*M_PI/(double(cols)-1);
    }
    r += 1.0/(double(rows)-1);
  }
  ofs << result;
}

template<typename RBF>
map<point, double> give_u1(const interpolator<RBF>& u,
			   const matrix& intr)
{
  map<point, double> u1;
  slice s(1,2);
  for(size_t i = 1; i <= intr.rows(); i++)
    u1[intr(i,s)] = u.at(intr(i,s));
  return u1;
}