/* the Wolf single cluster algorithm for the 2D Ising model */
/* a basic program, making just the dynamics and time-update
	for a fixed temperature */

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

void init();
void initconf();
void sw();
void flip();
void rec(int i, int j);
void pconf();

#define Dim  16        /* size of the square lattice */
#define MC_steps 200   /* number of the MC steps    */
#define ToTc 1.2       /* value of T / Tc           */


int S[Dim+2][Dim+2];  /* vector for the constructon of the clusters */
int W[Dim+2][Dim+2];  /* vector for the values of the spins */
double a, timp;
int p, NN, is, www;

/*-----------------------------------------------*/
/* generating integer random numbers 1,2,3 ....q */

int randint(int q)
{int d;
 d=(int)((double)(rand())/(RAND_MAX+1.0)*q)+1;
 return(d);
}


/*-----------------------------------------------*/
/* generating float random numbers in the [0,1) intervall */

double randfloat()
{
 double g;
 g=(double)(rand())/(RAND_MAX+1.0);
 return(g);
}


/* ----------------------------------------------*/

/* the main part of the program */

void main(void)
{
  NN=Dim*Dim;
  a=(1.0-exp(-2.0*.4409/ToTc));
  initconf();
  for (is=0; is<MC_steps; is++)
  {
	 init();
	 sw();
	 pconf();
	 getchar();
  }

}


/*-----------------------------------------*/
/* realizing the Wolf dynamics */

void sw()
{
  int lx,ly;
  lx=randint(Dim);
  ly=randint(Dim);
  S[lx][ly]=1;
  p=1;
  www=W[lx][ly];
  rec(lx,ly);
  timp=timp+(double)(p)/NN;
  return;
}
       

/*----------------------------------------------------*/
/* recursion subrutine for constracting the S-W clusters */

void rec( int i, int j)
{
p++;
	  if  (S[i+1][j]==0)
	  if ((www*W[i+1][j]==1) && (randfloat() < a ))
		{ S[i+1][j]=1;
		  W[i+1][j]=-W[i+1][j];
               	  rec(i+1,j);
		}

	if (S[i][j+1]==0)
	if ((www*W[i][j+1]==1) && ( randfloat() < a  ))
		 {S[i][j+1]=1;
		  W[i][j+1]=-W[i][j+1];
    	          rec(i,j+1);
                 }

        if (S[i-1][j]==0)
	if ((www*W[i-1][j]==1) &&  (randfloat() < a ))
		 { S[i-1][j]=1;
		   W[i-1][j]=-W[i-1][j];
	           rec(i-1,j);
                 }

   if (S[i][j-1]==0)
	if ((www*W[i][j-1]==1) &&  (randfloat() < a ))
		 {  S[i][j-1]=1;
		    W[i][j-1]=-W[i][j-1];
	            rec(i,j-1);
	         }
 return;
}


/* --------------------------------------------------------------------*/
/* initializing the vector for the clusters */

void init()
{
  int lx,ly;
  for (lx=1 ; lx<=Dim ;lx++)
	  for (ly=1 ;ly<=Dim ;ly++)
		 S[lx][ly]=0;
 return;
 }

/* -------------------------------------------- */
/* initial configuration of the spins + making the boundaries */

void initconf()
{int lx, ly;
for (lx=1 ; lx<=Dim ;lx++)
	  for (ly=1 ;ly<=Dim ;ly++)
		 W[lx][ly]=1;
for (lx=1; lx<=Dim+1 ; lx++)
		 {S[0][lx]=2;
		  S[lx][0]=2;
		  S[Dim+1][lx]=2;
		  S[lx][Dim+1]=2;
		 }
return;
}

/* -------------------------------------- */
/* printing the configuration on the screen */

void pconf()
{int ix, iy;
 printf("         %f\n", timp);
 for (iy=1; iy<=Dim; iy++)
	 {for(ix=1; ix<=Dim; ix++)
		 printf("%d ", (W[ix][iy]+1)/2);
		 printf("\n");
	 }
return;
}


/* --------------------------------------- */