#include <windows.h>
#include <gl/glut.h>
#include <cstdlib>
#include <climits>
#include <fstream>
#include <iostream>
#include <cstdio>
#include <cmath>

using namespace std;

// Soft simulation gives more precise results.
#define SOFT_SIMULATION
//#define DEBUG


// Real number can be either double or float; I prefer float for
// more precission.
typedef double real;

// 2D vector structure.
struct Vec2 {
  real x;
  real y;
};

// A structure representing one molecule in defined space.
struct Molecule {
  // Current position of molecule.
  Vec2 position;
  
#ifdef SOFT_SIMULATION
  // Potencial of molecule.
  real potencial;
#endif
};

// Finds the nearest radius.
real nearRadius2(Vec2 p1, const Vec2& p2, const real& a) {
  real a_2 = a * 0.5;
  real xdif = p2.x - p1.x;
  real ydif = p2.y - p1.y;
  
  // We check if more then a/2 away in x.
  if(xdif > a_2) p1.x += a;
  else if(xdif < -a_2) p1.x -= a;
  
  // We check if more then a/2 away in y.
  if(ydif > a_2) p1.y += a;
  else if(ydif < -a_2) p1.y -= a;
  
  // We return the radius squared.
  Vec2 res;
  res.x = p2.x - p1.x;
  res.y = p2.y - p1.y;
  return res.x*res.x + res.y*res.y;
}

#ifdef SOFT_SIMULATION
real calcPotencial(unsigned int n, unsigned int i, Molecule* m, 
                   const Vec2& newPos, real maxr, real a, real epsilon, real sigma) {
  unsigned int j;
  real potencial = 0.0;
  real maxr2 = maxr*maxr;
  
  // We calculate potencial (before i)
  for(j = 0; j < i; j++) {
    // Get radius squared. Ignore all not in maxr range.
    real radSq = nearRadius2(m[j].position, newPos, a);
    if(radSq > maxr2) continue;
    
    // Get real radius
    real rad = sqrt(radSq);
    
    // We calculate potencial
    real p = 4.0*epsilon * (pow(sigma/rad, 12.0) - pow(sigma/rad, 6));
    potencial += p;
  }
  
  // We calculate potencial (after i)
  for(j = i+1; j < n; j++) {
    // Get radius squared. Ignore all not in maxr range.
    real radSq = nearRadius2(m[j].position, newPos, a);
    if(radSq > maxr2) continue;
    
    // Get real radius
    real rad = sqrt(radSq);
    
    // We calculate potencial
    real p = 4.0*epsilon * (pow(sigma/rad, 12.0) - pow(sigma/rad, 6));
    //printf("Epsilon is %f, radius %f, sigma %f and probability is %f\n", epsilon, rad, sigma, p);
    potencial += p;
  } 
    
  // We return a potencial.  
  return potencial; 
}

#endif

// Returns random number in interval [0,1].
inline real frand() {
  return (real)rand() / (real)RAND_MAX;  
}

// Sets the position in range.
inline void toRange(Vec2& pos, real a) {
  if(pos.x < 0) pos.x += a;
  if(pos.x > a) pos.x -= a;
  if(pos.y < 0) pos.y += a;
  if(pos.y > a) pos.y -= a;
}

// Simulates one step. Returns 1 if succeeded and 0 if not.
int simulate(unsigned int n, Molecule* m, real radius, real a, real maxr
#ifdef SOFT_SIMULATION
  , real eps, real sigma, real T
#endif
  ) {
  // This asumption is ok, if 2^n particles are used and
  // n is less then RAND_MAX.
  unsigned int i = rand() % n;
  Molecule& mol = m[i];
  real radius2 = radius*radius;
  
  // We have a random molecule, we make a random move.
  Vec2 newPos;
  newPos.x = mol.position.x + (frand()*2.0-1.0)*maxr;
  newPos.y = mol.position.y + (frand()*2.0-1.0)*maxr;
  toRange(newPos, a);
  
#ifdef SOFT_SIMULATION
  // We check if molecule is likely to move there.
  real newPotencial = calcPotencial(n, i, m, newPos, maxr, a, eps, sigma);
  
  // If delta < 0.0, then the step will happen because molecule want to go
  // to lower potencial energies.
  real delta = newPotencial - mol.potencial;
  mol.potencial = newPotencial;
  if(delta > 0.0) {
    // We allow jumping to higher potencials based on probability
    // of e^(delta/kT) in range [0,1]. We let k = 1.
    return 0;
    /*real prob = exp(-delta/T);
    printf("Probability is %f, based on dela %f and temperature %f\n", prob, delta, T);
    if(frand() > prob) return 0;*/
  }
#else
  // We use two loops to avoid ifs in loop.
  for(unsigned int j = 0; j < i; j++) {
    real radSq = nearRadius2(m[j].position, newPos, a);
    if(radSq < radius2) {
#ifdef DEBUG
      printf("Not able to move\n");
#endif
      return 0;      
    }    
  }
  
  for(unsigned int z = i+1; z < n; z++) {
    real radSq = nearRadius2(m[z].position, newPos, a);
    if(radSq < radius2) {
#ifdef DEBUG
      printf("Not able to move\n");
#endif
      return 0;      
    }        
  }
#endif

  // Move to new position.
  mol.position = newPos;
#ifdef DEBUG
  printf("Moved to new position (%f, %f)\n", newPos.x, newPos.y);
#endif
  return 1;
}

// Function initializes space, given the number of particles n,
// the space defined in cartesian coordinares from ([0,a], [0,a])
// and collection of molecules that is allocated into array.
void initSpace(unsigned int n, real a, Molecule*& m, real radius
#ifdef SOFT_SIMULATION
  , real epsilon, real sigma, real maxr
#endif
  ) {
  // We dynamically initialize the array.
  m = new Molecule[n];
  real radius2 = radius*radius;
  
  // We init the molecules.
  for(unsigned int i = 0; i < n; i++) {
    // We process the i-th molecule by setting it's position somewhere
    // in defined range
    for(;;) {
      // We create random position
      Vec2 randpos;
      randpos.x = frand()*a;
      randpos.y = frand()*a;
      
      // We check if position is ok; we repeat it if not ok.
      for(unsigned int j = 0; j < i; j++) {
        Molecule& mol = m[j];
        // We calculate dot product.
        real dot = randpos.x*mol.position.x + randpos.y*mol.position.y;
        
        // Sample is not ok, if the molecules are too close.
        if(dot <= radius2) continue;
                   
      } 
      
      m[i].position = randpos;
#ifdef DEBUG
      printf("Setting molecule position to (%f, %f)\n", randpos.x, randpos.y);
#endif
      
#ifdef SOFT_SIMULATION
      m[i].potencial = calcPotencial(n, i, m, m[i].position, maxr, a, epsilon, sigma);
#endif
      break;      
    }           
  }  
}

real reevaluateMaxR(real maxR, unsigned int M, unsigned int S) {
  if(M*2 > S) return maxR*0.9;
  return maxR*1.1;    
}

// Prepares transformations for rendering.
void renderPrepare(real a) {
  glMatrixMode(GL_PROJECTION);
  glLoadIdentity();
  glScalef(2.0f/a, 2.0f/a, 1.0f);
  glTranslatef(-1.0f, -1.0f, 0.0f);
  glMatrixMode(GL_MODELVIEW);
}


// Functions renders molecules in space. It is given molecules, their radius
// and also the repeat rate; this means how many times the same image is
// rendered.
void render(unsigned int n, Molecule* m, real radius, real a, unsigned int repeat) {
  for(unsigned int rx = 0; rx < repeat; rx++) {
    for(unsigned int ry = 0; ry < repeat; ry++) {
      // We setup the matrix.
      glLoadIdentity();
      
      // Translate to position.
      glTranslatef((float)rx*a/(real)repeat, (float)ry*a/(real)repeat, 0.0f);
      
      // We shrink by repeat rate to enable rendering image many times.
      glScalef(1.0f/((float)repeat), 1.0f/((float)repeat), 1.0f);
      glBegin(GL_POINTS);
      
      for(unsigned int i = 0; i < n; i++) {
        // We render all the molecules as points with their coordinates; hardware
        // takes care of scaling.
        glVertex2d(m[i].position.x, m[i].position.y);
      }
      glEnd();
    }
  }
}

// Globals used by main() and loop(). They are not accessible
// to any of the functions above directly.
unsigned int n;
real radius, a, maxr;
unsigned int M, S;

#ifdef SOFT_SIMULATION
real epsilon, sigma, T;
#endif

// Where does simulation run.
Molecule* m;

// The loop function.
void loop() {
  glClear(GL_COLOR_BUFFER_BIT);
  
  for(unsigned int k = 0; k < 100; k++) {
#ifndef SOFT_SIMULATION
  // We simulate the molecule.
    M += simulate(n, m, radius, a, maxr);
#else
    M += simulate(n, m, radius, a, maxr, epsilon, sigma, T);
#endif
    S++;
  }
  
  // we may need to reevaluate
  if(S > 1000) {
    maxr = reevaluateMaxR(maxr, M, S);
    maxr = maxr > 0.01 ? maxr : 0.01;
    printf("Max radius reevaluated to %f, %d, %d\n", maxr, M, S); 
    
    S = M = 0;
    
  }
  
  render(n, m, radius, a, 1);
  
  glFinish();
  glutSwapBuffers();
  
  // We need to trigger new rendering
  glutPostRedisplay();
  Sleep(100);
}

// Input handling function.
void input(unsigned char key, int x, int y) {
  /*if(key == '+') T *= 1.1;
  if(key == '-') T *= 0.9;
  printf("Temperature is %f\n", T);*/
}


// We enter main function.
int main(int argc, char** argv) {
    printf("Entering Monte Carlo simulation ...\n");
    ifstream in("parameters.in");
    
    // Randomize seed.
    srand(time(NULL));
    
    // We read input.
    in >> n >> radius >> a >> maxr;
#ifdef SOFT_SIMULATION
    in >> epsilon >> sigma >> T;
#endif
    
    // We initialize.
#ifndef SOFT_SIMULATION
    initSpace(n, a, m, radius);
#else
    initSpace(n, a, m, radius, epsilon, sigma, maxr);
#endif
    
    // Initialialize GLUT.
    glutInit(&argc, argv);
	glutInitDisplayMode(GLUT_RGBA | GLUT_DOUBLE);

	glutInitWindowPosition(0,0);
	glutInitWindowSize(800, 640);
	glutCreateWindow("Monte carlo simulacija molekul");

	glutDisplayFunc(loop);
	glutKeyboardFunc(input);

    // Setup color.
    glColor3f(1.0f, 1.0f, 1.0f);
    glPointSize(5.0);

    renderPrepare(a);
	glutMainLoop();    
	delete [] m;
    
}