#ifndef PICCARD_HPP
#define PICCARD_HPP

#include <cstring>
#include "solution.hpp"
#include "horizontal_problem.hpp"
#include "all_vertical_richards.hpp"
#include "overland.hpp"

namespace Kernel{
  //! A class implementing Piccard algorithm
  class Piccard{
  private:
    static constexpr double tolerence_Piccard=1e-6;
    static constexpr double oscilation_Piccard=1e-4;
    static constexpr size_t max_iterations_Piccard=50;

    //! Geometry of the simulation
    const Geometry* geometry;

    //! A pointer to an HorizontalProblem object
    HorizontalProblem horizontal_problem;

    //! A pointer to an Overland object
    Overland overland;

    //! A pointer to an AllVerticalRichards object
    AllVerticalRichards all_vertical_richards;

    double *l;
    double *Pl;

    //! Compute l from h and hsat
    void compute_l(double* h,double* hsat,double& error);

    //! Compute Pl from l and P
    void compute_Pl(double** P);
  public:
    double dt;
    Solution* previous_solution;
    Solution* new_solution;
    bool has_converged;

    //! Construct an empty Piccard object
    Piccard();

    //! Init from a Geometry
    void init(const Geometry* geometry);

    //! Compute a new Solution form a previous once.
    //! The delta time depends of the number of subdivisions required by Solver::space_solution.
    void run();

  };

}


#endif