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

/* function prototypes */
void IsingStep(int N, int **State, double JoverT, int *pM);
void PrintStep(int N, int **State);

int main(int argc, char **argv)
{
  int i,j;              /* loop variables                 */
  int N;                /* system size                    */
  int Nrounds;          /* number of iterations           */
  int M;                /* magnetization                  */
  int **State;          /* state of all the spins         */
  double JoverT;        /* J/T                            */
  FILE *gfp;            /* file pointer to gnuplot        */

  /* interpret command line arguments */
  if (argc!=4) {
    fprintf(stderr,"ising usage: ising <N> <Nrounds> <J/T>\n");
    return(1);
  }
  N=atoi(argv[1]);
  Nrounds=atoi(argv[2]);
  JoverT=atof(argv[3]);

  /* get memory for state array */
  State=(int **)malloc((N+2)*sizeof(int *));
  if (State) {
    State[0]=(int *)malloc((N+2)*(N+2)*sizeof(int));
    if (!State[0]) {
      free((void *)State);
      State=NULL;
    }
  }
  if (!State) {
    fprintf(stderr,"Not enough memory !\n");
    return(2);
  }
  for(i=1;i<N+2;i++)
    State[i]=State[0]+i*(N+2);

  /* initialize state */
  /* boundary */
  for(i=0;i<N+2;i++)
    State[0][i]=State[i][0]=State[N+1][i]=State[i][N+1]=0;
  /* and interior */
  M=0;
  for(i=1;i<=N;i++)
    for(j=1;j<=N;j++) {
      if (drand48()<0.5)
        State[i][j]=1;
      else
        State[i][j]=-1;
      M+=State[i][j];
    }

  /* get gnuplot file pointer */
  gfp=popen("gnuplot","w");
  if (!gfp) {
    fprintf(stderr,"Could not open gnuplot !\n");
    free((void *)State[0]);
    free((void *)State);
    return(3);
  }

  /* do simulation */
  for(i=0;i<Nrounds;i++) {
    IsingStep(N,State,JoverT,&M);
    if (!(i%1000)) {
      PrintStep(N,State);
      fprintf(gfp,"set title \"M=%d %d Steps\"\n",M,i);
      fprintf(gfp,"plot [0:%d] [0:%d] \"isingup.dat\" t \"\"\n",N+1,N+1);
      fflush(gfp);
    }
  }

  /* free memory */
  fclose(gfp);
  free((void *)State[0]);
  free((void *)State);

  return(0);
}

void IsingStep(int N, int **State, double JoverT, int *pM)
{
  int x,y;    /* coordinates of random position     */
  int S;      /* sum over spins around current spin */
  double dE;  /* change in energy                   */

  /* choose a random position between 1 and N */
  x=(int)(N*drand48())+1;
  y=(int)(N*drand48())+1;
  /* sum over spins around current spin (note that there is an edge of
     zero spins around the lattice to avoid having to do something special
     at the edge */
  S=State[x+1][y]+State[x][y+1]+State[x-1][y]+State[x][y-1];
  /* calculate the change in energy if we flip the spin at (x,y) */
  dE=JoverT*S*State[x][y];
  /* if it is advantageous, we do it right away, otherwise we flip a coin */
  if (dE<0.0 || drand48()<exp(-dE)) {
    *pM-=2*State[x][y];
    State[x][y]=-State[x][y];
  }
}

void PrintStep(int N, int **State)
{
  int i,j;   /* loop variables          */
  FILE *fp;  /* file pointer for output */

  fp=fopen("isingup.dat","w");
  if (fp) {
    for(i=1;i<=N;i++)
      for(j=1;j<=N;j++)
        if (State[i][j]>0)
          fprintf(fp,"%d %d\n",i,j);
    fclose(fp);
  }
  fp=fopen("isingdown.dat","w");
  if (fp) {
    for(i=1;i<=N;i++)
      for(j=1;j<=N;j++)
        if (State[i][j]<0)
          fprintf(fp,"%d %d\n",i,j);
    fclose(fp);
  }

  for(j=0;j<100;j++)
    for(i=0;i<100000;i++);
}