#ifndef PDE_1D_H
#define PDE_1D_H

#include <cmath>
#include <time.h>
#include <cstring>
#include <cstdio>
#include <string>

class grf_pde_1d;

class pde_1d {
public:
  pde_1d(double xmin,double xmax,int NoEq,int NoDens,int NoScalars,int nx,
           std::string file_name,char const *Fct_Names[],
           char const *Density_Names[],char const *Scalars_Names[]);
  virtual ~pde_1d();
  void init_grf(grf_pde_1d *grfp) { grfp_ = grfp; };

  // Method to Override for specific equation
  virtual void Eqn(double *F)=0;           // Equations 
  virtual double Dens(double *D)=0;        // Density function  
  virtual void Boundary(double *fct) {};   // Boundary function 
  virtual void Init(double x,double *F)=0; // Initial condition function
  virtual void Scalars() {};               // Compute scalars 

  void set_dt(double dt) { dt_ = dt;}

  // Get some parameters
  double dx()   { return(dx_);}
  double xmin() { return(Xmin_);}
  double xmax() { return(Xmax_);}
  int nx()      { return(nx_);}

  double *fcts()      { return(Fcts_);}
  double *densities() { return(Densities_);}
  double *intDens()   { return(IntDens_);}
  double *scalars()   { return(Scalars_);}

  /************* METHODS ****************/
  virtual void   initial_value();   
  virtual void   usefulVars();
  virtual void   computeDensities();
  virtual void   set_neighbour_fcts(double *f,int i);

  /************* OUTPUT DATA ****************/
  void write_header_scalars();
  void write_scalars();
  void write_grf();

  /************* INTEGRATION ****************/
  virtual double integrate_1step() {return(0.0);};
  virtual double integrate_by_steps(int n);
  virtual double integrate_until(double t);
  virtual double integrate_grf_by_steps(int n,int scalars_step,int grf_step);
  virtual double integrate_grf_until(double tmax,double scalar_dt,
                                     double grf_dt);

  /************* EVAL DERIVATIVES ****************/
  virtual void eval_derivatives();            // Evaluates fx_ and fxx_

protected:
  double dt_;
  double t_;
  double Xmin_;
  double Xmax_;
  double dx_;

  int NoEq_;
  int NoDens_;
  int NoScalars_;
  int nx_;

  // Useful variables
  double dxs_;       // dx*dx
  double tdx_;       // 2*dx
  int lastPoint_;    // nx_-1;
  int f_array_size_; // NoEq_*nx_

  // allocate arrays for the fcts and densities
  double *Fcts_;       // a vector of size nx* NoEq 
  double *Densities_;  // a vector of size nx* NoDensity 
  double *IntDens_;    // a vector of size NoDensity 
  double *Scalars_;    // a vector of size NoScalars 
  double *Force_;      // a vector of size NoEq 

  // Allocated arrays for Evaluated quantities 
  double *fx_,*fxx_,*fxxx_,*fxxxx_;  // df/fx, d^2f/dx^2 
  double *f_,*fl_,*fll_,*fr_,*frr_;  // local fcts ptrs.
  double x_;                 // local value of x = Xmin_+index*dx_
  int i_;                    // local value of grid index

  // Graphic variables
  grf_pde_1d *grfp_;
  int grf_index_;
  std::string file_name_;

  // Array of names for the Fcts, Densities and Data
  char const **Fct_Names_;     // [NoEq]
  char const **Density_Names_; // [NoDens]
  char const **Scalars_Names_; // [NoScalars]
};

/******************************** Euler ************************************/
class pde_1d_Euler : public pde_1d {

public:
  pde_1d_Euler(double xmin,double xmax,
               int NoEq,int NoDens,int NoScalars,int nx,std::string f_name,
               char const *Fct_Names[],char const *Density_Names[],
               char const *Scalars_Names[]);
  ~pde_1d_Euler();
  double integrate_1step();

protected:
  double *tmp1_;        // same as Fcts 
};

/******************************** RK4 *************************************/
class pde_1d_RK4 : public pde_1d {

public:
  pde_1d_RK4(double xmin,double xmax,
             int NoEq,int NoDens,int NoScalars,int nx,std::string f_name,
             char const *Fct_Names[],char const *Density_Names[],
             char const *Scalars_Names[]);
  ~pde_1d_RK4();
  double integrate_1step();
  void sub_RK4_1step(double *f,double *ef,double *mf,double *nef,
                     double Cnef,double Cmf);

protected:
  double *tmp1_;        // same as Fcts 
  double *tmp2_;        // same as Fcts 
  double *tmp3_;        // same as Fcts 
};

/******************************** PERIODIC PDE ***************************/

/******************************** Periodic Euler ***************************/
class pde_1d_Euler_periodic : public pde_1d_Euler {

public:
  pde_1d_Euler_periodic(double xmin,double xmax,
               int NoEq,int NoDens,int NoScalars,int nx,std::string f_name,
               char const *Fct_Names[],char const *Density_Names[],
               char const *Scalars_Names[]);
  virtual void  set_neighbour_fcts(double *f,int i);
};

/******************************** Periodic RK4 ***************************/
class pde_1d_RK4_periodic : public pde_1d_RK4 {

public:
  pde_1d_RK4_periodic(double xmin,double xmax,
             int NoEq,int NoDens,int NoScalars,int nx,std::string f_name,
             char const *Fct_Names[],char const *Density_Names[],
             char const *Scalars_Names[]);
  virtual void  set_neighbour_fcts(double *f,int i);
};
#endif