#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;
  int cell;
} particle;

typedef struct {
  int N;
  particle *p;
  double LX, LY;
  double delta;
  double temp;
  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 temp, double gravity,
    double delta, double cell_size, double rc);
void evolve(simulation *s);
void random_init(simulation *s); 
void visualize(FILE * gp, simulation * s, char * title);
double compute_energy(simulation *s);
double compute_energy_particle(simulation *s, int k);
int save(simulation *s, char * filename);
simulation * load(char * filename);
void destroy_simulation(simulation * s);
void init_visualize(FILE *gp, simulation * s);


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];
  char title[100];
  clock_t start, end;
  double cpu_time_used;
  double gravity;
  double cell_size;
  double rc;
 
  
  T = 100000;
  VIS = 1000;
  temp = 1; // temperature
  delta = 0.1;
  LX=20;
  LY=20;
  cell_size = 2.5; rc = 2.5;
  N=100;
  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 delta=%f \n", N, LX, LY, delta);
    printf(" gravity=%g temp=%f T=%d VIS=%d\n", 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");
#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, &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, cell_size, rc);
          init_visualize(gp, s);
          random_init(s);    
          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->delta); 
          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++) {
            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 relative temperature=%f\n", s->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;
    }
  }
}

// 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) {
  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;
    // inserisco le particelle nelle celle
    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 temp, double gravity, 
    double delta, double cell_size, double rc) {
  simulation *s;
  
  s=calloc(1, sizeof(simulation));
  s->N = N;
  s->LX = LX;
  s->LY = LY;
  s->delta = delta;
  s->temp=temp;
  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 += s->gravity * s->p[i].m * s->p[i].y;
  }
  return(energy);
}

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


double compute_energy_particle(simulation *s, int i) {
  int  j;
  double energy;
  double rx, ry, r2, rm2, rm6;
  
  //calcolo dell'energia
  int ii, ix, iy, jx, jy, ci, cj;
  node *nj;  

  energy = 0;
  
  ci = s->p[i].cell;
  ix = ci % s->Mx;
  iy = ci / s->My;

  for (ii=0; ii<9; ii++) { // loop sulle 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 (nj = s->cell[cj]; nj != NULL; nj=nj->next) {
        j = nj->v; // indice seconda particella
        if (!(ii==0 && j == i)) {     // loop sulle altre  particelle 
          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;
            energy += 4*rm6*(rm6-1) - s->shift;
          }
        }
      }
    } 
  }
  energy += s->p[i].y * s->p[i].m * s->gravity;
  return(energy);
}



void evolve(simulation *s) {
  int i, j;
  double energy, oldenergy, oldx, oldy;
  double deltaE;
  
  for (i=0; i<s->N; i++) { // passo mont-carlo
    j = (int) (drand48() * s->N); //prendo una particella a caso
    oldenergy = compute_energy_particle(s, j);
    oldx = s->p[j].x; // salvo la posizione
    oldy = s->p[j].y;
    s->p[j].x += s->delta * (2 * drand48()-1); // sposto
    s->p[j].y += s->delta * (2 * drand48()-1);
    // condizioni al bordo
    if (s->p[j].x < 0) {
      s->p[j].x *= -1;
    }
    if (s->p[j].y < 0) {
      s->p[j].y *= -1;
    } 
    if (s->p[j].x > s->LX) {
      s->p[j].x = 2 * s->LX - s->p[j].x;
    }
    if (s->p[j].y > s->LX) {
      s->p[j].y = 2 * s->LY - s->p[j].y;
    } 
    energy = compute_energy_particle(s, j);
    deltaE = energy - oldenergy;
    if (deltaE <= 0 || drand48() < exp(-deltaE / s->temp)) { 
      //accetto 
      // sposto la particella nella cella (se cambiata)
      int ix, iy, k;
      ix = (int) (s->p[j].x / s->cell_size);
      iy = (int) (s->p[j].y / s->cell_size);
      k = iy * s->Mx + ix; // numero cella
      if (k != s->p[j].cell) { // cambia cella
        find_remove(s->cell + s->p[j].cell, j); // si rimuove dalla ll
        push(s->cell + k, j); // si inserisce nella nuova ll
        s->p[j].cell=k; // memorizzo la cella nella struttura.  
      }
    } else {
      //revert
      s->p[j].x = oldx;
      s->p[j].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 %f %f\n",s->N, s->LX, s->LY,s->temp, s->gravity, s->delta, s->cell_size, s->rc);
  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, cell_size, rc;
  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 %lf %lf\n", &N, &LX, &LY, &temp, &gravity, &delta, &cell_size, &rc);
  s = new_simulation(N, LX, LY, temp, gravity, delta, cell_size, rc);
  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;
     // 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 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);
}