#include <stdio.h>
#include <math.h>
#include <stdlib.h>

#define PNUM   200
#define Temp   1.0 
#define Dens   1.1
#define MAXLEN 10000
#define RMAX   3.3
#define RCUT   2.5
#define Dt     0.004
#define Tt     20.0
#define At     10.0

double pos[PNUM][2];
double vel[PNUM][2];
double force[PNUM][2];

double dist[MAXLEN][3];
int    pair[MAXLEN][2];
int    pairnum;

double a, virial;

double gauss_rand()
{
  double r;

  while(1){
    r = (double)rand()/RAND_MAX;
    if(exp(-(r-0.5)*(r-0.5)) >= (double)rand()/RAND_MAX)
      break;
  }
  return(r);
}

void print_lattice()
{
  int i;
  for (i = 0; i < PNUM; ++i)
    printf("%11.3le %11.3le %11.3le %11.3le\n", pos[i][0], pos[i][1], vel[i][0]*vel[i][0]+vel[i][1]*vel[i][1], force[i][0]*force[i][0]+force[i][1]*force[i][1]);
}

void init_lattice()
{
  int i, j, k, n;
  double vol, lconst;
  
  vol =  (double)PNUM/Dens;
  n = (int)sqrt((double)PNUM/2.0);
  a = sqrt(vol);
  lconst = sqrt(vol)/n;

  k = 0;
  for (i = 0; i < n ; i++){
    for (j = 0; j < n ; j++){
      pos[k][0] = (i+0.25)*lconst;
      pos[k][1] = (j+0.25)*lconst;
      k++;
      pos[k][0] = (i+0.75)*lconst;
      pos[k][1] = (j+0.75)*lconst;
      k++;
    }
  }
}

void init_mom()
{
  int i;
  double xsum = 0.0, ysum = 0.0;


  for (i = 0; i < PNUM; ++i){
    vel[i][0] = gauss_rand();
    vel[i][1] = gauss_rand();
  }

  for (i = 0; i < PNUM; ++i){
    xsum = xsum + vel[i][0];
    ysum = ysum + vel[i][1];
  }
  
  xsum = xsum / PNUM;
  ysum = ysum / PNUM;
  
  for (i = 0; i < PNUM; ++i){
    vel[i][0] = vel[i][0] - xsum;
    vel[i][1] = vel[i][1] - ysum;
  }
}

void rescale()
{
  int i;
  double exsum = 0, eysum = 0, ene, factor;
  for (i = 0; i < PNUM; ++i){
    exsum = exsum + vel[i][0]*vel[i][0];
    eysum = eysum + vel[i][1]*vel[i][1];
  }
  ene = exsum + eysum;
  factor = sqrt((PNUM-1)*2.0*Temp/ene);
  for (i = 0; i < PNUM; ++i){
    vel[i][0]=factor*vel[i][0];
    vel[i][1]=factor*vel[i][1];
  }
}

void calc_prop()
{
  int i;
  double pxsum = 0.0, pysum = 0.0, exsum = 0.0, eysum = 0.0; 

  for (i = 0; i < PNUM; ++i){
    pxsum = pxsum + vel[i][0];
    pysum = pysum + vel[i][1];
    exsum = exsum + vel[i][0]*vel[i][0];
    eysum = eysum + vel[i][1]*vel[i][1];
  }
}

void find_pair()
{
  int i, j;
  double dx, dy;

  pairnum = 0;
  for(i = 0; i < PNUM; ++i)
    for(j = i + 1; j < PNUM; ++j){
      dx = pos[i][0]-pos[j][0];
      dy = pos[i][1]-pos[j][1];
      if(dx >  a/2.0) dx = dx - a;
      if(dx < -a/2.0) dx = dx + a;
      if(dy >  a/2.0) dy = dy - a;
      if(dy < -a/2.0) dy = dy + a;

      if(dx*dx + dy*dy < RMAX*RMAX){
	pair[pairnum][0] = i;
	pair[pairnum][1] = j;
	pairnum++;
	if(pairnum >= MAXLEN){
	  printf("OUT OF MEMORY\n");
	  exit(0);
	}
      }
    }
}

void calc_pair()
{
  int i;
  double dx, dy;
  
  for(i = 0; i < pairnum; ++i){
    dx = pos[pair[i][0]][0] - pos[pair[i][1]][0];
    dy = pos[pair[i][0]][1] - pos[pair[i][1]][1];
    if(dx >  a/2.0) dx = dx - a;
    if(dx < -a/2.0) dx = dx + a;
    if(dy >  a/2.0) dy = dy - a;
    if(dy < -a/2.0) dy = dy + a;
    
    dist[i][0] = dx;
    dist[i][1] = dy;
    dist[i][2] = dx*dx + dy*dy;
  }
}

void calc_force()
{
  int i, j, pnum1, pnum2;
  double fc, r2, r4, r8, r14;

  for(i = 0; i < PNUM; ++i){
    force[i][0] = 0.0;
    force[i][1] = 0.0;
  }
  
  virial = 0.0;
  for(i = 0; i < pairnum; ++i){
    pnum1 = pair[i][0];
    pnum2 = pair[i][1];
    if(dist[i][2] <  RCUT*RCUT){
      r2 = 1.0 / dist[i][2];
      r4 = r2*r2;
      r8 = r4*r4;
      r14 = r8*r4*r2;
      fc = 12.0*(r14-r8);
      
      force[pnum1][0] = force[pnum1][0] + dist[i][0]*fc;
      force[pnum1][1] = force[pnum1][1] + dist[i][1]*fc;
      force[pnum2][0] = force[pnum2][0] - dist[i][0]*fc;
      force[pnum2][1] = force[pnum2][1] - dist[i][1]*fc;
      virial = virial + fc*dist[i][2];
    }
  }
}

void integration()
{
  int i;
  double dt2;
  
  dt2 = Dt*Dt/2.0;
  
  for(i = 0; i < PNUM; ++i){
    pos[i][0]  = pos[i][0] + Dt*vel[i][0] + dt2*force[i][0];
    pos[i][1]  = pos[i][1] + Dt*vel[i][1] + dt2*force[i][1];

    if(pos[i][0] > a) pos[i][0] = pos[i][0] - a;
    if(pos[i][0] < 0) pos[i][0] = pos[i][0] + a;
    if(pos[i][1] > a) pos[i][1] = pos[i][1] - a;
    if(pos[i][1] < 0) pos[i][1] = pos[i][1] + a;
  }
  calc_pair();
  for(i = 0; i < PNUM; ++i){
    vel[i][0]  = vel[i][0] + 0.5*Dt*force[i][0];
    vel[i][1]  = vel[i][1] + 0.5*Dt*force[i][1];
  }
  calc_force();
  for(i = 0; i < PNUM; ++i){
    vel[i][0]  = vel[i][0] + 0.5*Dt*force[i][0];
    vel[i][1]  = vel[i][1] + 0.5*Dt*force[i][1];
  }
}

double calc_temp()
{
  int i;
  double exsum = 0, eysum = 0, ene;
  for (i = 0; i < PNUM; ++i){
    exsum = exsum + vel[i][0]*vel[i][0];
    eysum = eysum + vel[i][1]*vel[i][1];
  }
  ene = (exsum + eysum)/2.0;
  return(ene/PNUM);
}

double calc_pot()
{
  int i;
  double pot, r2, r4, r6, r12;
  
  pot = 0.0;
  for(i = 0; i < pairnum; ++i){
    if(dist[i][2] < RCUT*RCUT){
      r2 = 1/dist[i][2];
      r4 = r2*r2;
      r6 = r4*r2;
      r12= r6*r6;
      pot = pot + (r12 - 2.0*r6);
    }     
  }
  return(pot/PNUM);
}

main()
{
  int i;

  srand(55);

  init_lattice();
  init_mom();
  rescale();
  find_pair();
  calc_pair();
  calc_force();
  
  for(i = 0; i < Tt/Dt; ++i){
    integration();
    if(i < (int)(At/Dt) && i%10 == 0) rescale();
    if(i%25 == 0) {
      find_pair(); 
      calc_pair(); 
      calc_force();
    }
    if(i%5 == 0){
      printf("%lf  %lf  %lf  %lf\n", Dt*i, calc_temp(), calc_pot(), calc_temp()+calc_pot());
      print_lattice();
    }
  }
}