% PILES
% copyright by Robert J. Marks II
% "My code isn't pretty - but it works!"
% A pile of sand is emulated.  Pieces of sand are dropped in the center of a sandbox.  They
% can only stack so high.  Above a certain threshold, some of the sand falls onto an adjacent
% pixel where there is less sand.  This new pixel may exceed the stack height and is force to
% shed some of its sand onto adjacent pixels, etc.  An avalanche can occur.  This process follows
% the power law described in Mark Buchanan's book, "Ubiquity".
%
% PARAMETERS (must be changed in the m-file)
% NX, XY = dimensions of the sandbox
% generations = number of times sand is dropped
% grains = number of drops before imaging.  Used to quickly reach steady state.
% stack = how high the sand is stacked before it tumbles
% drop = number of grains of sand dropped each iteration
% flicks = number of times end movie is played.  If zero, no movie is recorded.
% After the iterations are over, you can plot the following:
%		plot(sandboxsize)   gives the number of grains of sand on the table as a function of iteration
% 		plot(fell) gives the number of grains that fall as a function of generation
%		plot(his) gives a log-log histogram of fell
close all 
clear all
%
% PARAMETERS
NX=25; NY=25;			% This is the size of the sand box 
generations=26101;		% The number of times sand is dropped
grains=1;				% Sand drops before display.  Should be less than 'generations'
stack=5;					% How big the stack gets before falling onto neighbors
drop=100;				% Number of grains dropped at a time
flicks=0;				% When the program is over, this is how many times the movie is played.
%							% If zero, no movie is made
mesher=1;				% If =1, display will be a mesh (looks best).  Otherwise, imagesc.
tries=500;				% Max iterations to avoid infinite loop
bins=100;				% Bins between low and high
genstart=100;			% Deletes points earlier during transient
%
% INITIALIZE
frame=0;
pile=zeros(NX,NY);
xo=ceil(NX/2);    yo=ceil(NY/2);
fell=ones(1,generations);
sandboxsize=zeros(1,generations);
%
% GENERATIONS
for g=2:generations
   x=xo;y=yo;
   pile(x,y)=pile(x,y)+drop;
   %%%%
   %Image loop
   if(g>grains)
      a=pile(1,1); b=pile(NX,NY); pile(1,1)=1; pile(NX,NY)= 6*(stack+drop);
      if(mesher==1)
         mesh(pile)
      else
         imagesc(pile)
      end
      title(g);xlabel(fell(g-1));ylabel(sandboxsize(g-1));
      frame=frame+1;
      if(flicks>0)
         M(frame)=getframe;
      else
         M=getframe;
      end
      pile(1,1)=a; pile(NX,NY)=b;
   end 
   iter=0;
   while(and(pile(xo,yo)>=stack,iter<tries))
      tryy=0;
      while(and(pile(x,y)>=stack,tryy<tries))
         xx=x; rr=rand; if(rr<1/3); xx=x+1; end; if(rr>2/3); xx=x-1; end
         yy=y; rr=rand; if(rr<1/3); yy=y+1; end; if(rr>2/3); yy=y-1; end
         if(or(or(xx>NX,xx<1),or(yy>NY,yy<1)))
            old=pile(x,y);
            pile(x,y)= max(pile(x,y)-drop,0); 
            fell(g)=fell(g)+pile(x,y)-old;
            xx=x; yy=y;
         elseif(and(pile(x,y)>pile(xx,yy),pile(x,y)>=stack))
            moved=max(0,floor((pile(x,y)-pile(xx,yy))/2)-1);
            pile(x,y)=pile(x,y)-moved;
            pile(xx,yy)=pile(xx,yy)+moved;
            fell(g)=fell(g)+moved;
            if(moved~=0)
               %%% Image loop
               if(g>grains)
                  a=pile(1,1); b=pile(NX,NY); pile(1,1)=1; pile(NX,NY)= 6*(stack+drop);
                  if(mesher==1)
                     mesh(pile)
                  else
                     imagesc(pile)
                  end
                  title(g);xlabel(fell(g-1));ylabel(sandboxsize(g-1));
                  frame=frame+1;
                  if(flicks>0)
                     M(frame)=getframe;
                  else
                     M=getframe;
                  end
                  pile(1,1)=a; pile(NX,NY)=b;
               end 
               %%%
            end
         end
         tryy=tryy+1;
         x=xx;y=yy;
      end
      iter=iter+1;
   end
   sandboxsize(g)=sum(sum(pile));
end
if(flicks>=1)  
   movie(M,flicks)
end
% Generating log-log histogram
histlow=min(fell)+1;				% Low Point of log-log histogram
histhi=max(fell)+1;				% High Point of log-log histogram
rati=histhi/histlow;
histog=ones(bins+1,1);
if(generations>=genstart)
   for g=genstart:generations
      for b=1:bins+1 
         if(and(histlow*(rati^((b-1)/bins))<fell(g),fell(g)<=histlow*(rati^(b/bins))))
            histog(b)=histog(b)+1;
         break
         end
      end
   end
end
figure
semilogy(histog); title('log-log histogram');
figure
plot(sandboxsize); title('sandbox size');
figure
plot(fell); title('fallen grains');