#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 fVISDEBUG

// linked lists

typedef struct _node {
  int v;
  struct _node * next;
} node;


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;
  node ** cell;
  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);


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;
 
  
  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;
  
  // 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("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);
          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);
          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);
            }
          }
          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;
    }
  }
}

// linked lists
void push(node ** pstart, int v);
void print_ll(node * pstart);
int pop(node ** pstart);
node * remove_ll(node ** pstart);
void empty_ll(node ** pstart);
void find_remove(node ** pstart, int v);


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;
    // inserisco le particelle nelle celle (se cambiate)
    int ix, iy, k;
    ix = (int) (s->p[i].x / s->cell_size);
    iy = (int) (s->p[i].y / s->cell_size);
    k = iy * s->Mx + ix; // numero cella
    push(s->cell + k, i); // si inserisce nella nuova ll
    s->p[i].cell=k; // memorizzo la cella nella struttura.  
  }
#ifdef DEBUG
printf("\ncells:\n");
int ci, ix, iy;
node * ni;
for (ci=0; ci< s->Mx*s->My; ci++) {
  ix = ci % s->Mx; // recupero gli indici x, y della cella
  iy = ci / s->Mx;
  printf("cell[%d] (%d %d): ", ci, ix, iy); 
  for (ni = s->cell[ci]; ni != NULL;  ni=ni->next) { 
    printf("part[%d] ", ni->v);
  }
  printf("\n");
}
#endif

}


node * new_node(int v) {
  node * x;
  x = calloc(1, sizeof(node));
  x->v = v;
  x->next = NULL;
  return (x);
}

void push(node ** pstart, int v) {
  node * x;
  x = new_node(v);
  x->next = * pstart;
  * pstart = x;
}

void empty_ll(node ** pstart) {
  while (remove_ll(pstart));
}

node * remove_ll(node ** pstart) {
  node * temp;
  if (*pstart == NULL) {
    return(NULL);
  }
  temp = *pstart;
  *pstart = (*pstart)->next;
  free(temp);
  return(*pstart);
}

int pop(node ** pstart) {
  node * temp;
  if (*pstart == NULL) {
    printf("error! empty stack!\n");
    return(-1);
  }
  int v = (*pstart)->v;
  temp = *pstart;
  *pstart = (*pstart)->next;
  free(temp);
  return(v);
}

void find_remove(node ** pstart, int v) {
  node * temp;
  // find
  if (*pstart == NULL) {
    printf("error! empty stack!\n");
  }
  while ((*pstart)->v != v) {
    pstart = &((*pstart)->next);
    if (*pstart == NULL) {
      printf("error! not found!\n");
    }
  }
  // remove
  temp = *pstart;
  *pstart = (*pstart)->next;
  free(temp);
}
    

void print_ll(node * start) {
  while (start != NULL) {
    printf("%d\n", start->v);
    start = start->next;
  }
}


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=calloc(s->Mx * s->My, sizeof(node *)); // vedi evolve
  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);
  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;
      } 
    }
    // inserisco le particelle nelle celle (se cambiate)
    for (i=0; i< s->N; i++) {
      int ix, iy, k;
      ix = (int) (s->p[i].x / s->cell_size);
      iy = (int) (s->p[i].y / s->cell_size);
      k = iy * s->Mx + ix; // numero cella
      if (k != s->p[i].cell) {
        find_remove(s->cell + s->p[i].cell, i); // si rimuove dalla ll
        push(s->cell + k, i); // si inserisce nella nuova ll
        s->p[i].cell=k; // memorizzo la cella nella struttura.  
      }
    }



    // 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
    // tanto vale procedere per celle
    int ii, ix, iy, jx, jy, ci, cj;
    node * ni, *nj;  

#ifdef DEBUG
printf("\ncells:\n");
for (ci=0; ci< s->Mx*s->My; ci++) {
  ix = ci % s->Mx; // recupero gli indici x, y della cella
  iy = ci / s->Mx;
  printf("cell[%d] (%d %d): ", ci, ix, iy); 
  for (ni = s->cell[ci]; ni != NULL;  ni=ni->next) { 
    printf("part[%d] ", ni->v);
  }
  printf("\n");
}


 
printf("\nneighbors:\n");
for (i=0; i<s->N; i++) {
  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 (nj = s->cell[cj]; nj != NULL;  nj=nj->next) { 
          j = nj->v; // indice particella
          printf("%d (%d: %d %d) ", j, cj, jx, jy);
        }
      }
    }
  }
  printf("\n");
}
getchar();
printf("\nInteractions:\n");
#endif


    
    for (ix=0; ix < s->Mx; ix++) { 
      for (iy=0; iy < s->My; iy++) {
        ci = iy * s->Mx + ix;
        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
            // questa è la parte che cambia 
            for (ni = s->cell[ci]; ni != NULL; ni=ni->next) {
              i = ni->v; // indice prima particella
              for (nj = s->cell[cj]; nj != NULL; nj=nj->next) {
                j = nj->v; // indice seconda particella
            // fine cambiamento
                if (!(ii==0 && j <= i)) {     // loop su metà delle particelle della cella 0 
#ifdef DEBUG
printf("interaction (ii=%d ix=%d iy=%d - %d [%d %d]) %d - %d\n", ii,ix,iy, jy*s->Mx+jx, jx,jy, i, j);
getchar();
#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;
    // inserisco le particelle nelle celle (se cambiate)
    int ix, iy, k;
    ix = (int) (s->p[i].x / s->cell_size);
    iy = (int) (s->p[i].y / s->cell_size);
    k = iy * s->Mx + ix; // numero cella
    push(s->cell + k, i); // si inserisce nella nuova ll
    s->p[i].cell=k; // memorizzo la cella nella struttura.  
  }
  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);
}