#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;
  double vx, vy;
  double fx, fy;
} particle;

typedef struct {
  int N;
  particle *p;
  double LX, LY;
  double deltat, deltat2;
  double gravity;
} simulation;

  
simulation * new_simulation(int N, double LX, double LY, double gravity, double deltat);
void evolve(simulation *s, double *energy, double *K);
void random_init(simulation *s, double temp); 
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);
void increment_temperature(simulation * s, double deltaT);



int main() {
  int t, T;
  simulation *s;
  FILE * gp; 
  double temp;
  int VIS;
  double energy, K;
  double deltat;
  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 = 1000000;
  VIS = 1000;
  temp = 4; // temperature
  deltat = 0.0001;
  LX=100;
  LY=100;
  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 deltat=%f \n", N, LX, LY, deltat);
    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 deltat 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 change (deltaT)\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, &deltat, &temp);
          printf("Random init N=%d LX=%f LY=%f delta=%f temp=%f\n", N, LX, LY, deltat,temp);
          if (s) destroy_simulation(s);
          s = new_simulation(N, LX, LY, gravity, deltat);
          random_init(s, temp);          
          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, deltat=%lf\n", s->N, s->LX, s->LY, s->deltat); 
          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++) {
            increment_temperature(s, deltaT);
            evolve(s, &energy, &K); // evolution
            //printf("%d) en=%f\n",t, compute_energy(s));
            if (t % VIS == 0) {  // visualization
              sprintf(title, "time=%d, energy=%f, K = %f", t, energy, K);
              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':
          sscanf(str+i, "%lf", &deltaT);
          if (deltaT < -1e-2 || deltaT>1e-2) {
            printf("|deltaT| should be less than 1E-2\n");
            break;
          }
          printf("Setting relative change of kinetic energy deltaT=%f\n", deltaT);
           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 \n", s->N);
  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 = 0.2; // minimal distance for initial conf

void random_init(simulation *s, double temp) {
  int i,j;
    
  for (i=0; i < s->N; i++) {
    double x, y, r2;
    int flag;
    
    s->p[i].m = 1;
    s->p[i].vx = temp * (2* drand48() - 1);
    s->p[i].vy = temp * (2* drand48() - 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 gravity, double deltat) {
  simulation *s;
  
  s=calloc(1, sizeof(simulation));
  s->N = N;
  s->LX = LX;
  s->LY = LY;
  s->gravity = gravity;
  s->deltat = deltat;
  s->deltat2=deltat*deltat;
  s->p = calloc(N, sizeof(particle));
  return (s);
}

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

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);
    }
    energy += 0.5*s->p[i].m * (s->p[i].vx*s->p[i].vx + s->p[i].vy*s->p[i].vy);
  }
  return(energy);
}
   

void evolve(simulation *s, double *energy, double *K) {
    int i, j;
    double rx, ry, r2, rm2, rm6, ff;

    *energy = *K = 0;
    //avanzo le posizioni
    for (i=0; i<s->N; i++) {
      s->p[i].x += s->p[i].vx * s->deltat + 0.5*s->p[i].fx * s->deltat2 / s->p[i].m;
      s->p[i].y += s->p[i].vy * s->deltat + 0.5*s->p[i].fy * s->deltat2 / s->p[i].m;
      // condizioni al bordo
      if (s->p[i].x < 0) {
        s->p[i].x *= -1;
        s->p[i].vx *= -1;
      }
      if (s->p[i].y < 0) {
        s->p[i].y *= -1;
        s->p[i].vy *= -1;
      } 
      if (s->p[i].x > s->LX) {
        s->p[i].x = 2 * s->LX - s->p[i].x;
        s->p[i].vx *= -1;
      }
      if (s->p[i].y > s->LX) {
        s->p[i].y = 2 * s->LY - s->p[i].y;
        s->p[i].vy *= -1;
      } 
   }
    // prima metà della velocità  
    for (i=0; i<s->N; i++) {
      s->p[i].vx += 0.5 * s->p[i].fx*s->deltat/s->p[i].m;
      s->p[i].vy += 0.5 * s->p[i].fy*s->deltat/s->p[i].m;
      s->p[i].fx = 0;
      s->p[i].fy = 0;
    }
    
    //calcolo delle forze
    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;
        ff = -24*rm6*(2*rm6-1)*rm2;
        s->p[i].fx += ff * rx;
        s->p[i].fy += ff * ry;
        s->p[j].fx -= ff * rx;
        s->p[j].fy -= ff * ry;
        *energy += 4*rm6*(rm6-1);
      }
    }
    // gravity 
    if (s->gravity > 0) {
      for (i=0; i< s->N; i++) {
        s->p[i].fy -= s->gravity * s->p[i].m;
      }
    }
    //seconda metà della velocità 
    for (i=0; i<s->N; i++) {
      s->p[i].vx += 0.5 * s->p[i].fx*s->deltat/s->p[i].m;
      s->p[i].vy += 0.5 * s->p[i].fy*s->deltat/s->p[i].m;
      *K += 0.5*(s->p[i].vx*s->p[i].vx + s->p[i].vy*s->p[i].vy)*s->p[i].m;
    }
    *energy += *K;
    *energy /= s->N; // energia per particella
    *K /= s->N;
}


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\n",s->N, s->LX, s->LY, s->gravity, s->deltat);
  for (i=0; i< s->N; i++) {
    fprintf(f, "%f %f %f %f %f %f %f\n", 
      s->p[i].m,s->p[i].x,s->p[i].y,s->p[i].vx,s->p[i].vy,s->p[i].fx,s->p[i].fy);
  }
  fclose(f);
  return(0);
}
  
simulation * load(char * filename) {
  simulation *s;
  int i, N;
  double LX, LY, gravity, deltat;
  double m, x, y, vx, vy, fx, fy;
  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\n", &N, &LX, &LY, &gravity, &deltat);
  s = new_simulation(N, LX, LY, gravity, deltat);
  for (i=0; i<N; i++) {
    fscanf(f, "%lf %lf %lf %lf %lf %lf %lf\n", &m, &x, &y, &vx, &vy, &fx, &fy);
    s->p[i].m=m;
    s->p[i].x=x;
    s->p[i].y=y;
    s->p[i].vx=vx;
    s->p[i].vy=vy;
    s->p[i].fx=fx;
    s->p[i].fy=fy;
  }
  return(s);
}
    
void increment_temperature(simulation * s, double deltaT) {
  int i;
  
  if (deltaT == 0.0) {
    return;
  }
  for (i=0; i< s->N; i++) {
    s->p[i].vx *= 1+deltaT;
    s->p[i].vy *= 1+deltaT;
  }
}