function [nblock,block]=coarsening(bond)
% COARSENING Find clusters according to frozen bonds on a lattice
%      [nblock,block]=coarsening(bond) divides the points into
%      clusters according to frozen bonds in a 2D sqaure lattice
%      with nearest neighbor connectivity
% 
%      bond - N x 4 matrix:  
%           index of spins connected by frozen bonds to spin i
%           negative index means no bond
%
%      Example: 
%            2 10 11 -1  - spin 1 is connected to spins 2,10,11
%            1 3 -1  -1  - spin 2 is connected to spins 1,3
%      note that Bonds matrix can state each bond either once or twice
%
%      Results:
%      nblock: the number of clusters 
%      block: the index of the cluster to which spin i belongs
%


nblock=0;
N=size(bond,1);
stack=zeros(N,1);
block=zeros(N,1);

for i=1:N
  if block(i) == 0      % point i was not labeled yet
    ns=1;               % set stack index to 1
    nblock=nblock+1;
    block(i)=nblock;
    stack(ns)=i;
    while ns >= 1
      i0=stack(ns);
      ns = ns -1;
      for k=1:4
	if ( bond(i0,k) > 0 ) &  (bond(i0,k) <= N) & ...
	   ( block(bond(i0,k)) == 0 )
	   block( bond(i0,k) ) = nblock;
	   ns = ns +1;
	   stack(ns)=bond(i0,k);
	end
      end
    end
  end
end