#ifndef RICHARDS_HPP
#define RICHARDS_HPP

#include <iostream>
#include <limits>
#include <cassert>
#include "physics.hpp"
#include "math/algo.hpp"
#include "debug.hpp"

using namespace std;

namespace Kernel{
  class Richards{
  private:

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

    double* pumps;

    double* sup_A;
    double* diag_A;
    double* sub_A;
    double* sup_B;
    double* diag_B;
    double* sub_B;
    double* sup_C;
    double* diag_C;
    double* sub_C;
    double* F;
    double* G;
    double* S;
    double* H;
    double* R;
    double* I;
    double* BB;

    size_t ix;
    size_t nZ;
    double dZ;

    double norm_previous_P;
    bool weighted_run(double w);
    void solve_system();
    double* temp_P[2];
    double* in_P;
    double* out_P;
  public:
    double dt;
    double flux_bot;
    double* div_w;

    double* near_P; //for Newton algorithm
    double* previous_P; //previous relatively to n2
    double* P;
    bool has_converged;

    Richards();
    void init(size_t ix,size_t nZ,double dZ,double* pumps);

    void run();
  };


}
#endif