/* 
   
   Copyright (C) Bernard Piette 
                 University of Durham (UK) 
                 Department of Mathematical Sciences
   
   This program is free software; you can redistribute it and/or modify
   it under the terms of the GNU General Public License as published by
   the Free Software Foundation; either version 2 of the License, or
   (at your option) any later version.
   
   You can modify or use any part of the source code but only under the  
   condition that this Copyright notice is kept in the resulting source file.
*/
#include "input_par.h"
#include <string>
#include <pde_1d.h>
#include <grf_pde_plain.h>
#include <iostream>

using namespace std;

const int NO_MAX  = 2;

const int NO_EQ   = 1;
const int NO_DENS = 1;
const int NO_SCALARS = 2*NO_MAX;

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

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

char const *Densities_Names[] =
{ "E", 0 
};

char const *Scalars_Names[] =
{ "Max1" , "Pos1", "Max2",  "Pos2", 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;
  double A;
  double l;
  int d;
} par;

char *help[] = {
 "radialHeat PAR_FILE: ",
 "radialHeat -h|-help: ouput this help message",
 "radialHeat -v|-version: ouput the program version number.",
 "  Solve the radial heat equation:",
 "   df/dt = d^2f/dr^2 + (d-1)/r df/dr - l/r^2 f",
 " 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",
 "   A DOUBLE",
 "   l DOUBLE",
 "   d INTEGER",
 0
};

double pi;

/* Radial heat equation */
class radialHeat_pde : public pde_1d_RK4
{ 
public:
  radialHeat_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] : fxx+(d-1)*fx/r-l*f/(x*x)*/
  /* 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] = fxx_[0]+(par.d-1)*fx_[0]/x_-par.l*f_[0]/(x_*x_);
  }

  /* Compute Densities                        */
  /* D[0] :E  : Energy density                */
  /********************************************/
  double Dens(double *D)
  { if(!fl_ || !fr_) // Desnisties are 0 at the edge of the box
    { for(int i=0; i < NoDens_; i++) D[i] = 0; 
      return(1.0); 
    }
    
    eval_derivatives();
    D[0] = fx_[0]*fx_[0]+par.l*f_[0]*f_[0]/(x_*x_);

    return(2.0*pi*x_);
  }

  /* Compute the scalar data:        */
  /* Max : local maximum values of E */
  /* MaxPos : position of maximum    */
  /***********************************/
  void Scalars() 
  { double fm[NO_MAX],xm[NO_MAX],dm,x;
    int i,j,m,k = 0;

    for(int i = 1; i < nx_-1; i++)
    { // D[i] > D[i-1] && D[i] > D[i+1] => Local maximum
      if((Densities_[i*NoEq_] > Densities_[(i-1)*NoEq_])&&
         (Densities_[i*NoEq_] > Densities_[(i+1)*NoEq_]))
      { if (k < NO_MAX) { j = k++;}
        else 
        { // Find the smalest maxima
          dm = fm[0]; j = 0; 
          for(m = 1 ; m< NO_MAX; m++)
          { if(fm[m] < dm) { dm = fm[m]; j = m;}
          }
          if(Densities_[i*NoEq_] < fm[j]) j = -1; // Local max too small
        }
        if(j >= 0)
        { fm[j] = Densities_[i*NoEq_]; 
          xm[j] = Xmin_+i*dx_;
        }
      }
    }

    // If k < NO_MAX : use edge value:
    if (Densities_[0] > Densities_[nx_-1]) { dm = Densities_[0]; x = Xmin_; }
    else                                  { dm = Densities_[nx_-1]; x = Xmax_;}
    for(i = k ; i< NO_MAX; i++) 
    { fm[i] = dm;
      xm[i] = x;
    }

    for(m = 0 ; m < NO_MAX; m++) 
    { Scalars_[2*m]   = fm[m];
      Scalars_[2*m+1] = xm[m];
    }
  };

  /* Implement some boundary conditions       */
  /********************************************/
  void Boundary(double *fct)
  { // df/dr is 0 at the origin if l=0
    if (par.l > 0) fct[0] = 0;
    else           fct[0] = fct[NoEq_];
  }

  /* Compute the initial condition. */
  /* Fcts_[i*NO_EQ]   : f           */
  /* Fcts_[i*NO_EQ+1] : df          */
  /**********************************/
  void Init(double x,double *F)
  { F[0] = par.A*exp(-par.k*x*x)*pow(x,sqrt(par.l));
  }
};

/*****************************************/
/* 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;

  pi = 4*atan(1.0);

  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);
    ip.get("A",par.A);
    ip.get("l",par.l);
    ip.get("d",par.d);

          // Build a pde for sineGordon
    radialHeat_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()*pde.dx()*par.dt_r);

          // Create a new graphic object
    grfp = new grf_pde_plain(&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;
  }
}