#include <pde_1d.h>
#include <grf_pde_1d.h>
#include <string>
#include <iostream>
#include <fstream>
#include <cstdio>
#include <cmath>

using namespace std;

/**************************************************************************/
/* Constructor :                                                          */
/* size : initial array size                                              */
/* xmin,xmax : grid range                                                 */
/* NoEq,NoDens,NoScalars : No of equations, densities and scalars         */
/* data_prefix : prefix name for the data files.                          */
/* nx : No of grid points                                                 */
/* Fct_Names: const array of names for the fcts                           */
/* Density_Names: const array of names for the densities                  */
/* Fct_Names[] : array of str: name of each fct (do not delete!)          */
/* Density_Names[] : array of str: name of each density (do not delete!)  */
/* Scalars_Names[] : array of str: name of each Scalar (do not delete!)   */
/* Allocates all the fcts and densities arrays.                           */
/* gtype : type of grid: GRID_1D_REGULAR or GRID_1D_PERIODIC              */
/* throw string on error                                                  */
/**************************************************************************/
pde_1d::pde_1d(double xmin,double xmax,int NoEq,int NoDens,int NoScalars,
               int nx,string file_name,char const *Fct_Names[],
               char const *Density_Names[],char const *Scalars_Names[]) :
  t_(0.0),Xmin_(xmin),Xmax_(xmax),NoEq_(NoEq),NoDens_(NoDens),
  NoScalars_(NoScalars),nx_(nx),file_name_(file_name), 
  Fct_Names_(Fct_Names),
  Density_Names_(Density_Names), Scalars_Names_(Scalars_Names)
{ 
  dx_ = (Xmax_-Xmin_)/(nx_-1);
  dt_ = dx_*0.25;
  grf_index_ = 0;
    

  usefulVars();

  try 
  { Fcts_     = new double[nx_*NoEq_];
    Densities_  = new double[nx_*NoDens_];
    IntDens_    = new double[NoDens_];
    Scalars_    = new double[NoScalars_];
    Force_      = new double[NoEq_];
    fx_         = new double[NoEq_];
    fxx_        = new double[NoEq_];
    fxxx_       = new double[NoEq_];
    fxxxx_      = new double[NoEq_];
  } 
  catch(...) { throw((string)"No more memory!"); }
}

/**************/
/* Destructor */
/**************/
pde_1d::~pde_1d()
{ delete [] Fcts_;
  delete [] Densities_;
  delete [] IntDens_;
  delete [] Scalars_;
  delete [] Force_;
  delete [] fx_;
  delete [] fxx_;
  delete [] fxxx_;
  delete [] fxxxx_;
}

/************************ Evaluates various quantities  *****************/
/************************************/
/* Initialise the Fcts_             */
/* Calls Init() for every value x_i */
/************************************/
void pde_1d::initial_value()
{ int i;
  double x;

  x = Xmin_;

  for(i=0; i < nx_; i++)
  { Init(x,&Fcts_[i*NoEq_]);
     x += dx_;
  }
}

/*****************************/
/* Eval the useful variables */
/*****************************/
void pde_1d::usefulVars()
{ dxs_ = dx_*dx_;
  tdx_ = 2.0*dx_;
  lastPoint_ = nx_-1;
  f_array_size_ = NoEq_*nx_;
}

/*************************************************/
/* Eval the fcts derivatives                     */
/* On the edge of the grid use the inside points */
/* Uses the following variables:                 */
/* f_, fl_, fr_ and tdx and dxs.                 */
/* fl_ and fr_ can be null                       */
/*************************************************/
void pde_1d::eval_derivatives()
{ int i;

  if(!fl_) 
  { for(i=0; i< NoEq_; i++)
    { fx_[i] = (fr_[i]-f_[i])/dx_;
      fxx_[i] = (f_[i]+fr_[NoEq_+i]-2*fr_[i])/dxs_; 
    }
    return; 
  }
  if(!fr_) 
  { for(i=0; i< NoEq_; i++)
    { fx_[i] = (f_[i]-fl_[i])/dx_;
      fxx_[i] = (fl_[-NoEq_+i]+f_[i]-2*fl_[i])/dxs_; 
    }
    return; 
  }

  for(i=0; i< NoEq_; i++) 
  { fx_[i] = (fr_[i]-fl_[i])/tdx_;
    fxx_[i] = (fl_[i]+fr_[i]-2*f_[i])/dxs_;
  }
}

/***********************************************************/
/* Computes the densities and store them in Densities_     */
/* Dens(Vec) : a function which computes all               */
/*     the densities.                                      */
/*       use f_,fl_,fr_ : pointers to the vector of fcts   */
/*       Vec : vector where to store the results.          */
/*       i_  : lattice index                               */   
/*     return the integration mesure.                      */
/***********************************************************/
void pde_1d::computeDensities()
{ int i,d;
  double mesure;
  
     /* Initialise the Integration variables to 0 */
  for(d = 0; d < NoDens_; d++)
  { IntDens_[d] = 0.0;
  } 

        /* Scan the whole grid, including the first and the last point */
  for(i=0; i <= lastPoint_; i++)
  {     /* fct_vec : the fct at the lattice point i */
    set_neighbour_fcts(Fcts_,i);

        /* Computes the densities and store them in pde->Densities */
    mesure = dx_*Dens(&Densities_[i*NoDens_]); 

        /* Integrate each density */
    for(d = 0; d < NoDens_; d++)
    { IntDens_[d] += mesure*Densities_[i*NoDens_+d];
    } 
  }
}

/*******************************************/
/* Set the pointers to the neighbour fcts  */
/* f : fct array.                          */
/* i : index position of the grid          */
/*******************************************/
void pde_1d::set_neighbour_fcts(double *f,int i)
{ 
  x_ = Xmin_+i*dx_;
  i_ = i;

  f_   = &f[i*NoEq_];
  fl_  = &f[(i-1)*NoEq_];
  fll_ = &f[(i-2)*NoEq_];
  fr_  = &f[(i+1)*NoEq_];
  frr_ = &f[(i+2)*NoEq_];

  if(i==0)                  { fl_ =  0; fll_ = 0; }
  else if(i==1)             { fll_ = 0; }
  else if(i==(lastPoint_-1)){ frr_ = 0; }
  else if(i==lastPoint_)    { fr_ =  0; frr_ = 0; }
}

/************************ Data files  ****************************/
/**********************************************************/
/* Write a header describing the scalars columns:         */
/* # intDens1  intDens2 ...  scalars1 ...                 */
/* Append data to the existing file                       */
/* scalars_name : name of the ouput file                  */
/* throw m_except on error                                */
/**********************************************************/
void pde_1d::write_header_scalars()
{ int i;
  string scalars_name = file_name_+".s";
  ofstream ofs(scalars_name.c_str());

  if(!ofs)
  { throw((string)" Error: Can't open file : "+scalars_name);
  }

  // Write Fcts and Densities titles as comments
  ofs <<"# t";
  for(i = 0; i < NoDens_; i++)
  { ofs << " "<<Density_Names_[i];
  }
  for(i = 0; i < NoScalars_; i++)
  { ofs << " "<<Scalars_Names_[i];
  }
  ofs <<"\n";

}
/************************************************************/
/* Write the integrated densities and the scalars in a file */
/* Append data to the existing file                         */
/* scalars_name : name of the ouput file                    */
/* Output one line of the type:                             */
/* t intDens1  intDens2 ...  scalars1 ...                   */
/* throw m_except on error                                  */
/************************************************************/
void pde_1d::write_scalars()
{ int i;
  string scalars_name = file_name_+".s";
  ofstream ofs(scalars_name.c_str(),ios::app);

  if(!ofs)
  { throw((string)" Error: Can't open file : "+scalars_name);
  }

  ofs <<t_;

  // Write Integrated densities:
  for(i=0; i < NoDens_; i++)
  { ofs <<" "<< intDens()[i];
  }

  // Write Scalar data
  for(i=0; i < NoScalars_; i++)
  { ofs <<" "<< scalars()[i];
  }
  ofs <<"\n";
}

/***************************************************************/
/* Create a graphic file for the selected fcts and densities   */
/* Throw a string on error                                     */
/***************************************************************/
void pde_1d::write_grf()
{ int i;

  // Write fcts graphs
  for(i=0; i < NoEq_; i++)
  { grfp_->write_1_grf(Fct_Names_[i],grf_index_,&Fcts_[i],NoEq_,t_);
  }

  // Write Density graphs
  for(i=0; i < NoDens_; i++)
  { grfp_->write_1_grf(Density_Names_[i],grf_index_,&Densities_[i],NoDens_,t_);
  }

  grf_index_++;
}

/************************ Integration ****************************/
/**********************************************************/
/* Perform n integration steps on the grid                */
/* Update the time t_                                     */
/* pde : pde parameters                                   */ 
/* Return : the updated time                              */     
/**********************************************************/
double pde_1d::integrate_by_steps(int n)
{ while(n--)
  { integrate_1step();
  }
  return(t_);
}

/**********************************************************/
/* Perform n integration steps on the grid                */
/* Update the time t_                                     */
/* pde : pde parameters                                   */ 
/* Return : the updated time                              */     
/**********************************************************/
double pde_1d::integrate_until(double t)
{ 
  while(integrate_1step() < t);
  return(t_);
}

/****************************************************************/
/* Perform n integration steps on the grid with graphics        */
/* Output scalars every scalars_step                            */
/* Output graphs every grf_step                                 */
/* Update the time t_                                           */
/* n : No of integration steps                                  */
/* scalars_step : no of steps between scalars outputs           */
/* grf_step : no of steps between graphics                      */
/* Return : the updated time                                    */     
/****************************************************************/
double pde_1d::integrate_grf_by_steps(int n,int scalars_step,int grf_step)
{ int i;
  string grf_name;

  write_header_scalars();

  i = 0;
  while(i < n)
  { if(((i%scalars_step) == 0) || ((i%grf_step) == 0))
    { computeDensities();
      Scalars();
      if((i%grf_step) == 0)   write_grf();
      if((i%scalars_step) == 0) write_scalars();
    }
    integrate_1step();
    i++;
  }

  if(((i%scalars_step) == 0) || ((i%grf_step) == 0))
  { computeDensities();
    Scalars();
    if((i%grf_step) == 0)   write_grf();
    if((i%scalars_step) == 0) write_scalars();    
  }
  return(t_);
}

/****************************************************************/
/* Perform n integration step on the grid with graphics         */
/* Output scalars every scalars_dt                              */
/* Output graphs every grf_dt                                   */
/* Update the time t_                                           */
/* tmax : integration bound.                                    */
/* scalars_dt : time interval between scalars outputs           */
/* grf_dt : time interval between graphics                      */
/* Return : the updated time                                    */     
/****************************************************************/
double pde_1d::integrate_grf_until(double tmax, double scalar_dt,double grf_dt)
{ int i;
  string grf_name;
  double scalar_t,grf_t;
  write_header_scalars();

  scalar_t = t_ - dt_*0.01;
  grf_t = t_ - dt_*0.01;

  while(t_ < tmax)
  { if((t_ >= scalar_t) || (t_ > grf_t))
    { computeDensities();
      Scalars();
      if(t_ >= grf_t)   { write_grf(); grf_t += grf_dt; }
      if(t_ >= scalar_t){ write_scalars(); scalar_t += scalar_dt; }
    }
    integrate_1step();
    i++;
  }

  if((t_ > scalar_t) || (t_ > grf_t))
  { computeDensities();
    Scalars();
    if((t_ >= grf_t) == 0)  write_grf();
    if((t_ >= scalar_t) == 0) write_scalars(); 
  }
 
  return(t_);
}

/********************* Euler ***********************************/
/********************************/
/* Constructor for Euler pde_1d */
/* Allocated 2 extra arrays     */
/********************************/
pde_1d_Euler::pde_1d_Euler(double xmin,double xmax,int NoEq,int NoDens,
          int NoScalars,int nx,string data_prefix,char const *Fct_Names[],
          char const *Density_Names[],char const *Scalars_Names[]) :
    pde_1d(xmin,xmax,NoEq,NoDens,NoScalars,nx,data_prefix,Fct_Names,
              Density_Names,Scalars_Names)
{ tmp1_ = new double[nx_*NoEq_];    
}

pde_1d_Euler::~pde_1d_Euler()
{ delete [] tmp1_;
}

/**********************************************************/
/* Perform a single euler integration step on the grid    */
/* Update the time t_                                     */
/* pde : pde parameters                                   */ 
/* Eqn(Force) : computes the right hand side of           */
/*     the equation.                                      */
/*       uses f_,fl_,fr_: local fct pointer               */
/*       Force : vector where to store the result.        */
/*       i_    : lattice index                            */   
/* Return : the updated time                              */     
/**********************************************************/
double pde_1d_Euler::integrate_1step()
{ int i,j,grid_size;
  
  grid_size = NoEq_*nx_;

      /* Scan the whole grid */
  for(i=0; i <= lastPoint_; i++)
  {     /* fct_vec : the fct at the lattice point i */
    set_neighbour_fcts(Fcts_,i);

        /* Computes the force for point i and  */
        /* stores it in Force                  */
    Eqn(Force_); 

    /* Update each fct using Force : f(n+1) = f(n) + Force*dt */
    /* tmp1 is used to store the fct temporarily              */
    for (j = 0; j < NoEq_; j++)
    { tmp1_[i*NoEq_+j] = f_[j]  + dt_*Force_[j];
    }
  }

  /* Update Fct with the new fct */
  for(i=0; i < grid_size; i++)
  { Fcts_[i] = tmp1_[i];
  }

  Boundary(Fcts_);
  t_ += dt_;  /* Udate the time */
  return(t_);
}



/********************* RK4 ***********************************/
/******************************/
/* Constructor for RK4 pde_1d */
/* Allocated 2 extra arrays   */
/******************************/
pde_1d_RK4::pde_1d_RK4(double xmin,double xmax,int NoEq,int NoDens,
          int NoScalars,int nx,string data_prefix,char const *Fct_Names[],
          char const *Density_Names[],char const *Scalars_Names[]) :
       pde_1d(xmin,xmax,NoEq,NoDens,NoScalars,nx,data_prefix,Fct_Names,
              Density_Names,Scalars_Names)
{ tmp1_ = new double[nx_*NoEq_];    
  tmp2_ = new double[nx_*NoEq_];    
  tmp3_ = new double[nx_*NoEq_];
}

pde_1d_RK4::~pde_1d_RK4()
{ delete [] tmp1_;
  delete [] tmp2_;
  delete [] tmp3_;
}

/**********************************************************/
/* Perform a 4th order Runge-Kutta integration step on    */
/* f(n+1) = fn + 1/6(k1 +2 k2 + 2 k3 + k4)                */
/* k1 = dt F(t,fn)                                        */
/* k2 = dt F(t+dt/2,fn+k1/2)                              */
/* k3 = dt F(t+dt/2,fn+k2/2)                              */
/* k4 = dt F(t+dt,fn+k3)                                  */
/* the grid. Update the time in pde.                      */
/* Eqn(Force) : a function which computes the right       */
/*     hand side of the equation.                         */
/*       Force : vector where to store the result.        */
/*       i_    : lattice index                            */   
/* Return : the updated time                              */     
/**********************************************************/
double pde_1d_RK4::integrate_1step()
{ int i;

    /* tmp1 is used to compute fn + 1/6(k1 + 2 k2 + 2 k3 + k4) */
    /* We initialise it to 0                                */
  for(i=0; i < f_array_size_; i++)
  { tmp1_[i] = 0;
  }

    /* tmp1 = 0 */
  sub_RK4_1step(Fcts_,Fcts_,tmp1_,tmp2_,0.5,1.0);

    /* tmp1 = k1 */
    /* tmp2 = fn + 1/2 k1 */
  t_ += 0.5*dt_; 
  sub_RK4_1step(Fcts_,tmp2_,tmp1_,tmp3_,0.5,2.0);

    /* tmp1 = fn + 1/6 (k1 + 2 k2 ) */
    /* tmp3 = fn + 1/2 k2 */
  sub_RK4_1step(Fcts_,tmp3_,tmp1_,tmp2_,1.0,2.0);

    /* tmp1 = fn + 1/6 (k1 + 2 k2 + 2 k3) */
    /* tmp3 = fn + k3 */
  t_ += 0.5*dt_; 
  sub_RK4_1step(Fcts_,tmp2_,tmp1_,0,0.0,1.0);

    /* tmp1 = (k1 + 2 k2 + 2 k3 + k4) */
    /* we update fn => fn + tmp1 / 6  */
  for(i=0; i < f_array_size_; i++)
  { Fcts_[i] += tmp1_[i]/6.0;
  }
  Boundary(Fcts_);

  /* Time has been updated before the last step */
  return(t_);
}

/**************************************************************/
/* Perform one RK intermediate step                           */
/* f : f(n) the fcts at time tn                               */
/* ef : estimated f for this step.                            */
/* mf : modified f : use add the contribution of each step    */
/* nef : estimated f for the next step, must be computed here */
/* Cnef : coeficiant used to computed nef                     */ 
/* Cmf  : coeficiant used to computed mf                      */
/**************************************************************/
void pde_1d_RK4::sub_RK4_1step(double *f,double *ef,double *mf,double *nef,
                             double Cnef,double Cmf)
{ int index,i,j;

  for(i=0; i <= lastPoint_; i++)
  {     /* fct_vec : the fct at the lattice point i */

    index = i*NoEq_;

    set_neighbour_fcts(ef,i);

        /* Computes the force for point i and  */
        /* stores it in Force                  */
    Eqn(Force_); 

    /* Update each fct using Force : f(n+1) = f(n) + Force*dt */
    for (j = 0; j < NoEq_; j++)
    { mf[index+j] += dt_ * Force_[j] * Cmf;

      if(nef) /* nef will be 0 for the last step */ 
      { nef[index+j] = f[index+j] + dt_ * Force_[j] * Cnef;
      }
    }
  }

    /* Set the boundary condition for the estimated fct */
  if(nef) Boundary(nef);
}

/******************** PERIODIC PDE ********************************/
/*****************************************/
/* Constructor for periodic Euler pde_1d */
/* Allocated 2 extra arrays              */
/*****************************************/
pde_1d_Euler_periodic::pde_1d_Euler_periodic(double xmin,double xmax,int NoEq,
          int NoDens,
          int NoScalars,int nx,string data_prefix,char const *Fct_Names[],
          char const *Density_Names[],char const *Scalars_Names[]) :
    pde_1d_Euler(xmin,xmax,NoEq,NoDens,NoScalars,nx,data_prefix,Fct_Names,
              Density_Names,Scalars_Names)
{ dx_ = (Xmax_-Xmin_)/nx_;
}

/*******************************************/
/* Set the pointers to the neighbour fcts  */
/* When using periodic boundray conditions */
/* f : fct array.                          */
/* i : index position of the grid          */
/*******************************************/
void pde_1d_Euler_periodic::set_neighbour_fcts(double *f,int i)
{ 
  x_ = Xmin_+i*dx_;
  i_ = i;

  f_   = &f[i*NoEq_];
  fl_  = &f[(i-1)*NoEq_];
  fll_ = &f[(i-2)*NoEq_];
  fr_  = &f[(i+1)*NoEq_];
  frr_ = &f[(i+2)*NoEq_];

  if(i==0)            { fl_ = &f[(nx_-1)*NoEq_] ; fll_ = &f[(nx_-2)*NoEq_]; }
  else if(i==1)             { fll_ = &f[(nx_-1)*NoEq_]; }
  else if(i==(lastPoint_-1)){ frr_ = &f[0]; }
  else if(i==lastPoint_)    { fr_ =  &f[0]; frr_ = &f[NoEq_]; }
}

/***************************************/
/* Constructor for periodic RK4 pde_1d */
/* Allocated 2 extra arrays            */
/***************************************/
pde_1d_RK4_periodic::pde_1d_RK4_periodic(double xmin,double xmax,int NoEq,
          int NoDens,int NoScalars,int nx,string data_prefix,
          char const *Fct_Names[],
          char const *Density_Names[],char const *Scalars_Names[]) :
       pde_1d_RK4(xmin,xmax,NoEq,NoDens,NoScalars,nx,data_prefix,Fct_Names,
              Density_Names,Scalars_Names)
{ dx_ = (Xmax_-Xmin_)/nx_;
}

/*******************************************/
/* Set the pointers to the neighbour fcts  */
/* When using periodic boundray conditions */
/* f : fct array.                          */
/* i : index position of the grid          */
/*******************************************/
void pde_1d_RK4_periodic::set_neighbour_fcts(double *f,int i)
{ 
  x_ = Xmin_+i*dx_;
  i_ = i;

  f_   = &f[i*NoEq_];
  fl_  = &f[(i-1)*NoEq_];
  fll_ = &f[(i-2)*NoEq_];
  fr_  = &f[(i+1)*NoEq_];
  frr_ = &f[(i+2)*NoEq_];

  if(i==0)            { fl_ = &f[(nx_-1)*NoEq_] ; fll_ = &f[(nx_-2)*NoEq_]; }
  else if(i==1)             { fll_ = &f[(nx_-1)*NoEq_]; }
  else if(i==(lastPoint_-1)){ frr_ = &f[0]; }
  else if(i==lastPoint_)    { fr_ =  &f[0]; frr_ = &f[NoEq_]; }
}