function sand1
clf; hold on;
h = zeros(1,20);

% iterate for a while to get to steady state (almost)
[h hist] = iterate (h, 500);

% iterate and collect a histogram of landslides
[h hist] = iterate (h, 100);

% plot the histogram
y = 0:length(h);
bar (y, hist);


function [h, hist] = iterate (h, n)
% simulate n time steps.  Return the updates vector h
% and a histogram of landslide sizes
hist = zeros (1,length(h)+1);
for i=1:n
    h(1) = h(1) + 1;
    [h count] = topple (h);
    index = count + 1;
    hist(index) = hist(index) + 1;
end


function [h, count] = topple (h)
% simulate a single time step by making two passes
% through the vector h.  On the first pass, identify
% locations that are unstable.  On the second pass,
% topple each unstable location.  Return the updated vector.
i1 = [2:length(h) 1];
m = h - h(i1);

count = 0;
for i=1:length(m)-1
    if m(i) > 1
        h(i) = h(i) - 1;
        h(i+1) = h(i+1) + 1;
        count = count + 1;
    end
end
h(length(h)) = 0;