#include "input_par.h"
#include <string>
#include <pde_1d.h>
#include <grf_pde_1d.h>
#include <iostream>

#define NO_EQ      2
#define NO_DENS    3
#define NO_SCALARS 2

using namespace std;

char *VERSION = "sineGordon: v1.0; 13/05/03; Author: B. Piette";

char const *Fcts_Names[] =
{ "f", "df", 0
};

char const *Densities_Names[] =
{ "En", "Kin", "Pot", 0 
};

char const *Scalars_Names[] =
{ "Max" , "MaxPos", 0 
};

struct {
  double x_min;
  double x_max;
  int nx;
  double dt_r;
  double t_max;
  double grf_dt;
  double scalars_dt;
  char grf_name[256];

  double k;
} par;

char *help[] = {
 "sineGordon PAR_FILE: ",
 "sineGordon -h|-help: ouput this help message",
 "sineGordon -v|-version: ouput the program version number.",
 "  Solve the sine-Gordon equation.",
 " PAR_FILE:",
 "   x_min DOUBLE",
 "   x_max DOUBLE",
 "   nx INT",
 "   dt_r DOUBLE",
 "   t_max DOUBLE",
 "   grf_dt DOUBLE",
 "   scalars_dt DOUBLE",
 "   grf_name STRING",
 "   k DOUBLE",
 0
};


/* sine-Gordon Equation */
class sg_pde : public pde_1d_RK4 
{ 
public:
  sg_pde(double xmin,double xmax,
          int NoEq,int NoDens,int NoScalars,int nx,string f_name) :
    pde_1d_RK4(xmin,xmax,NoEq,NoDens,NoScalars,nx,f_name,
          &Fcts_Names[0],&Densities_Names[0],&Scalars_Names[0])
  {
  }

  /* Evaluate the Equation to solve */
  /* functions:                     */
  /* f[0] : f : sine Gordond fct    */
  /* f[1] : df : df/dt              */
  /* F : array for result           */
  /**********************************/
  void Eqn(double *F) 
  { 
    if(!fl_ || !fr_)  // Do nothing at the edge of the box
    { for(int i=0; i< NoEq_; i++) F[i] = 0; 
      return; 
    }

    eval_derivatives();

    F[0] = f_[1];
    F[1] = fxx_[0] -par.k*sin(f_[0]);
  }

  /* Compute Densities                        */
  /* D[0] :En  : Total energy density         */
  /* D[1] :Kin : Kinetic energy density       */
  /* D[2] :Pot : Potential energy density     */
  /********************************************/
  double Dens(double *D)
  { double fx,K,P;
    int i;

    if(!fl_ || !fr_) // Desnisties are 0 at the edge of the box
    { for(i=0; i < NoDens_; i++) D[i] = 0; 
      return(1.0); 
    }
    
    fx = (fr_[0]-fl_[0])/tdx_;
    K = 0.5*f_[1]*f_[1];
    P = 0.5*fx*fx+ par.k*(1.0-cos(f_[0]));
    D[0] = K+P;
    D[1] = K;
    D[2] = P;

    return(1.0);
  }

  /* Compute the scalar data:       */
  /* Max : maximum value of f       */
  /* MaxPos : position of maximum   */
  /**********************************/
  void Scalars() 
  { double fmax = -HUGE_VAL,xmax=Xmin_;
    
    for(int i = 0; i < nx_; i++)
    { if(Fcts_[i*NoEq_] > fmax) 
      { fmax = Fcts_[i*NoEq_]; 
        xmax = Xmin_+i*dx_;
      }
    }
    Scalars_[0] = fmax;
    Scalars_[1] = xmax;
  };

  /* Implement some boundary conditions       */
  /********************************************/
  void Boundary(double *fct)
  { // Nothing Here
  }

  /* Compute the initial condition. */
  /* Fcts_[i*NO_EQ]   : f           */
  /* Fcts_[i*NO_EQ+1] : df          */
  /**********************************/
  void Init(double x,double *F)
  { double xm;

    xm = (Xmin_+Xmax_)*0.5;

    F[0] = exp(-(x-xm)*(x-xm));
    F[1] = 0.0;
  }
};

/*****************************************/
/* Ouput a help message on cerr and exit */
/*****************************************/
void Help()
{ int i=0;
  while(help[i]) cerr << help[i++] << std::endl;
  exit(1);
}

int main(int argc, char **argv)
{ grf_pde_1d *grfp;

  try
  { // One argument expected.
    if((argc != 2) || !strcmp(argv[1],"-help") || !strcmp(argv[1],"-h")) 
    { Help();
    }
    if(!strcmp(argv[1],"-version") || !strcmp(argv[1],"-v")) 
    { cerr << VERSION<<"\n";
      exit(0);
    }


        // Read the parameters for the input file
        // The parameter filename is the 1st argument
    input_par ip(argv[1]);

    ip.get("x_min",par.x_min);
    ip.get("x_max",par.x_max);
    ip.get("nx",par.nx);
    ip.get("dt_r",par.dt_r);
    ip.get("t_max",par.t_max);
    ip.get("grf_dt",par.grf_dt);
    ip.get("scalars_dt",par.scalars_dt);
    ip.get("grf_name",par.grf_name);

    ip.get("k",par.k);

          // Build a pde for sineGordon
    sg_pde pde(par.x_min,par.x_max,NO_EQ,NO_DENS,NO_SCALARS,par.nx,
               par.grf_name);

          // Set dt by hand : dx * par.dt_r
    pde.set_dt(pde.dx()*par.dt_r);

          // Create a new graphic object
    grfp = new grf_pde_1d(&pde,par.grf_name);
    pde.init_grf(grfp);

    pde.initial_value();

    pde.integrate_grf_until(par.t_max,par.scalars_dt,par.grf_dt);
  }
  catch(string e) 
  { cerr << e << std::endl;
  }
}