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

using namespace std;

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

// Number of particles in one dimensions and in two.
unsigned int n, n2;
float delta;
double T, E, J, B, m;
ofstream out("out.txt");

unsigned int displayCount;

// Charge (either 1 or -1).
double*     particles;


// Prints string.
void glutPrint(float x, float y, char* c) {
  glColor3f(1,1,1);
  glRasterPos2f(x,y);
  for(;*c;c++) {
    glutBitmapCharacter(GLUT_BITMAP_TIMES_ROMAN_24, *c);            
  }     
}

double calcE(unsigned int x, unsigned int y, double pe) {
  double r = 0.0f;
  if(x > 0) r += -J*particles[x-1+y*n]*pe - B*pe;
  if(x < (n-1)) r += -J*particles[x+1+y*n]*pe - B*pe; 
  if(y < (n-1)) r += -J*particles[x+(y+1)*n]*pe - B*pe; 
  if(y > 0) r += -J*particles[x+(y-1)*n]*pe - B*pe; 
  return r;
}

void simulate() {
  unsigned int i = rand()%n2;
  unsigned x = i/n, y = i%n;
  
  // We now calculate E
  double pe = particles[x+y*n];
  double zE = calcE(x,y, pe);
  
  // We switch and recalculate E
  double kE = calcE(x,y, pe*-1.0);
  
  // We check dE
  double dE = kE - zE;
  if(dE > 0.0) {
    double n = exp(-dE*100/T);
    if(frand() > n) return; 
  }
  
  particles[x+y*n] *= -1.0;    
}

void recalcE() {
  E = 0.0;
  for(unsigned int x = 0; x < n; x++) {
    for(unsigned int y = 0; y < n; y++) {
      double pe = particles[x+y*n];
      E += calcE(x,y,pe);           
    }             
  }
}

double calcMag() {
  double r = 0.0f;
  for(unsigned int i = 0; i < n2; i++) {
    r += particles[i];             
  }       
  return r;
}

void display() {
  glClear(GL_COLOR_BUFFER_BIT);  
  
  simulate();
  
  // We render particles
  glBegin(GL_POINTS);
  unsigned int x = 0, y = 0;
  for(float fx = delta/2; fx < 1.0f; fx += delta, x++) {
    y = 0;
    for(float fy = delta/2; fy < 1.0f; fy += delta, y++) {
      if(particles[x+y*n] > 0.0) glColor3f(1.0f, 0.0f, 0.0f);
      else glColor3f(0.0f, 1.0f, 0.0f);
      glVertex2f(fx, fy);
    }                     
  }
  glEnd();
  
  if(++displayCount > 100) {
    // recalc E
    recalcE();
    
    m = calcMag();
    out << T << " " << B << " " << m << "\n";
    
    displayCount = 0;
  } 
  
  // Output T
  char _ttext[100];
  sprintf(_ttext, "T = %f xK", T);
  glutPrint(0.01, 0.01, _ttext);
  
  // Output E
  char _etext[100];
  sprintf(_etext, "E = %f xJ", E);
  glutPrint(0.01, 0.07, _etext); 
  
  // Output B
  char _btext[100];
  sprintf(_btext, "B = %f xT", B);
  glutPrint(0.01, 0.13, _btext); 
  
  // Output J
  char _jtext[100];
  sprintf(_jtext, "J = %f xJ", J);
  glutPrint(0.01, 0.19, _jtext); 
  
  // Output J
  char _emtext[100];
  sprintf(_emtext, "e = %f", m);
  glutPrint(0.01, 0.25, _emtext); 
  
  glFinish();
  glutSwapBuffers();
  glutPostRedisplay();
}

void key(unsigned char key, int x, int y) {
  if(key == '+') T++;
  if(key == '-') T--;
  if(key == '1') B-=0.01;
  if(key == '2') B+=0.01;
  if(key == '3') J*=-1;
}
     
double randtable[] = { -1.0, 1.0 };

int main(int argc, char **argv) {
    ifstream in("parameters.in");
    in >> n >> T >> J >> B;
    
    // generate data
    n2 = n*n;
    delta = 1.0f / (float)n;
    particles = new double[n2];
    
    // initialize data
    for(unsigned int i = 0; i < n2; i++) {
      particles[i] = randtable[rand()&0x1];             
    }
    
    glutInit(&argc, argv);
    glutInitWindowSize(640,480);
    glutInitWindowPosition(10,10);
    glutInitDisplayMode(GLUT_RGB | GLUT_DOUBLE);

    glutCreateWindow("Magnetism - Monte Carlo approach");
    
    // Setups GL
    glPointSize(2.0f);
    glMatrixMode(GL_PROJECTION);
    glLoadIdentity();
    glTranslatef(-1.0f, -1.0f, 0.0f);
    glScalef(2.0f, 2.0f, 1.0f);
    glMatrixMode(GL_MODELVIEW);
    
    glutDisplayFunc(display);
    glutKeyboardFunc(key);

    glutMainLoop();

    delete [] particles;
    return 0;
}