% sand_piles
% C Bishop 
% We have a grid on a n by n torus and place a
% random number of grains of sand on each loation. 
% if there are four or more grains at one place, 
% four are moved, one each, to each of the neighboring 
% sites: up, down, left right. 
close all;
num=20;top=6;
delay=0;skip=1;limit=1000;
old=randi(3,num); % uniform 1,2,3
%old=randi(4,num); % uniform 1,2,3,4
p=.1;
%old= 3 *ones(num) + (rand(num)<p); % 3 + 1 with prob p
old= 2 *ones(num) + 2*(rand(num)<p); % 2 + 2 with prob p
%q=.1;
%old= 3*ones(num) + (rand(num)<p) - (rand(num)<q) ;% 
figure;axis equal; axis off; hold on;
image(old);
colormap(prism(top));
colorbar;
title(['\fontsize{30}step = 0']);
pause(5);
f=@(n) mod(n-1,num)+1;
new=zeros(num);
counter=0;
limit=100;
while counter < limit
    counter=counter+1;
    new=zeros(num);
   for k=1:num
      for j=1:num
          if old(k,j)<4
            new(k,j)=old(k,j)+new(k,j);
          else
            new(k,j)=old(k,j)-4+new(k,j);
            new(f(k-1),j)=new(f(k-1),j)+1;
            new(f(k+1),j)=new(f(k+1),j)+1;
            new(k,f(j+1))=new(k,f(j+1))+1;
            new(k,f(j-1))=new(k,f(j-1))+1;
          end 
       end % for j
   end  % for k 
   counter;
   sum(sum(new));
   if mod(counter,skip)==0
      image(new);
      colormap(prism(top));
      colorbar;
      title(['\fontsize{30}step = ',num2str(counter)]);
      pause(delay)
   end % if skip
   if old==new
     break
   end % if
   old=new;
 end % while 

   image(new);
   colormap(prism(top));
   colorbar;
   title(['\fontsize{30}step = ',num2str(counter)]);

   return;
