#include "ClusterUpdate.h"

void ClusterUpdate(ScalarField &phi){
LatticeVector x,active_sites[SITES],new_active_sites[SITES],e_mu;
int mu,ACTIVE,newcount=0,count=0,cluster_count;
ScalarField cluster=ScalarField();

//cout << "In cluster update\n" << flush;

// random site
for(mu=0;mu<D;mu++){
x.set(mu,(int)(RandomNumber()*L));
}

// flip sign of field on site
phi.set(x,-phi.get(x));


// mark site in cluster
cluster.set(x,1.0);
cluster_count=1;

// add x to array of active sites
active_sites[0]=x;

ACTIVE=1;

// iterate to grow cluster
do{

// loop over active sites
for(count=0;count<ACTIVE;count++){
for(mu=0;mu<D;mu++){

e_mu=LatticeVector(mu);
x=active_sites[count];
//cout << "active site\n";
//x.print();
//cout << flush;

if(cluster.get(x+e_mu)==0.0){
   if(phi.get(x)*phi.get(x+e_mu)<0.0){
      if(RandomNumber()<(1.0-exp(2*kappa*phi.get(x)*phi.get(x+e_mu)))){
//      cout << "linked neighbor\n";
//      (x+mu).print();
//      cout << flush;
      cluster.set(x+e_mu,1.0);
      cluster_count++;
      new_active_sites[newcount]=x+e_mu;
      newcount++;
      phi.set(x+e_mu,-phi.get(x+e_mu));
      }
   }
}

if(cluster.get(x-e_mu)==0.0){
   if(phi.get(x)*phi.get(x-e_mu)<0.0){
      if(RandomNumber()<(1.0-exp(2*kappa*phi.get(x)*phi.get(x-e_mu)))){
//      cout << "linked neighbor\n"; 
//      (x-mu).print();
//      cout << flush;
      cluster.set(x-e_mu,1.0);
      cluster_count++;
      new_active_sites[newcount]=x-e_mu;
      newcount++;
      phi.set(x-e_mu,-phi.get(x-e_mu));
      }
   }
}


}
}

// copy back and iterate from new active list


ACTIVE=newcount;
newcount=0;
for(count=0;count<ACTIVE;count++){
active_sites[count]=new_active_sites[count];
}

}while(ACTIVE>0);

//cout << "size of cluster " << cluster_count << "\n" << flush;
return;
}