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

#define sigma(i,j) _sigma[((j+L) % L)  * L + ((i+L) % L)]

void visualize(FILE * gp, int * conf, int L, char * title);
double magnetization(int * x, int L);

typedef struct {
  int * data;
  int nbin;
  double hmin, hmax;
  int count;
  double * distr;
  double * x;
} histogram;

histogram * new_histogram (int nbin, double hmin, double hmax);
int add_histogram(histogram * h, double x);
void histogram2distr (histogram * h);
int zero_histogram (histogram * h);

int main(){
  int L;
  int *_sigma;
  int i, j;
  double deltaH, Temp;
  int t, T, TRANS, l;
  double m, summ;
  double J;
  FILE *gp;
  FILE *gp1;
  FILE *gp2;
  char prompt[100];
  histogram * h;

  L = 50; 
  T = 30;
  TRANS= 10;
  _sigma = calloc(L*L, sizeof(int));
  J = 1;
  Temp = .01;
  
  h = new_histogram(100, -1., 1.);

  gp = popen("gnuplot", "w");
  fprintf(gp, "set size square; unset key; set cbrange [-1:1], set xrange [0:%d]; set yrange [0:%d]\n", L-1, L-1);
  gp1 = popen("gnuplot", "w");
  fprintf(gp1, "unset key; set xlabel 'Temp'; set ylabel 'm' norotate\n");
  gp2 = popen("gnuplot", "w");
  fprintf(gp2, "unset key; set xlabel 'm'; set ylabel 'P' norotate\n");

  srand48(time(NULL));

  // inizializziamo a caso
  for (i=0; i<L; i++) {
    for (j=0; j<L; j++) {
      sigma(i,j) = drand48() < 0.5 ? -1 : 1;
    }
  }
  fprintf(gp1, "plot '-' w l\n");
  for (Temp = 3; Temp >=0.01; Temp -= 0.1) {
    summ = 0;
    zero_histogram (h);
    for (t=0; t<TRANS+T; t++) { // passo MC
      for (l=0; l<L*L; l++) { //tentativi di cambio
        i = (int) (drand48() * L);
        j = (int) (drand48() * L);
        deltaH = sigma(i,j)* 2* J * (sigma(i-1,j)+sigma(i, j-1)+sigma(i+1,j) + sigma(i, j+1));
        if (deltaH < 0 || drand48() < exp(-deltaH/Temp)) { // metropolis
          sigma(i,j) *= -1;
        }
      }
      if (t>TRANS) {
        m = magnetization(_sigma, L);
        summ += m;
        add_histogram(h, m);
      }
    }
    fprintf(gp1, "%f %f\n", Temp, summ/T);
    sprintf(prompt, "Temp=%f", Temp);
    visualize(gp, _sigma, L, prompt);
    histogram2distr(h);
    fprintf(gp2, "plot '-' w l\n");
    for (i=0; i< h->nbin; i++) {
      fprintf(gp2, "%f %f\n", h->x[i], h->distr[i]);
    }
    fprintf(gp2, "end\n");
    fflush(gp2);
  }
  fprintf(gp1, "end\n");
  fflush(gp1);
  getchar();
}

double magnetization(int * x, int L) {
  int i, j;
  double m;
  
  m = 0;
  for (i=0; i<L; i++) {
    for (j=0; j<L; j++) {
      m += x[j*L+i];
    }
  }
  return(m/(L*L));
}

void visualize(FILE * gp, int * c, int L, char * prompt) {
  fprintf(gp, "set title '%s'\n", prompt);
  fprintf(gp, "plot '-' binary array=%dx%d w image \n", L,L);
  fwrite(c, sizeof(int),L*L, gp);
  fflush(gp);
}

histogram * new_histogram (int nbin, double hmin, double hmax) {
  histogram * h;
  
  h = calloc(1, sizeof(histogram));
  h->nbin=nbin;
  h->hmin=hmin;
  h->hmax=hmax;
  h->data = calloc(nbin, sizeof(int));
  h->distr = calloc(nbin, sizeof(double));
  h->x = calloc(nbin, sizeof(double));
  h->count = 0;
  
  return(h);
}
  
  
int add_histogram(histogram * h, double x) {
  int i;
  
  if (x >= h->hmin && x < h->hmax) {
    i = (int) ((x - h->hmin)/ (h->hmax - h->hmin) * h->nbin);
    h->data[i] ++;
    h->count ++ ;
    return(h->count);
  } else {
    return(-1);
  }
}

void histogram2distr (histogram * h) {
  int i;
  
  for (i=0; i< h->nbin; i++) {
    h->distr[i] = (double) h->data [i] * (h->hmax - h->hmin) / (h->count * h->nbin);
    h->x[i] = h->hmin + (i+0.5)*(h->hmax - h->hmin)/h->nbin;
  }
}
  
int zero_histogram (histogram * h) {
  int i;
  
  if (h) {
    h->count = 0;
    for (i=0; i < h->nbin; i++) {
      h->data[i]=0;
      h->x[i] = 0;
    }
    return(0);
  } else {
    return(1);
  }
}