#include <fstream>
#include <string>
#include <random>
#include <iostream>
int main(int argc, char** argv) {
  const int     n = 1000,                     //lattice size
                half_n = .5*n,                //lattice size halved
                n_sq = n*n,                   //total number of sites
                n_min1 = n - 1,               //index of edge elements
                rid = std::stoi(argv[1]),     //index of simulation run
                t_max = 1000,                 //maximum time value
                dir_i[] = {-1, 0, 1, 0},      //set of neighbour index offsets in i direction
                dir_j[] = {0, -1, 0, 1},      //set of neighbour index offsets in j direction
                ipos_mod[] = {0, 0, 1, 1},    //arrays for calculating edge coordinates
                ipos_off[] = {0, n - 1, 0, 0},
                jpos_mod[] = {1, 1, 0, 0},
                jpos_off[] = {0, 0, 0, n - 1};

  int           lat[n][n];                    //lattice

  for(int i = 0; i < n; i++) {
    for(int j = 0; j < n; j++) {
      lat[i][j] = 0;
    }
  }

  lat[half_n][half_n] = 1;                          //particle in the middle of the lattice
  for(int i = 0; i < 4; i++) {                      //loop over the neighbours of the partcle
    lat[half_n + dir_i[i]][half_n + dir_j[i]] = 2;  //make sticky
  }  

  std::random_device rd;
  std::mt19937 mt(rd());
  std::uniform_int_distribution<int> U4(0, 3), Un(0, n - 1);

  for(int t = 0; t < t_max; t++) {
    int edge_id = U4(mt),  //random edge
        rnd = Un(mt),      //random start
        ipos = rnd*ipos_mod[edge_id] + ipos_off[edge_id], //adjust i
        jpos = rnd*jpos_mod[edge_id] + jpos_off[edge_id]; //adjust j

    while(lat[ipos][jpos] != 2) {                 //while the particle hasn't encountered a sticky patch
      int rnd_dir_id = U4(mt);                    //choose a random neighbour
      ipos = (ipos + dir_i[rnd_dir_id] + n) % n;  //move particle there
      jpos = (jpos + dir_j[rnd_dir_id] + n) % n;  //note the use of % for PBC (this is expensive, there are better ways to implement PBC)
    }

    lat[ipos][jpos] = 1;                            //by now the particle is at a sticky patch, so update lattice
    for(int i = 0; i < 4; i++) {                    //go over the neighbours of the stuck particle
      int nbr_i = (ipos + dir_i[i] + n) % n,        //calculate the neighbour coordinate with PBC
          nbr_j = (jpos + dir_j[i] + n) % n;
      if(lat[nbr_i][nbr_j] != 1) {                  //if the neighbour is not another particle
        lat[nbr_i][nbr_j] = 2;                      //make a sticky patch
      }
    }
    if(t % 100 == 0) {
      std::cout << "Timestep " << t << " of " << t_max << "\n";
    }
  }
  std::string fname = std::string("dlaframe_") + std::to_string(rid) + ".bin";
  std::ofstream ofile(fname.c_str(), std::ios::binary);             //open a binary file for output
  ofile.write((const char *) &lat[0][0], n_sq*sizeof(int));     //write the lattice there
  ofile.close();
  return 0;
}