#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#define L0 32
#define L1 32
#define V L0*L1
#define D 2
int spin[V];
int nlist[V][2*D];


double p[9]; 
double pdelete;
FILE *fp1;

void metro();
void sw_cluster();
double drand48(); 

main()
{
int i0,i1,i0p,i0m,i1p,i1m;
int i,d;
int iu,im,imax,order,times;
int mag,ene; 
double beta; 
unsigned short int seed16v[3];
printf("give a seed for the random number generator \n"); 
scanf("%d",&seed16v[0]);
seed48(seed16v);

fp1=fopen("output","w");

printf("give beta, the number of measurements and the number of \n");
printf("sweeps per measurement \n");
scanf("%lf",&beta); 
scanf("%d %d",&imax,&times); 
printf("random start =0 ; ordered start=1 \n");
scanf("%d",&order); 

fprintf(fp1,"%d %d %d %f \n",L0,L1,times,beta); 

/* compute the acceptance probabilities */
p[0]=1.;
p[2]=1.;
p[4]=1.;
p[6]=exp(-4.*beta);
p[8]=exp(-8.*beta);

/* compute exp(-2*beta) for the cluster */
pdelete=exp(-2.*beta);

/* define the neighbours of a lattice point */
for(i0=0;i0 < L0;i0++)
  {
  /*take periodic boundary conditions into account*/
  i0p=i0+1; if(i0p==L0) i0p=0;  i0m=i0-1; if(i0==0) i0m=L0-1;

  for(i1=0;i1 < L1;i1++)
    {
    /*take periodic boundary conditions into account*/
    i1p=i1+1; if(i1p==L1) i1p=0; i1m=i1-1; if(i1==0)   i1m=L1-1;

    nlist[i0*L1+i1][0]=i0p*L1+i1;
    nlist[i0*L1+i1][1]=i0*L1+i1p;
    nlist[i0*L1+i1][2]=i0m*L1+i1;
    nlist[i0*L1+i1][3]=i0*L1+i1m;
    }
  }

/*initialize the field variables */
if(order==1)
  {
  for(i=0; i < V; i++)  spin[i]=1;
  }
else
  {
  for(i=0;i < V;i++) if(drand48()>0.5) spin[i]=1; else spin[i]=-1;
  }

for(iu=1;iu<=imax;iu++)
  {
  for(im=1;im<=times;im++) sw_cluster(); 
/*  for(im=1;im<=times;im++) metro(); */

/*measure the magnetisation and the energy*/
  ene=0;
  mag=0; 
  for(i=0;i<V;i++)
    {
    mag+=spin[i];
    for(d=0; d<D; d++) ene+=spin[i]*spin[nlist[i][d]];
    }
    fprintf(fp1,"%d %d \n",ene,mag); 
  }
}

void metro() 
{
int i,d,neigh; 

/*sweep over the lattice*/
for(i=0;i < V;i++)
  {
  neigh=0;
  /*take the sum over the nearest neighbours*/
  for(d=0; d<2*D; d++) neigh+=spin[nlist[i][d]];
  neigh*=spin[i];
  /* decide whether to accept the proposal spin to -spin */
  if(p[4+neigh]>drand48()) 
    {
    spin[i]=-spin[i];
    }
  }
}

void sw_cluster() 
{
/* for the d's of the cluster algorithm  */
int d_array[V][D];
/* for the Hoshen-Kopelman */
int label_array[V];
int cluster_sign[V];
int i,i2,im,in,in1,in2,in_alt;
/* compute the d's for fixed spins : */
for(i=0;i < V;i++)
  {
  for(im=0;im<D;im++)
    {
    if(spin[i]==spin[nlist[i][im]])
      {
      if(drand48()>pdelete) d_array[i][im]=1; else d_array[i][im]=0;
      }
    else
      {
      d_array[i][im]=0;
      }
    }
  label_array[i]=i;  /* initialize the label array for the cluster search*/
  }
/*  determine the clusters using the Hoshen-Kopelman algorithm */

for(i=0;i < V;i++)
  {
  for(im=0;im<D;im++)
    {
    if(d_array[i][im])
      {
      i2=nlist[i][im];
      in1=i;  do {in_alt=in1; in1=label_array[in1];}  while(in1<in_alt);
      in2=i2; do {in_alt=in2; in2=label_array[in2];}  while(in2<in_alt);
      /* set the labels to the minimum of the labels on the edge */
      if(in1>in2) 
        {label_array[in1]=in2; label_array[i]=in2; label_array[i2]=in2;}
      else if(in1<in2) 
	{label_array[in2]=in1; label_array[i]=in1; label_array[i2]=in1;}
      }
    }
  }
/* go through to the fixed point of the mapping im=label_array[im]; */
for(i=0;i<V;i++) {in=i;
do {in_alt=in; in=label_array[in];}  while(in<in_alt); label_array[i]=in;}

/*  give a new sign to the clusters */
for(i=0;i < V;i++)
  {
  cluster_sign[i]=1; if(drand48()>0.5) cluster_sign[i]=-1;
  }
for(i=0;i < V;i++)
  {
  spin[i]=cluster_sign[label_array[i]];
  }
}