#include <numeric>
#include <iostream>
#include <vector>
#include <cstdlib>
#include <cmath>
#include <time.h>
#include <sys/time.h>
#include <sstream>
#include <string>

const double ARmass = 39.948; //A.U.s
const double ARsigma = 3.405; // Angstroms
const double AReps   = 119.8; // Kelvins
const double CellDim = 12.0; // Angstroms
const int    NPartPerCell = 10;
double WallTime(void);
//vector <vector <int> > ComputeAtomsPerCell(int,int,int,int);

using namespace std;
// Class for keeping track of the properties for a particle
class Particle{
  public:
  double x;		// position in x axis
  double y;		// position in y axis
  double fx;		// total forces on x axis
  double fy;		// total forces on y axis
  double ax;		// acceleration on x axis
  double ay;		// acceleration on y axis
  
  double vx;		// velocity on x axis
  double vy;		// velocity on y axis
  // Default constructor
  Particle() 
    : x(0.0),y(0.0),fx(0.0),fy(0.0),ax(0.0),
      ay(0.0),vx(0.0),vy(0.0){
  }
  
  double update(){
    // We are using a 1.0 fs timestep, this is converted
    double DT = 0.000911633;
    ax = fx/ARmass;
    ay = fy/ARmass;
    vx += ax*0.5*DT;
    vy += ay*0.5*DT;
    x += DT*vx;
    y += DT*vy;
    fx = 0.0;
    fy = 0.0;
    return 0.5*ARmass*(vx*vx+vy*vy);
  }
};

class Cell {
public:
  vector<Particle> particles;
  Cell(int x, int y, int nParticles) : particles(nParticles){
    //  We will be filling the cells, making sure than
    //  No atom is less than 2 Angstroms from another
    double rcelldim = double(CellDim);
    double centerx = rcelldim*double(x) + 0.5*rcelldim;
    double centery = rcelldim*double(y) + 0.5*rcelldim;
   
    // Randomly initialize particles to be some distance 
    // from the center of the square
    // place first atom in cell
    particles[0].x = centerx + ((drand48()-0.5)*(CellDim-2));
    particles[0].y = centery + ((drand48()-0.5)*(CellDim-2));
    
    double r,rx,ry;
    for (int i = 1; i < particles.size(); i++) {
      r = 0;
      while(r < 4.0) {   // square of 2
	r = 4.00001;
	rx = centerx + ((drand48()-0.5)*(CellDim-2)); 
	ry = centery + ((drand48()-0.5)*(CellDim-2));
      	//loop over all other current atoms
	for(int j = 0; j<i; j++){
	  double rxj = rx - particles[j].x;
	  double ryj = ry - particles[j].y;
	  double rt = rxj*rxj+ryj*ryj;
	  if(rt < r) r=rt;
	}
      }
      particles[i].x = rx;
      particles[i].y = ry;
    }
  }
};

//Force calculation
//------------------------------------
double interact(Particle &atom1, Particle &atom2){
  double rx,ry,rz,r,fx,fy,fz,f;
  double sigma6,sigma12;

  // computing base values
  rx = atom1.x - atom2.x;
  ry = atom1.y - atom2.y;
  r = rx*rx + ry*ry;

  if(r < 0.000001)
    return 0.0;
  
  r=sqrt(r);
  sigma6 = pow((ARsigma/r),6);
  sigma12 = sigma6*sigma6;
  f = ((sigma12-0.5*sigma6)*48.0*AReps)/r;
  fx = f * rx;
  fy = f * ry;
  // updating particle properties
  atom1.fx += fx;
  atom1.fy += fy;
  return 4.0*AReps*(sigma12-sigma6);
}

int ComputeAtomsPerCell(int **sizes, 
			int NRows,int NCols, 
			int NParts){

  
  int max = NPartPerCell;
  for(int i=0;i<NRows;i++)
    for(int j=0;j<NCols;j++)
      sizes[i][j] = NPartPerCell;
  
  int molsum = NRows*NCols*NPartPerCell;
  while(molsum < NParts){
    max++;
    for(int i=0;i<NRows;i++){
      for(int j=0;j<NCols;j++){
	sizes[i][j]++;
	molsum++;
	if(molsum >= NParts) {
	  return max;
	}
      }
    }
  }
  return max;
}
			 

int main(int argc,char* argv[]){
  
  // Get Command Line Arguments
  if(argc!=3){
    cerr << "Incorrect syntax: should be two arguments";
    exit(2);
  }
  
  int NumParticles = atoi(argv[1]);
  int NumIterations  = atoi(argv[2]);
  
  int Dimension = int(sqrt(NumParticles/double(NPartPerCell)));
  int TotalCells = Dimension*Dimension;
  int NCellRows = Dimension;
  int NCellCols = Dimension;
  int MaxCellSize;

  int **CellSizes = (int **)malloc(sizeof(int*)*NCellRows);
  for(int i=0;i<NCellRows;i++)
    CellSizes[i] = (int *)calloc(sizeof(int),NCellRows);
  
  MaxCellSize = ComputeAtomsPerCell(CellSizes,NCellRows,
				    NCellCols,NumParticles);
  
  cout << "\nThe Total Number of Cells is " << TotalCells
       << " With a maximum of " << MaxCellSize << " particles per cell,"
       << "\nand "  << NumParticles << " particles total in system\n";
  
  // Allocate Cell, space for Particles couched in the Cell creator
  vector< vector< Cell > > cells(NCellRows);
  for (int i = 0; i < NCellRows; i++) {
    for (int j = 0; j < NCellCols; j++) {
      cells[i].push_back(Cell(i,j,CellSizes[i][j]));
    }
  }
  
  double TimeStart = WallTime();
  for (int t = 0; t < NumIterations; t++) { // For each timestep
    double TotPot=0.0;
    double TotKin=0.0;
    
    // For each cell...
    for (int cx1 = 0; cx1 < NCellRows; cx1++) {
      for (int cy1 = 0; cy1 < NCellCols; cy1++) {
        // Consider interactions between particles within the cell
        for (int i = 0; i < cells[cx1][cy1].particles.size(); i++) {
          for (int j = 0; j < cells[cx1][cy1].particles.size(); j++) {
	    if(i!=j)
	      TotPot += interact(cells[cx1][cy1].particles[i],
				 cells[cx1][cy1].particles[j]);
	  }
	}
	
	
	// Consider each other cell... 
	for (int cx2 = 0; cx2 < NCellRows; cx2++) {
	  for (int cy2 = 0; cy2 < NCellCols; cy2++) {
	    // if condition is true if cell[cy1][cy2] is neighbor (N,NE,E,SE,S,SW,W, or NW) of cell[cx1][cy1] 
	    if ((abs(cx1-cx2) < 2) && (abs(cy1-cy2) < 2) && 
		(cx1 != cx2 || cy1 != cy2)) {
	      // Consider all interactions between particles
	      for (int i = 0; i < cells[cx1][cy1].particles.size(); i++) {
		for (int j = 0; j < cells[cx2][cy2].particles.size(); j++) {
		  TotPot += interact(cells[cx1][cy1].particles[i],
				     cells[cx2][cy2].particles[j]);
		}
	      }
	    }
	  }
	} 
      }
    }
    
    // End iteration over cells    
    // Apply the accumulated forces; update accelerations, velocities, positions
    for (int cx1 = 0; cx1 < NCellRows; cx1++) {
      for (int cy1 = 0; cy1 < NCellCols; cy1++) {
	for(int i=0; i < cells[cx1][cy1].particles.size(); i++) {
          TotKin += cells[cx1][cy1].particles[i].update();
        }
      }
    }
    printf("\nIteration#%d with Total Energy %e per Atom",
	   t,(TotPot+TotKin)/NumParticles);
  } // iterate time steps
  double TimeEnd = WallTime();
  cout << "\nTime for " << NumIterations << " is "<< TimeEnd-TimeStart;
  
  return 0;
}

// A Simple timer for measuring the walltime
double WallTime(void){
  static long zsec =0;
  static long zusec = 0;
  double esec;
  struct timeval tp;
  struct timezone tzp;
  gettimeofday(&tp, &tzp);
  if(zsec ==0) zsec  = tp.tv_sec;
  if(zusec==0) zusec = tp.tv_usec;
  esec = (tp.tv_sec-zsec)+(tp.tv_usec-zusec)*0.000001;
  return esec;
}