// Compile using g++ -O or cpp -O or CC -O depending on architecture.
// random numbers might give trouble on some architectures, though...
//
// Main parameters are:         gxmax, gymax -  grid size
//                              nparticles   -  no of particles to deposit
//                              p_solid      -  sticking probabilities
//                              seed         -  random number seed
//
// Output files are:    field.dat (dump of the whole grid, x y color format)
//                      ordered.dat (chronological list of particle deposits.
//                                      x y no. format).
//                      fractal dimension as function of box size.
//
// Probably ordered.dat is most useful for fractal dimensions...
//
//
#include<iostream>
#include<fstream>
#include<math.h>
#include<stdlib.h>

//-------------------------------------------------------------------

void fileoutput(int world[], int xmax, int ymax);
void iterate(int world[], int neighbourx[][8], int neighboury[][8],
             int xmax, int ymax, int nparticles, double p_move,
             double p_solid, int colstep, int colrange );
void fractal(int world[], int xmax, int ymax, int nparticles);
void init(int world[], int neighbourx[][8], int neighboury[][8], int xmax,
          int ymax);

   
const int gxmax=300;    // size of the grid
const int gymax=300;

int gworld[gxmax*gymax];
int gneighbourx[gxmax*gymax][8];
int gneighboury[gxmax*gymax][8];
  
//-------------------------------------------------------------------

int main()
{
   double p_move=0.99;          // less than one to guarantee ergodicity

   double p_solid=0.9;          // sticking probability
   int seed=2;                  // random seed
   int nparticles=8000;         // number of particles
   
   int colstep=50;
   int colrange=10;
      
   srand(seed);
   init(gworld, gneighbourx, gneighboury, gxmax, gymax);  
   iterate(gworld,gneighbourx, gneighboury,gxmax,gymax,nparticles,
           p_move,p_solid,colstep,colrange);
   fractal(gworld,gxmax,gymax,nparticles);
   fileoutput(gworld,gxmax,gymax);

   return 0; 
}

//-------------------------------------------------------------------

void init(int world[], int neighbourx[][8], int neighboury[][8], int xmax,
          int ymax)
{
   int x,y,yy;
   int x1,x2,y1,y2;             // neighbour coordinates
   
   for(y=0;y<xmax*ymax;y+=xmax)         // zero everywhere
      for(x=0;x<xmax;x++)
         world[y+x]=0;

   y=(ymax/2)*xmax;                     // create condensation core
   x=xmax/2;
   world[y+x]=1;

   
   for(y=0;y<ymax;y++)         // initialize neighbour list
      for(x=0;x<xmax;x++)
      {
         x1= (x>0 ? x-1 : xmax-1); // find the neighbour x and y coordinates
         x2= (x<xmax-1 ? x+1 : 0);
         y1= (y>0 ? y-1 : ymax-1);
         y1*=xmax;
         y2= (y<ymax-1 ? y+1 : 0);
         y2*=xmax;
         yy=y*xmax;
         neighbourx[x+yy][0] = x1;   neighboury[x+yy][0] = y1;
         neighbourx[x+yy][1] = x1;   neighboury[x+yy][1] = yy;
         neighbourx[x+yy][2] = x1;   neighboury[x+yy][2] = y2;
         neighbourx[x+yy][3] = x;    neighboury[x+yy][3] = y1;
         neighbourx[x+yy][4] = x;    neighboury[x+yy][4] = y2;
         neighbourx[x+yy][5] = x2;   neighboury[x+yy][5] = y1;
         neighbourx[x+yy][6] = x2;   neighboury[x+yy][6] = yy;
         neighbourx[x+yy][7] = x2;   neighboury[x+yy][7] = y2;
      }
}

//-------------------------------------------------------------------

void iterate(int world[], int neighbourx[][8], int neighboury[][8],
             int xmax, int ymax, int nparticles, double p_move,
             double p_solid, int colstep, int colrange )
{
   int x,y;                     // walker coordinates
   int particle,dir,pos;
   bool walker;
   const double invmax=1.0/(RAND_MAX+1.0);

   std::ofstream myfile("ordered.dat");

   for(particle=1;particle<=nparticles;particle++)
   {  
      if (particle%100==0) cout<<particle<<endl;
      walker=true;
      if(rand()*invmax<0.5)
      {
         x=(int) floor(rand()*invmax*xmax); 
         y=0;
      }
      else
      {
         x=0;
         y=(int) floor(rand()*invmax*(ymax-1)); 
         y=(y+1)*xmax;
      }           
      
      do
      { 
         if(rand()*invmax<p_move)   // try to move
         {

            dir=(int) floor(rand()*invmax*8);               // pick direction
            pos=neighbourx[x+y][dir]+neighboury[x+y][dir];  // get position

            if(world[pos])                              // occupied location
            {
               if(rand()*invmax<p_solid)                // try to solidify
               {
                  world[x+y]=(particle/colstep)%colrange+1; // deposit particle
                  walker=false;                         // eliminate walker
                  myfile<<x<<" "<<y/xmax<<" "<<particle<<endl;
               }
            }
            else                                        // empty location
            {
               x=neighbourx[x+y][dir];                  // move walker
               y=neighboury[x+y][dir];
            }
 
         } /* p_move */
         
      } while(walker);

   }

   myfile.close();
}

//-------------------------------------------------------------------

void fractal(int world[], int xmax, int ymax, int nparticles)
{
   int xlow,ylow,x,y,bl,blmax;
   int sum,n=1;

   std::ofstream myfile("fractal.dat");
      
   xlow=xmax/2-1;                              // lower left corner
   ylow=ymax/2-1;       
   
   if(xmax>ymax) blmax=xmax;
   else blmax=ymax;
   
   for(bl=2;bl<blmax;bl+=2,xlow--,ylow--)
   {
      sum=0;
      y=ylow;
      for(x=xlow;x<=xlow+bl;x++)      
         if( world[y*xmax+x] ) sum++;
      y=ylow+bl;
      for(x=xlow;x<=xlow+bl;x++)      
         if( world[y*xmax+x] ) sum++;
      x=xlow;
      for(y=ylow+1;y<ylow+bl;y++)      
         if( world[y*xmax+x] ) sum++;
      x=xlow+bl;
      for(y=ylow+1;y<ylow+bl;y++)      
         if( world[y*xmax+x] ) sum++;
      if(sum==0) break;
         
      n+=sum;
      myfile << bl+1 <<" "<< n <<" "<< log(n)/log(bl+1) <<endl;
   }
   
   myfile.close();
}

//-------------------------------------------------------------------

void fileoutput(int world[], int xmax, int ymax)
{
   int x,y;
   
   std::ofstream myfile("field.dat");

   for(y=0;y<ymax;y++)
      for(x=0;x<xmax;x++)
         if( world[y*xmax+x] )
            myfile<<x<<" "<<y<<" "<<world[y*xmax+x]<<endl;

   myfile.close();
   
}

//-------------------------------------------------------------------