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

// handling control-backslash
#include <signal.h>
static int keepRunning = 1;
void intHandler(int sig) {
    keepRunning = 0;
}

#ifdef READLINE
#include <readline/readline.h>
#include <readline/history.h>
#endif

#define noDEBUG
#define noVISDEBUG

// linked lists


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

typedef struct {
  int N;
  particle *p;
  double LX, LY;
  double deltat, deltat2;
  double gravity;
  double cell_size;
  int * cell_count;
  int * cell_ptr;
  int Mx, My;
  double rc;  // truncation of potential, normally 2.5
  double rc2; // rc^2
  double shift; // V(rc)
} simulation;


  
simulation * new_simulation(int N, double LX, double LY, double gravity,
    double deltat, double cell_size, double rc);
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);
void init_visualize(FILE *gp, simulation * s);

typedef struct {
  int n;
  int nbin;
  double max, min;
  int * c;
} histo;

typedef struct  {
  FILE * gp;
  histo * h;
  char * descr;
} vhisto;

histo * new_histo(int nbin, double min, double max);
void destroy_histo (histo * h);
void clear_histo(histo * h);
int add_histo(histo * h, double x);
vhisto * new_vhisto(int nbin, double min, double max, char * descr);
void destroy_vhisto(vhisto * v);
void add_distr(vhisto * v, vhisto * y, simulation *s);
void visualize_distr(vhisto * v);


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;
  double cell_size;
  double rc;
  vhisto * vdistr, * ydistr;
  
  T = 1000000;
  VIS = 1000;
  temp = 1; // temperature
  deltat = 0.0001;
  LX=100;
  LY=100;
  cell_size = 2.5; rc = 2.5;
  N=100;
  deltaT = 0;
  gravity = 0;
  s=NULL;
  vdistr = NULL;
  ydistr = NULL;
  
  // handling control-backslash
  struct sigaction act;
  act.sa_handler = intHandler;
  sigaction(SIGQUIT, &act, NULL);
  
  gp = popen("gnuplot", "w");
   // 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("C: Clear histograms\n");
    printf("Q: Quit\n");
    printf("\n");
#ifdef READLINE
    char * tmp;
    tmp = readline("# ");
    add_history(tmp);
    strcpy(str,tmp);
    free(tmp);
#else
    printf("> ");
    fgets(str, 100, stdin);
    str[strcspn(str, "\n")]='\0'; // remove trailing newlines
#endif
    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, cell_size, rc);
          
          if (ydistr) destroy_vhisto(ydistr);
          ydistr = new_vhisto(20, 0, s->LY, "vertical distribution");
          if (vdistr) destroy_vhisto(vdistr);
          vdistr = new_vhisto(20, 0, 4, "velocity distribution");
         
          init_visualize(gp, s);
          random_init(s, temp);    
          visualize(gp, s, "init");      
          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);
          
          if (ydistr) destroy_vhisto(ydistr);
          ydistr = new_vhisto(20, 0, s->LY, "vertical distribution");
          if (vdistr) destroy_vhisto(vdistr);
          vdistr = new_vhisto(50, 0, 4, "velocity distribution");

          init_visualize(gp, s);
          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();
          keepRunning = 1;
          for (t=0; t<T && keepRunning; 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);
              add_distr(vdistr, ydistr, s);
              visualize_distr(vdistr);
              visualize_distr(ydistr);
            }
          }
          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 'C':
              clear_histo(vdistr->h);
              clear_histo(ydistr->h);
            printf("Cleared\n");
            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;
    }
  }
}


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, double cell_size, double rc) {
  simulation *s;
  
  s=calloc(1, sizeof(simulation));
  s->N = N;
  s->LX = LX;
  s->LY = LY;
  s->deltat = deltat;
  s->deltat2=deltat*deltat;
  s->p = calloc(N, sizeof(particle));
  s->Mx = ((int) (LX/cell_size)) + 1;  // for half-cells
  s->My = ((int) (LY/cell_size)) + 1;  
  s->cell_count=calloc(s->Mx * s->My + 1, sizeof(int)); // vedi evolve
  s->cell_ptr=calloc(N, sizeof(int));
  s->cell_size = cell_size;
  s->rc = rc;
  s->rc2 = s->rc*s->rc;
  s->shift = 4*pow(rc, -6)*(pow(rc, -6)-1);
  printf("rc=%f cell_size=%f Mx=%d My=%d, shift=%g\n", s->rc, s->cell_size, s->Mx, s->My, s->shift);

  return (s);
}

void destroy_simulation(simulation *s) {
  free(s->p);
  free(s->cell_count);
  free(s->cell_ptr);
  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) - s->shift;
    }
    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);
}
     

int Delta[][2]= {{0,0}, {1,0}, {1,1}, {0,1}, {-1,1}}; // celle vicine


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;
      } 
    }
    
    // azzero la lista delle celle
    for (i=1; i< s->My * s->Mx; i++) {
      s->cell_count[i] = 0;
    }
    // calcolo quante particelle stanno nelle celle 
    for (i=0; i< s->N; i++) {
      int ix, iy;
      ix = (int) (s->p[i].x / s->cell_size);
      iy = (int) (s->p[i].y / s->cell_size);
      s->p[i].cell=iy * s->Mx + ix; // memorizzo la cella nella struttura
      s->cell_count[s->p[i].cell] ++; // incremento il contatore
    }
    // a questo punto cell_count contiene per es. [5,3,3,4...]
    // supponiamo sia
    // cella    particelle
    //  0       2, 4, 6, 11, 12 (5 particelle)  
    //  1       1, 7, 10        (3 particelle)  
    //  2       5, 9, 14        (4 particelle)  
    //  3       0, 3, 8, 13     (4 particelle)  
    //  ...  
 
#ifdef DEBUG
printf("cell_count:\n");
for (i=0; i< s->My * s->Mx ; i++) {
  printf("cell_count[%d]=%d\n", i, s->cell_count[i]);
}
#endif
  
    // calcolo la cumulata
    for (i=1; i< s->My * s->Mx ; i++) {
      s->cell_count[i] += s->cell_count[i-1];
    }
    // a questo punto cell_count contiene [5,8,11,15...] 
   
#ifdef DEBUG
printf("\ncumulata:\n");
for (i=0; i< s->My * s->Mx ; i++) {
  printf("cell_count[%d]=%d\n", i, s->cell_count[i]);
}
#endif

    // inserisco i puntatori alle particelle nella lista
    // ma prima decremento il contatore, così va nella posizione giusta
    for (i=0; i< s->N; i++) {
      s->cell_ptr[-- s->cell_count[s->p[i].cell]] = i;
    }
    // a questo punto cell_count contiene [0,5,8,11...] 
    // e inserisco nell'ultima posizione N
    s->cell_count[s->My * s->Mx]=s->N;

    // cell_ptr è 
    // indice    0  | 1  | 2 | 3 | 4 | 5  | 6 | 7 | 8  | 9 | 10 | 11 | 12 | 13 | 14 
    // partic.   12 | 11 | 6 | 4 | 2 | 10 | 7 | 1 | 14 | 9 | 5  | 13 | 8  | 3  | 0  
    // (cella)   0  | 0  | 0 | 0 | 0 | 1  | 1 | 1 | 2  | 2 | 2  | 3  | 3  | 3  | 3 
    
    // e cell_count è tornato ad essere [0,5,8,11...] 
    // ogni numero è l'indice di inizio e la differenza con 
    // il successivo dà l'indice di fine (più uno)

    
#ifdef DEBUG
printf("\ncell_ptr:\n");
for (i=0; i< s->N ; i++) {
  printf("cell_ptr[%d]=%d\n", i, s->cell_ptr[i]);
}

printf("\ncell_count:\n");
for (i=0; i< s->My * s->Mx ; i++) {
  printf("cell_count[%d]=%d\n", i, s->cell_count[i]);
}
#endif

    //calcolo delle forze
    // tanto vale procedere per celle
    int ii, ix, iy, jx, jy, ki, kj, ci, cj; 

#ifdef DEBUG
printf("\nneighbors:\n");
for (i=0; i<s->N; i++) {
  int ix, iy, jx, jy, k, ci, cj;  
  ci = s->p[i].cell; 
  ix = ci % s->Mx; // recupero gli indici x, y della cella
  iy = ci / s->Mx;
  printf("part[%d] (%d: %d %d): ", i, ci, ix, iy); 
  for (jx=ix-1; jx<=ix+1; jx++) { // loop sulle celle adiacienti
    for (jy=iy-1; jy<=iy+1; jy++) {
      if (jx >=0 && jx < s->Mx && jy >=0 && jy < s->My) {
        cj = jy * s->Mx + jx; // indice cella
        for (k = s->cell_count[cj]; k < s->cell_count[cj+1]; k++) { //ecco perché mi serve cell_count[M+1]
          j = s->cell_ptr[k]; // indice particella
          printf("%d (%d: %d %d) ", j, cj, jx, jy);
        }
      }
    }
  }
  printf("\n");
}

printf("\n interactions:\n");
#endif

    

      
    // 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;
    }
    
    
    for (ix=0; ix < s->Mx; ix++) { 
      ci = iy * s->Mx + ix;
      for (iy=0; iy < s->My; iy++) {
        for (ii=0; ii<5; ii++) { // loop su "meta" celle adiacienti
          jx = ix + Delta[ii][0];
          jy = iy + Delta[ii][1];
          if (jx >=0 && jx < s->Mx && jy >=0 && jy < s->My) {
            cj = jy * s->Mx + jx; // indice cella
            for (ki = s->cell_count[ci]; ki < s->cell_count[ci+1]; ki++) {
              i = s->cell_ptr[ki]; // indice prima particella
              for (kj = s->cell_count[cj]; kj < s->cell_count[cj+1]; kj++) {
                j = s->cell_ptr[kj]; // indice seconda particella
                if (!(ii==0 && j <= i)) {     // loop su metà delle particelle della cella 0 
#ifdef DEBUG
printf("interaction %d - %d\n", i, j);
#endif
                  rx = s->p[j].x - s->p[i].x; 
                  ry = s->p[j].y - s->p[i].y; 
                  r2 = rx*rx+ry*ry;
                  if (r2 < s->rc2) {
                    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) - s->shift;
                  }
                }
              }
            }
          }
        }
      }
    }
    // gravity
    if (s->gravity > 0) {   
      for (i=0; i<s->N; i++) {
        s->p[i].fy -= s->gravity * s->p[i].m;
      }
    }

#ifdef DEBUG
getchar();
#endif
    //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 %f %f\n",s->N, s->LX, s->LY,s->gravity, s->deltat, s->cell_size, s->rc);
  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, cell_size, rc;
  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 %lf %lf\n", &N, &LX, &LY,&gravity, &deltat, &cell_size, &rc);
  s = new_simulation(N, LX, LY, gravity, deltat, cell_size, rc);
  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;
  }
}





void init_visualize(FILE *gp, simulation * s) {
  
  fprintf(gp, "set size square; unset key\n");
  fprintf(gp, "unset arrow; unset label\n");
  fprintf(gp, "set xrange [0:%f]; set yrange [0:%f]\n", s->LX, s->LY);
#ifdef VISDEBUG
  int i, j;
  for (i=0; i< s->Mx; i++) {
    fprintf(gp, "set arrow from %f,0 to %f,%f nohead lt 0 lc 'red'\n", i*s->cell_size,i*s->cell_size, s->LY);
  } 
  for (i=0; i< s->My; i++) {
    fprintf(gp, "set arrow from 0,%f to %f,%f nohead lt 0 lc 'red'\n", i*s->cell_size, s->LX,i*s->cell_size);
  } 
  for (j=0; j<s->My; j++) {
    for (i=0; i<s->Mx; i++) {
      fprintf(gp, "set label '(%d)' at %f,%f textcolor 'red'\n", j*s->Mx+i,
          (i+0.1)*s->cell_size, (j+0.2)*s->cell_size);
    }
  }
#endif
}

  
void visualize(FILE * gp, simulation * s, char * title) {
  int i;
   
  fprintf(gp, "set title '%s'\n", title);
#ifdef VISDEBUG
  fprintf(gp, "plot '-' w p pt 6 ps %f, '-' u 1:2:3 w labels offset 1.2\n", 50/s->LX);
  for (i=0; i<s->N; i++) {
    fprintf(gp, "%f %f\n", s->p[i].x, s->p[i].y);
  }
  fprintf(gp, "e\n");
  
  for (i=0; i<s->N; i++) {
    fprintf(gp, "%f %f %d\n", s->p[i].x, s->p[i].y, i);
  }
  fprintf(gp, "e\n");
#else

  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);
  }
#endif
  fflush (gp);
}
  

histo * new_histo(int nbin, double min, double max) {
  histo * h;
  
  h = calloc (1, sizeof(histo));
  h->nbin = nbin;
  h->max = max;
  h->min = min;
  h->c = calloc(nbin, sizeof(int));
  h->n=0;
  return h;
}

void destroy_histo (histo * h) {
  free(h->c);
  free(h);
}

void clear_histo(histo * h) {
  int i;
  
  h->n = 0;
  for (i=0; i < h->n; i++) {
    h->c[i] = 0;
  }
}


int add_histo(histo * h, double x) {
  int k;
  
  k = (int) (h->nbin * (x - h->min)/(h->max - h->min));
  if (k >= 0 && k < h->nbin) {
    h->c[k]++;
    h->n++;
    return(h->n);
  } else {
    return(0);
  }
}


vhisto * new_vhisto(int nbin, double min, double max, char * descr) {
  vhisto * v;
  
  v = calloc(1, sizeof(vhisto));
  v->gp = popen("gnuplot", "w");
  v->h = new_histo(nbin, min, max);
  v->descr = calloc(100, sizeof(char));
  strncpy(v->descr, descr, 100);
  return v;
}

void destroy_vhisto(vhisto * v) {
  fclose(v->gp);
  destroy_histo(v->h);
  free(v->descr);
  free(v);
}

void add_distr(vhisto *v, vhisto *y, simulation * s) {
  int i;
  double v2;
  
  for (i=0; i < s->N; i++) {
    v2 = pow(s->p[i].vx,2)+pow(s->p[i].vy,2);
    add_histo(v->h, v2);
    add_histo(y->h, s->p[i].y);
  }
}
  
void visualize_distr(vhisto * v) {
  int i;
  double x, y;

  fprintf(v->gp, "set title '%s'\n", v->descr);
  
  fprintf(v->gp, "plot '-' w l\n");
  
  for (i=0; i< v->h->nbin; i++) {
    x = (i+0.5)* (v->h->max - v->h->min)/v->h->nbin + v->h->min;
    y = (double) v->h->c[i]/v->h->n;
    fprintf(v->gp, "%f %f\n", x, y);
  }
  fprintf(v->gp, "e\n");
  fflush(v->gp);
}