#ifndef KERNEL_HPP
#define KERNEL_HPP
#include <cstdio>
#include <cstring>
#include <fstream>
#include <stack>
#include "physics.hpp"
#include "geometry.hpp"
#include "time.hpp"
#include "input_data.hpp"
#include "solution.hpp"
#include "piccard.hpp"

using namespace std;

class Kernel:public InputData{
private:
  static constexpr double tolerence_Piccard=1e-6;
  static constexpr double oscilation_Piccard=1e-4;
  static constexpr size_t max_iterations_Piccard=50;
  static constexpr size_t max_time_subdivisions=32;

  static constexpr size_t max_iterations_Richards=50;
  static constexpr double tolerence_Richards=1e-10;

  //For which colum we have to compute Richards
  bool* indice_x_Richards;
  int scheme_error;



  void initSolutions();

  void spaceSolution();

  double bottomFluxRichards(size_t ix,const Solution& S_k);
  Solution* Piccard(Solution* Sin);
  bool Richards1dEvolutiveTime(size_t ix,const Solution& S_in,const Solution& S_k,Solution& S_kp);
  bool allVerticalRichards(Solution& S_in,Solution& S_k,Solution& S_kp);
  void compute_l(Solution& S,double& error);
  void compute_Pl(Solution& S,const double* h,double** P);
  void compute_mass(Solution& S);
  void compute_div_w(size_t ix,const Solution& S);
  bool horizontalProblem(Solution& S,double& error);
  bool overlandProblem(const Solution& S_in,Solution& S_out,double& error);
  void saturationInterface(const Solution& S_in,Solution& S_out);
  bool Richards1d(size_t ix,const Solution& S_0,const Solution& S_s,Solution& S_sp,double dt,double flux_bot);
  bool weightedRichards1d(size_t ix,double* P0,double* Ps,double* Psp,double dt,double flux_bot,double r,double n_P0);
  void solveSystemRichards(size_t ix,double* P0,double* Pin,double* Pout,double flux_bot,double dt,size_t count);
//  double computeError(const Solution& S_in,const Solution& S_out,size_t sol_temp);
public:
  size_t step;
  //Number of sub time intervals of [step,step+1]
  size_t m;
  Solution* solution;
  double** div_w;
  double** P_temp;

  double** t_sub_A;
  double** t_sub_B;
  double** t_sub_C;
  double** t_diag_A;
  double** t_diag_B;
  double** t_diag_C;
  double** t_sup_A;
  double** t_sup_B;
  double** t_sup_C;
  double** t_F;
  double** t_G;
  double** t_S;
  double** t_H;
  double** t_R;
  double** t_I;
  double** t_BB;
  size_t nZ(size_t ix);
  double Z(size_t ix,size_t iz);

  Kernel(string filename);
  void init(fstream& file);
  bool next();

};

inline size_t
Kernel::nZ(size_t ix){
  return geometry.nZ[ix];
}

inline double
Kernel::Z(size_t ix,size_t iz){
  return geometry.Z[ix][iz];
}

#endif