/* The Swendsen and Wang MC for the 2-D Ising model */
/* A basic program, calculating the fluctuations in magnetization
	as a function of the temperature of the system */


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

void init();
void inW();
void sw();
void flip();
void anal();
void rec(int i, int j);


#define Dim  20
#define MC_steps 800
#define Start_nr 400

FILE *fp;

int S[Dim+2][Dim+2];  /* vector for labeling the clusters */
int W[Dim+2][Dim+2];  /* values of the spins */
int Q[Dim*Dim];  /* new values for the spins in a given cluster */
double a;
int p, NN, is;
double Mmed, Mmed2, ToTc ;


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

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

/* ----------------------------------------------*/
/* the main part of the program */

void main(void)
{
  fp=fopen("sw.dat","w");
  NN=Dim*Dim;
  for(ToTc=0.7; ToTc<1.4; ToTc=ToTc+0.05)
  {
  printf("%f\n", ToTc);
  fprintf(fp,"%f  ",ToTc);
  Mmed=0;
  Mmed2=0;
  a=(1-exp(-2*.4409/ToTc));
  init();
  Q[0]=1;
  inW();
  for (is=0; is<MC_steps; is++)
  {
	 init();
	 sw();
	 inW();
	 if (is>Start_nr-1) anal();
  }
  fprintf(fp," %f\n",Mmed2-Mmed*Mmed);
  }
  fclose(fp);
}

/* ---------------------------------------*/
/* redoing the values of the spins after a S-W MC step */

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

/*-----------------------------------------*/
/* the S-W dynamics */


void sw()
{

  int lx,ly;
  for(lx=1; lx<=Dim; lx++)
	 for(ly=1; ly<=Dim; ly++)
		{ if (S[lx][ly]==0)
			{ p=p+1 ;
			  S[lx][ly]=p;
			  rec(lx,ly);
			}
		}
  flip();
  return;
}


/*----------------------------------------------------*/
/* recursion subrutine for the construction of the S-W clusters */

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

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

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

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


/* -------------------------------------------------------- */
/* flipping the clusters */

void flip()
{ int lx;
  for (lx=1 ; lx<=p; lx++)
  { if ( randfloat() < 0.5 ) Q[lx]=1;
	 else Q[lx]=-1;
  }
  return;
}

/* ---------------------------------------- */
/* analyzing the spin configuration */


void anal()
{ int lx, ly;
  double m;
  m=0.0;
  for (lx=1; lx<=Dim ; lx++)
	 for (ly=1; ly<=Dim; ly++)
	m=m+(double)(W[lx][ly]);
	if (m<0) m=-m;
		m=m/(double)(NN);
  Mmed=Mmed*(is-Start_nr)/(is-Start_nr+1)+m/(is-Start_nr+1);
  Mmed2=Mmed2*(is-Start_nr)/(is-Start_nr+1)+m*m/(is-Start_nr+1);
  return;
}

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

/* initial configuratio of the cluster vector*/



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

}


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