#include <cstdlib>
#include <cmath>
#include <iostream>
#include <string>
#include <vector>
#include <fstream>
#include <plotter.h>

template <class T>
inline T sqr(const T& x){return x*x;}

class Lattice{
  int N;                                // size of the lattice is NxN
  std::vector<std::vector<bool> > site; // matrix of positions, if it is false, it is empty
  std::pair<int,int> tpos, ppos;        // temporary and previous position of the particle
  std::pair<double,double> r0;          // average position <vec{r}>
  double r2;                            // average <r^2>
  static const int MAXRGB = 65535;
  int Nclust;                          // number of points in the cluster
public:
  Lattice(int N_);
  void ReleaseNew();
  int MakeStep();
  void Draw(Plotter& plotter, int pixsize);
  void DrawAStep(Plotter& plotter, int pixsize);
  double r_average(){return sqrt(sqr(r0.first)+sqr(r0.second))/Nclust;}// size of |<r>|
  double gyration(){return r2/Nclust-sqr(r_average());}
  //The radius of gyration, rgyr, is the average distance from the centre of a cluster to any member of the cluster
  int size(){return Nclust;}
  void PrintMass();
};

Lattice::Lattice(int N_): N(N_), site(N), r0(0,0), r2(0), Nclust(0)
{
  for (int i=0; i<site.size(); i++){
    site[i].resize(N);
    for (int j=0; j<site[i].size(); j++){
      site[i][j]=false;                // sites are empty at the beginning
    }
  }
  site[N/2][N/2]=true;// We put first point in the middle of the system
};

void Lattice::ReleaseNew()
{ // Random walker is released from the boundary
  do{
    int istart = static_cast<int>(drand48()*4*N); // boundary size is 4*N
    switch(istart/N){ // which face of a square the particle is put to?
    case 0: tpos = std::make_pair(istart,0); break;
    case 1: tpos = std::make_pair(N-1,istart%N); break;
    case 2: tpos = std::make_pair(istart%N,N-1); break;
    case 3: tpos = std::make_pair(0,istart%N); break;
    }
    // we have random site on the bundary. But is it empty?
  } while (site[tpos.first][tpos.second]);// just make sure that the site is empty
}

int Lattice::MakeStep()
{// Return codes: 1 - particle glued
  //              0 - particle moved
  std::pair<int,int> newpos(tpos);
  int ipos = static_cast<int>(drand48()*4); // particle can go in 4 directions: up,down,left,right
  switch(ipos){// Here we use peridoci boundary conditions in order that particle can not escape too fast
  case 0: newpos.first  = (tpos.first+1)%N;   break; // right
  case 1: newpos.second = (tpos.second+1)%N;  break; // up
  case 2: newpos.first  = tpos.first>0  ? (tpos.first-1) : N-1;   break; // left
  case 3: newpos.second = tpos.second>0 ? (tpos.second-1) : N-1;  break; // down
  }
  // Checking weather particle can be glued!
  // Here we do not take periodic boundary conditions since they put too many particles on the boundary
  bool r_neighbor = newpos.first<N-1   && site[newpos.first+1][newpos.second];
  bool u_neighbor = newpos.second<N-1  && site[newpos.first][newpos.second+1];
  bool l_neighbor = newpos.first>0     && site[newpos.first-1][newpos.second];
  bool d_neighbor = newpos.second>0    && site[newpos.first][newpos.second-1];
  
  if (r_neighbor || u_neighbor || l_neighbor || d_neighbor){ // If particle can be glued
    // Particle just glued
    site[newpos.first][newpos.second]=true; // lattice is now occupied
    r0.first += newpos.first-N/2;
    r0.second += newpos.second-N/2;
    r2 += sqr(newpos.first-N/2)+sqr(newpos.second-N/2); // average r^2
    Nclust++;                                   // number of particles in DLA cluster increases
    return 1;
  }
  ppos = tpos;
  tpos = newpos; // update position of the walker
  return 0;
}

void Lattice::Draw(Plotter& plotter, int pixsize)
{
  plotter.filltype(1); // filled boxes
  plotter.fillcolor(0,0,0);
  for (int i=0; i<N; i++){
    for (int j=0; j<N; j++){
      if (site[i][j]){ // plot a box if particle is there
	plotter.fbox(i*pixsize,j*pixsize,(i+1)*pixsize,(j+1)*pixsize);
      }
    }
  }
}

void Lattice::DrawAStep(Plotter& plotter, int pixsize)
{
  int i = tpos.first, j = tpos.second; // plots a point at temporary position tpos
  plotter.fillcolor(0,0,MAXRGB);// in blue color
  plotter.box(i*pixsize,j*pixsize,(i+1)*pixsize,(j+1)*pixsize);//box at position (i,j)
}

void Lattice::PrintMass()
{
  std::ofstream out("mass.dat");
  double dr = 10;   // concentric rings of radium  [i*rd,(i+1)*dr] are checked for number of points inside (mass)
  double Rmax = sqrt(2)*site.size()/2.; // maximum radius the particle could be
  int Nr = static_cast<int>(Rmax/dr)+1;   // number of rings
  std::vector<double> mass(Nr);           // rings mass stored
  for (int i=0; i<mass.size(); i++) mass[i]=0; // initialization
  
  for (int i=0; i<site.size(); i++){
    for (int j=0; j<site[i].size(); j++){
      if (site[i][j]){ // this site is occupied
	double rx = i-N/2;
	double ry = j-N/2;
	double r = sqrt(rx*rx+ry*ry);   // distance from the origin
	mass[static_cast<int>(r/dr)]++; // the ring at r/dr containe one more particle
      }
    }
  }
  for (int i=1; i<mass.size(); i++) mass[i] += mass[i-1]; // sumulative sum -> mass inside the circle of radius dr*(i+1)
  for (int i=0; i<mass.size(); i++) out<<dr*(i+1.)<<" "<<mass[i]<<endl;
}

int main(int argc, char *argv[], char *env[])
{
  int N=200;
  int M=3000;
  int pixsize=10;
  bool debug=false;
  int i=0;
  while (++i<argc){
    std::string str(argv[i]);
    if (str=="-N" && i<argc-1)  N = atoi(argv[++i]);
    if (str=="-M" && i<argc-1)  M = atoi(argv[++i]);
    if (str=="-ps" && i<argc-1)  pixsize = atoi(argv[++i]);
    if (str=="-debug") debug=true;
    if (str=="-h" || str=="--help"){
      std::clog<<"************** Dendritic cluster grow **************\n";
      std::clog<<"**                                                **\n";
      std::clog<<"**      Copyright Kristjan Haule, 26.09.2005      **\n";
      std::clog<<"****************************************************\n";
      std::clog<<"\n";
      std::clog<<"dla [-N int] [-h]\n" ;
      std::clog<<"Options:   -h       Displays this help message\n";
      std::clog<<"           -N       Linear dimention of the lattice ("<<N<<")\n";
      std::clog<<"           -M       Number of random walkers ("<<M<<")\n";
      std::clog<<"           -ps      Pixel size ("<<pixsize<<")\n";
      std::clog<<"           -debug   Each step of simulation is plotted ("<<debug<<")\n";
      std::clog<<"*****************************************************\n";
    return 0;
    }
  }

  /************** Initialization of plotter **************************/
  PlotterParams params; // set a Plotter parameter
  params.setplparam ("PAGESIZE", (char *)"letter");
  XPlotter plotter(cin, cout, cerr, params); // declare Plotter
  if (plotter.openpl () < 0){ // open Plotter      
    cerr << "Couldn't open Plotter\n";
    return 1;
  }
  plotter.fspace (0.0, 0.0, pixsize*N, pixsize*N); // specify user coor system
  plotter.erase ();                // erase Plotter's graphics display
  /*************************************************************************/

  /************** Initialization of random number gen. *********************/
  int random_seed = time(0);
  srand48(random_seed);
  /*************************************************************************/
  
  /************** Initialization of main class  ***************************/
  Lattice lattice(N);
  /*************************************************************************/
  
  /************** Simulation  ***************************/
  cout<<"#"<<" "<<"|<r>|"<<" "<<"<r^2>-<r>^2"<<endl;
  for (int i=0; i<M; i++){
    lattice.ReleaseNew();
    while (lattice.MakeStep()==0) if (debug) lattice.DrawAStep(plotter, pixsize);
    if (debug) plotter.erase ();
    if (debug || i%20==0) lattice.Draw(plotter,pixsize);
    if (i%100==0){
      cout<<i<<" "<<lattice.r_average()<<" "<<lattice.gyration()<<endl;
    }
  }
  /******************************************************/
  cout<<M<<" "<<lattice.r_average()<<" "<<lattice.gyration()<<endl;
  
  lattice.PrintMass();

  clog<<"DONE"<<endl;
  if (plotter.closepl () < 0){ // close Plotter
    cerr << "Couldn't close Plotter\n";
    return 1;
  }
  return 0;
}