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

typedef struct {
  double m; 
  double x, y;
  double vx, vy;
  double fx, fy;
} particle;


int main() {
  int N; // numero di particelle
  particle * p;
  int t, T;
  int i, j;
  double deltat, deltat2;
  double LX, LY;
  double thre;
  FILE * gp; 
  double temp;
  int VIS;
  
  N = 100;
  T = 1000000;
  VIS = 100;
  deltat = 0.01;
  deltat2=deltat*deltat; 
  LX = 100; 
  LY = 100;
  thre = 1e-2;
  temp = 0.01;
  
  p = calloc(N, sizeof(particle));
  
  gp = popen("gnuplot", "w");
  fprintf(gp, "set xrange [0:%f]; set yrange [0:%f]\n", LX, LY);
  // initial conditions
  srand48(time(NULL));
  
  for (i=0; i<N; i++) {
    double x, y, r2;
    int flag;
    
    p[i].m = 1;
    p[i].vx = temp * (2* drand48() - 1);
    p[i].vy = temp * (2* drand48() - 1);
    
    do {
      flag = 0;
      x = LX * drand48();
      y = LY * drand48();
      for (j=0; j<i; j++) {
        r2 = (x-p[j].x)*(x-p[j].x) + (y-p[j].y)*(y-p[j].y);
        if (r2 < thre) {
          flag = 1;
          break;
        }
      }
    } while (flag == 1);
    p[i].x = x;
    p[i].y = y;
  }
    
  
  // main loop;
  for (t=0; t<T; t++) {
    double rx, ry, r2, rm2,rm6,ff;
    double energy, K;
   
    // va in una funzione
    
    energy = K = 0;
    //avanzo le posizioni
    for (i=0; i<N; i++) {
      p[i].x += p[i].vx * deltat + 0.5*p[i].fx*deltat2/p[i].m;
      p[i].y += p[i].vy * deltat + 0.5*p[i].fy*deltat2/p[i].m;
      // condizioni al bordo
      if (p[i].x < 0) {
        p[i].x *= -1;
        p[i].vx *= -1;
      }
      if (p[i].y < 0) {
        p[i].y *= -1;
        p[i].vy *= -1;
      } 
      if (p[i].x > LX) {
        p[i].x = 2 * LX - p[i].x;
        p[i].vx *= -1;
      }
      if (p[i].y > LX) {
        p[i].y = 2 * LY - p[i].y;
        p[i].vy *= -1;
      } 
   }
    // prima metà della velocità  
    for (i=0; i<N; i++) {
      p[i].vx += 0.5 * p[i].fx*deltat/p[i].m;
      p[i].vy += 0.5 * p[i].fy*deltat/p[i].m;
      p[i].fx = 0;
      p[i].fy = 0;
    }
    
    //calcolo delle forze
    for (i=0; i<N-1; i++) {    // loop su metà delle particelle
      for (j=i+1; j<N; j++) {
        rx = p[j].x - p[i].x; 
        ry = p[j].y - p[i].y; 
        r2 = rx*rx+ry*ry;
        rm2 = 1/r2;
        rm6 = rm2*rm2*rm2;
        ff = -24*rm6*(2*rm6-1)*rm2;
        p[i].fx += ff * rx;
        p[i].fy += ff * ry;
        p[j].fx -= ff * rx;
        p[j].fy -= ff * ry;
        energy += 4*rm6*(rm6-1);
      }
    }
    //seconda metà della velocità 
    for (i=0; i<N; i++) {
      p[i].vx += 0.5 * p[i].fx*deltat/p[i].m;
      p[i].vy += 0.5 * p[i].fy*deltat/p[i].m;
      K += 0.5*(p[i].vx*p[i].vx+p[i].vy*p[i].vy)*p[i].m;
    }
    energy += K;
    energy /= N; // energia per particella
    K /= N;
    printf("energy=%f, K = %f\n", energy, K);
    
    
    // visualizzazione
    if (t % VIS == 0) {
      fprintf(gp, "plot '-' w p pt 6 \n");
      for (i=0; i<N; i++) {
        fprintf(gp, "%f %f\n", p[i].x, p[i].y);
      }
      fprintf(gp, "e\n");
      fflush (gp);
    }
  }
 
  
}