/* a basic Glauber MC algorithm for the 2D Ising model */

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


#define Npoints 1000
#define Kiir 10
#define  DimX  14
#define  DimY  14
#define T    1.10

int M[DimX+2][DimY+2];
double co;

/* ==================================================== */
							/*generates random numbers between 0 and 1 */
double randf()
{return( (double)(rand())/(double)(RAND_MAX+1.0));
}



/* ======================================================*/

int ram(int wi)			 /* gives random numbers from 1.....wi */
{
 int wo;
 wo=(int)((double)(rand())/(double)(RAND_MAX+1.0)*wi)+1;
 return(wo);
}


/* =====================================================*/

void pconf()              /* prints out the spin configuration */
{int ix, iy;
 for (iy=1; iy<DimY; iy++)
	 {for(ix=1; ix<DimX; ix++)
		 printf("%d ", (M[ix][iy]+1)/2);
		 printf("\n");
	 }
return;
}

/* ==================================================== */

void boundary()                /* makes the free boundary conditions */
{ int ik;
  for(ik=1; ik<=DimX; ik++)
		{M[ik][0]=0;
		 M[ik][DimY+1]=0;
		}
  for(ik=1; ik<=DimY; ik++)
		{M[0][ik]=0;
		 M[DimX+1][ik]=0;
		}
return;
}


/* =================================================  */


void init()                /* initializes the spin configuration */
{ int ii, jj, pt;
  for (ii=1; ii<DimY+1; ii++)
	 {
		 for (jj=1; jj<DimX+1; jj++)
			{pt=2*ram(2)-3;
			 M[jj][ii]=pt;
			}
	 }
boundary();
return;
}

/*================================================== */

void flip(int fx, int fy)          /* flip the spins and drive the systems
											  towards the desired canonical distribution */
{int sum, spin, pr;
 double p;
 float r, DE;
 spin=M[fx][fy];
 sum=M[fx+1][fy]+M[fx-1][fy]+M[fx][fy+1]+M[fx][fy-1];
 DE=2*spin*sum;
		  p=exp(-co*DE);
		  p=p/(1.0+p);
		  r=randf();
		  if (r<=p)  M[fx][fy]=-spin;


return;
}

/* ===========================================================*/

main()
{
 int i, j , x , y, NN;
 int t;
  NN=DimX*DimY;
  co=1.0/T;
  init();
  for(i=1; i<Npoints+1; i++)
	 {
		 printf("%d\n",i);
		 for(j=1; j<Kiir+1; j++)
		 {
			  for(t=1; t<NN+1; t++)
			  {
			  x=ram(DimX);
			  y=ram(DimY);
			  flip(x,y);
			  }
		  }
      pconf();
     getchar();

	  }

}


/* =========================================================== */