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

typedef struct {
  double m; 
  double x, y;
} particle;

typedef struct {
  int N;
  particle *p;
  double LX, LY;
  double delta; // ampiezza spostamento
  double gravity;
  double temp; // temperature
} simulation;

  
simulation * new_simulation(int N, double LX, double LY, double temp, double gravity, double delta);
void evolve(simulation *s);
void random_init(simulation *s); 
void visualize(FILE * gp, simulation * s, char * title);
double compute_energy(simulation *s);
int save(simulation *s, char * filename);
simulation * load(char * filename);
void destroy_simulation(simulation * s);
double compute_energy_particle(simulation *s, int k);



int main() {
  int t, T;
  simulation *s;
  FILE * gp; 
  double temp;
  int VIS;
  double delta;
  double LX, LY;
  int i, N;
  char str[100];
  char filename[100];
  double deltaT;
  char title[100];
  clock_t start, end;
  double cpu_time_used;
  double gravity;
 
  
  T = 100000;
  VIS = 1000;
  temp = 4; // temperature
  delta = 0.1;
  LX=20;
  LY=20;
  N=100;
  deltaT = 0;
  gravity = 0;
  s=NULL;
  
  
  gp = popen("gnuplot", "w");
  fprintf(gp, "set size square; unset key; set xrange [0:%f]; set yrange [0:%f]\n", LX, LY);
  // initial conditions
  srand48(time(NULL));
  
  
  for (;;) {
    printf("\nMake your choice (character [+ parameters]):\n");
    printf("\nDefaults: N=%d LX=%f LY=%f delta=%f \n", N, LX, LY, delta);
    printf("deltaT=%g gravity=%g temp=%f T=%d VIS=%d\n", deltaT,gravity,temp, T, VIS); 
    printf("?: Help\n");
    printf("I: Random init (N LX LY delta temp)\n");
    printf("L: Load simulation (file)\n");
    printf("S: Save simulation (file)\n");
    printf("R: Run simulation (steps)\n");
    printf("V: Set visualization intervals (steps)\n");
    printf("T: Set temperature\n");
    printf("G: Set gravity (g)\n");
    printf("Q: Quit\n");
    printf("\n> ");
      fgets(str, 100, stdin); 
      str[strcspn(str, "\n")]='\0'; // remove trailing newlines
      i = 0;
      while(isblank(str[i])) i++; // skip initial blanks
      switch(toupper(str[i++])) {
        case '?':
          printf("help is missing\n");
          break;
        case 'I':
          sscanf(str+i, "%d %lf %lf %lf %lf", &N, &LX, &LY, &delta, &temp);
          printf("Random init N=%d LX=%f LY=%f delta=%f temp=%f\n", N, LX, LY, delta,temp);
          if (s) destroy_simulation(s);
          s = new_simulation(N, LX, LY, temp, gravity, delta);
          fprintf(gp, "set size square; unset key; set xrange [0:%f]; set yrange [0:%f]\n", LX, LY);
          random_init(s);          
          break;
        case 'L':
          while(isblank(str[i])) i++; // skip  blanks
          sscanf(str+i, "%s", filename); 
          printf("Loading from file '%s'\n", filename);
          s = load(filename);
          printf("%d particles, LX=%lf, LY=%lf, delta=%lf\n", s->N, s->LX, s->LY, s->delta); 
          fprintf(gp, "set size square; unset key; set xrange [0:%f]; set yrange [0:%f]\n", s->LX, s->LY);
          break;
        case 'S':
          if (!s) {
            printf("Simulation is not initialized\n");
            break;
          }
          while(isblank(str[i])) i++; // skip  blanks
          sscanf(str+i, "%s", filename); 
          printf("Saving to file '%s'\n", filename);
          save(s,filename);
          break;
        case 'R':
          if (!s) {
            printf("Simulation is not initialized\n");
            break;
          }
          sscanf(str+i, "%d", &T);
          T = abs(T);
          printf("Running for %d steps\n", T);
          start = clock();
          for (t=0; t<T; t++) {
            evolve(s); // evolution
            //printf("%d) en=%f\n",t, compute_energy(s));
            if (t % VIS == 0) {  // visualization
              sprintf(title, "time=%d temp=%f", t, s->temp);
              visualize(gp, s, title);
            }
          }
          end = clock();
          cpu_time_used = ((double) (end - start)) / CLOCKS_PER_SEC;
          printf("CPU time used: %lf s\n", cpu_time_used);
          break;
        case 'V':
          sscanf(str+i, "%d", &VIS);
          VIS = abs(VIS);
          printf("Setting VIS to %d steps\n", VIS);
          break;
        case 'T':
          if (!s) {
            printf("Simulation is not initialized\n");
            break;
          }
          sscanf(str+i, "%lf", &temp);
          s->temp = temp;
          printf("Setting temperature=%f\n", temp);
           break;
        case 'G':
          sscanf(str+i, "%lf",  &gravity);
          gravity = fabs(gravity);
          if (gravity > 1) {
            printf("|gravity| should be less than 1\n");
            break;
          }
          s->gravity = gravity;
          printf("Setting gravity to %f\n", gravity);
          break;
        case 'Q':
          printf("Are you sure [Y/N]\n");
          fgets(str, 100, stdin);
          i=0;
          while(isblank(str[i])) i++; // skip initial blanks
          if (toupper(str[i])=='Y') exit(0);
          break;
        default:
          printf("Choice not present\n");
          break;
    }
  }
}

void visualize(FILE * gp, simulation * s, char * title) {
  int i;
  
  fprintf(gp, "set title '%s'\n", title);
  fprintf(gp, "plot '-' binary record=%d format='%%double' w p pt 6 ps %f \n", s->N, 50/s->LX);//##
  for (i=0; i<s->N; i++) {
    fwrite(&(s->p[i].x), sizeof(double),2, gp);
   // fwrite(&(p[i].y), sizeof(double),1, gp);
  }
  fflush (gp);
}


double   thre = 1; // minimal distance for initial conf

void random_init(simulation *s) {
  int i,j;
    
  for (i=0; i < s->N; i++) {
    double x, y, r2;
    int flag;
    
    s->p[i].m = 1;
    
    do {
      flag = 0;
      x = s->LX * drand48();
      y = s->LY * drand48();
      for (j=0; j<i; j++) {
        r2 = (x - s->p[j].x)*(x - s->p[j].x) + (y - s->p[j].y)*(y - s->p[j].y);
        if (r2 < thre) {
          flag = 1;
          break;
        }
      }
    } while (flag == 1);
    s->p[i].x = x;
    s->p[i].y = y;
  }
}

simulation * new_simulation(int N, double LX, double LY, double temp, double gravity, double delta) {
  simulation *s;
  
  s=calloc(1, sizeof(simulation));
  s->N = N;
  s->LX = LX;
  s->LY = LY;
  s->gravity = gravity;
  s->delta = delta;
  s->temp = temp;
  s->p = calloc(N, sizeof(particle));
  return (s);
}

void destroy_simulation(simulation *s) {
  free(s->p);
  free(s);
}

double compute_energy_particle(simulation *s, int k) {
  int i;
  double energy;
  double rx, ry, r2, rm2, rm6;
  
  energy = 0;
  for (i=0; i < s->N; i++) {    // loop su tutte le particelle diverse da k
    if (i != k) {
      rx = s->p[k].x - s->p[i].x; 
      ry = s->p[k].y - s->p[i].y; 
      r2 = rx*rx+ry*ry;
      rm2 = 1/r2;
      rm6 = rm2*rm2*rm2;
      energy += 4*rm6*(rm6-1);
    }
  }
  energy += s->p[k].m * s->gravity * s->p[k].y; // gravità
  return(energy);
}

double compute_energy(simulation *s) {
  int i, j;
  double energy;
  double rx, ry, r2, rm2, rm6;
  
  energy = 0;
  for (i=0; i<s->N-1; i++) {    // loop su metà delle particelle
    for (j=i+1; j<s->N; j++) {
      rx = s->p[j].x - s->p[i].x; 
      ry = s->p[j].y - s->p[i].y; 
      r2 = rx*rx+ry*ry;
      rm2 = 1/r2;
      rm6 = rm2*rm2*rm2;
      energy += 4*rm6*(rm6-1);
    }
  }
  return(energy);
}
   

void evolve(simulation *s) {
  int i, k;
  double energy;
  double oldx, oldy, oldenergy, deltaE;
  

  for (i=0; i < s->N; i++) { // 1 passo monte carlo
    k = (int) (drand48()* s->N); // una particella a caso
    oldenergy=compute_energy_particle(s, k);
    oldx = s->p[k].x;
    oldy = s->p[k].y;
    s->p[k].x += s->delta*(2*drand48() -1);
    s->p[k].y += s->delta*(2*drand48() -1);
    // condizioni al bordo
    if (s->p[k].x < 0) {
      s->p[k].x *= -1;
    }
    if (s->p[k].y < 0) {
      s->p[k].y *= -1;
    } 
    if (s->p[k].x > s->LX) {
      s->p[k].x = 2 * s->LX - s->p[k].x;
    }
    if (s->p[k].y > s->LX) {
      s->p[k].y = 2 * s->LY - s->p[k].y;
    } 
    energy = compute_energy_particle(s, k);
    deltaE = energy-oldenergy;
    if (deltaE <= 0|| drand48() < exp(-deltaE/s->temp)) {
      // accept
    } else {
      //revert
      s->p[k].x = oldx;
      s->p[k].y = oldy;
    }
  }
}

int save(simulation *s, char * filename) {
  int i; 
  FILE * f;
  
  if ((f = fopen(filename, "w")) == NULL) {
    printf("file '%s' is not writable\n", filename);
    return (1);
  }
  fprintf(f, "%d %f %f %f %f %f\n",s->N, s->LX, s->LY, s->temp, s->gravity, s->delta);
  for (i=0; i< s->N; i++) {
    fprintf(f, "%f %f %f\n", 
      s->p[i].m,s->p[i].x,s->p[i].y);
  }
  fclose(f);
  return(0);
}
  
simulation * load(char * filename) {
  simulation *s;
  int i, N;
  double LX, LY, gravity, delta, temp;
  double m, x, y;
  FILE *f;
  
  if ((f = fopen(filename, "r")) == NULL) {
    printf("file '%s' does not exists or is not readable\n", filename);
    return (NULL);
  }
  fscanf(f, "%d %lf %lf %lf %lf %lf\n", &N, &LX, &LY, &temp, &gravity, &delta);
  s = new_simulation(N, LX, LY, temp, gravity, delta);
  for (i=0; i<N; i++) {
    fscanf(f, "%lf %lf %lf\n", &m, &x, &y);
    s->p[i].m=m;
    s->p[i].x=x;
    s->p[i].y=y;
   }
  return(s);
}