#include <GL/glut.h>
#include <cstdlib>
#include <cmath>
#include <iostream>
#include <climits>
#include <ctime>
#include <fstream>

using namespace std;

typedef double real;

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

// A molecule.
struct Mol {
  Vec2 pos;
  Vec2 vel;
  Vec2 force;       
};

// Global data.
unsigned int n, simn;
Mol* molecules;
real eps, sig, size, mass, dt, T, limit, g;

// Calculates force based on distance
real calcForce(real radius, real epsilon, real sigma) {
  real force = 24.0*epsilon * (-2.0*pow(sigma/radius, 13.0)/sigma +  pow(sigma/radius, 7.0)/sigma);
  if(force > limit) return limit;
  if(force < -limit) return -limit;
  return force; 
}

// Finds the nearest radius squared.
real nearRadius2(const Vec2& p1, const Vec2& p2, const real& a) {
  Vec2 res;
  res.x = p2.x - p1.x;
  res.y = p2.y - p1.y;
  return res.x*res.x + res.y*res.y;
}

// We simulate.
void simulate(real dt) {
     
  for(unsigned int z = 0; z < n; z++) { molecules[z].force.x = 0.0; molecules[z].force.y = g; }
  // We calculate forces.
  for(unsigned int i = 0; i < n; i++) {
    for(unsigned int j = i+1; j < n; j++) {
      real rad = sqrt(nearRadius2(molecules[i].pos, molecules[j].pos, size));
      real force_scalar = calcForce(rad, eps, sig);
      
      // We setup force and normalize it in direction.
      Vec2 force;
      force.x = molecules[i].pos.x - molecules[j].pos.x;
      force.y = molecules[i].pos.y - molecules[j].pos.y;
      real l_inv = 1.0 / sqrt(force.x * force.x + force.y * force.y);   
      force.x *=  l_inv * force_scalar; 
      force.y *=  l_inv * force_scalar;  
      
      // We set the force.
      molecules[i].force.x -= force.x;
      molecules[i].force.y -= force.y;
      molecules[j].force.x += force.x;
      molecules[j].force.y += force.y;
    }                
  }
  
  // We calculate new v and displacement.
  for(unsigned int k = 0; k < n; k++) {
    molecules[k].vel.x += molecules[k].force.x / mass * dt;   
    molecules[k].vel.y += molecules[k].force.y / mass * dt;
    
    // New position.
    Vec2 newPos;
    newPos.x = molecules[k].pos.x + molecules[k].vel.x*dt; 
    newPos.y = molecules[k].pos.y + molecules[k].vel.y*dt;
    if(newPos.x < 0.0 && molecules[k].vel.x < 0.0) molecules[k].vel.x *= -1.0;
    if(newPos.x > size && molecules[k].vel.x > 0.0) molecules[k].vel.x *= -1.0;
    if(newPos.y < 0.0 && molecules[k].vel.y < 0.0) molecules[k].vel.y *= -1.0;
    if(newPos.y > size && molecules[k].vel.y > 0.0) molecules[k].vel.y *= -1.0;
    // New position, clampfed.
    molecules[k].pos = newPos;
  }
}

unsigned int index=0, molindex;
Vec2 pos[100] = { {0,0} };


void render() {
  glColor3f(1,0,0);
  
  glPointSize(4.0);
  glBegin(GL_POINTS);     
  for(unsigned int i = 0; i < n; i++) {
    glVertex2f(molecules[i].pos.x, molecules[i].pos.y);             
  }
     
  glEnd();
  glPointSize(1.0);
  glBegin(GL_LINE_STRIP);
  
  index++;
  if(index > 999) index = 0;
  pos[index].x = molecules[molindex].pos.x;
  pos[index].y = molecules[molindex].pos.y;
  
  for(unsigned int j = index; j < 100; j++) {
    glVertex2f(pos[j].x, pos[j].y);             
  }
  
  for(unsigned int k = 0; k< index; k++) {
    glVertex2f(pos[k].x, pos[k].y);              
  }
  
  glEnd();
     
}


void display() {
  glClear(GL_COLOR_BUFFER_BIT);
  
  for(unsigned int i = 0; i < simn; i++) simulate(dt);
  render();
  
  glutSwapBuffers();
  glutPostRedisplay();
}

void key(unsigned char key, int x, int y) {
  if(key == '1') {
    T*=1.1;
    for(unsigned int i = 0; i < n; i++) {
      molecules[i].vel.x *= 1.1;
      molecules[i].vel.y *= 1.1;
    }    
  }
  if(key == '2') {
    T*=0.9;
    for(unsigned int i = 0; i < n; i++) {
      molecules[i].vel.x *= 0.9;
      molecules[i].vel.y *= 0.9;
    }
  } 
  if(key == '3') molindex = rand()%n;  
}

inline real frand() { return (real)rand() / (real)RAND_MAX; }

int main(int argc, char **argv)
{
    ifstream in("parameters.in");
    srand(time(NULL));
    in >> n >> eps >> sig >> size >> mass >> dt >> limit >> simn >> g;
    T = 300;
    molecules = new Mol[n];
    
    // Select random positions.
    for(unsigned int i = 0; i < n; i++) {
      molecules[i].force.x = molecules[i].force.y = 0.0;       
      molecules[i].vel.x = (frand()*2-1.0)*size;
      molecules[i].vel.y = (frand()*2-1.0)*size;
      molecules[i].pos.x = frand()*size; 
      molecules[i].pos.y = frand()*size;        
    }
    
    glutInit(&argc, argv);
    glutInitWindowSize(800,600);
    glutInitWindowPosition(10,10);
    glutInitDisplayMode(GLUT_RGB | GLUT_DOUBLE);

    glutCreateWindow("Moleculer dynamics");

    glMatrixMode(GL_PROJECTION);
    glTranslatef(-1,-1,0);
    glScalef(2,2,0);
    glMatrixMode(GL_MODELVIEW);
    glScalef(1.0/(real)size, 1.0/(real)size, 1.0);


    glPointSize(3);
    glutDisplayFunc(display);
    glutKeyboardFunc(key);

    glutMainLoop();

    delete [] molecules;
    return 0;
}