/* calculating <m(T)> , <Cv(T)> and <Chi(T)> for the 2D Ising model
    by using the Metropolis algoritm */

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


#define Nmax 1000
#define Ntrans 1000
#define  DimX  14
#define  DimY  14

FILE *fp;

int M[DimX+2][DimY+2];
long int Mag, Ener;
double co, T, Emed, Mmed, Emed2, Mmed2, MMM;

/* ==================================================== */
							/*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;
  Mag=0;

  for (ii=1; ii<DimY+1; ii++)
	 {
		 for (jj=1; jj<DimX+1; jj++)
			{pt=1;
			 M[jj][ii]=pt;
			 Mag+=pt;
			}
	 }
boundary();

	Ener=0;
	for (ii=1; ii<DimY+1; ii++)
	 {
		 for (jj=1; jj<DimX+1; jj++)
			{
			  Ener+=-M[ii][jj]*(M[ii+1][jj]+M[ii-1][jj]+M[ii][jj+1]+M[ii][jj-1]);
			}
	 }
	 Ener=Ener/2;
return;
}

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

void flip(int fx, int fy)          /* flip the spins and drive the systems
											  towards the desired canonical distribution */
{int sum, spin, pr, DE;
 double p;
 double r;
 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;
		  if (DE>0) p=exp(-co*DE);
		  else p=1;
		  r=randf();
		  if (r<=p)  {M[fx][fy]=-spin;
						  Mag+=-2*spin;
						  Ener+=DE;
						  }


return;
}

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

main()
{
 int i, j , x , y, NN;
 int t;
 double Cv, chi;
  NN=DimX*DimY;
  for(T=0.5; T<=4.0; T+=0.1)
  {
  Emed=0; Emed2=0; Mmed=0; Mmed2=0; MMM=0;
  printf("%lf\n",T);
  co=1.0/T;
  init();
  for(i=1; i<=Ntrans; i++)
	 {
			  for(t=1; t<NN+1; t++)
			  {
			  x=ram(DimX);
			  y=ram(DimY);
			  flip(x,y);
			  }

	  }

  for(i=1; i<=Nmax; i++)
	 {
			  for(t=1; t<NN+1; t++)
			  {
			  x=ram(DimX);
			  y=ram(DimY);
			  flip(x,y);
			  }

			  Emed+=Ener;
			  Mmed+=Mag;
			  if (Mag>0) MMM+=Mag;
			  else MMM+=-Mag;
			  Emed2+=Ener*Ener;
			  Mmed2+=Mag*Mag;


	  }

	  Mmed=Mmed/(double)(Nmax);
	  MMM=MMM/(double)(Nmax)/(double)(NN);
	  Mmed2=Mmed2/(double)(Nmax);
	  Emed=Emed/(double)(Nmax);
	  Emed2=Emed2/(double)(Nmax);
	  Cv=(Emed2-Emed*Emed)/T/T/(double)(NN);
	  chi=(Mmed2-Mmed*Mmed)/T/(double)(NN);

	  fp=fopen("magn.dat","a");
	  fprintf(fp,"%lf      %lf        %lf        %lf\n", T,  MMM, Cv, chi);
	  fclose(fp);

  }

}


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