function [S,H,p] = sand(L,N,varargin)

% SAND
% AUTHOR
% INPUTS
% OPTIONS
% OUTPUTS
% REF
% See also

% initialise pile
% relax
% for Nnew {
% 	add grain
% 	relax
% 	}

addlib({'lib-misc'});

% DEFAULTS ---------------------------------------------------------------- %
Zc = 6;		% critical gradient
vis = 1;	% show pile
bnd = 1;	% add grains only to boundary
% ------------------------------------------------------------------------- %

% variable inputs
if nargin > 4, bnd = varargin{3}; end
if nargin > 3, vis = varargin{2}; end
if nargin > 2
	p = varargin{1};
else
	p = round(Zc*rand(L,L))+1;
end

[p,s] = relax(p,Zc,L,vis);
S = zeros(N,1);
H = zeros(N,1);

% boundary
bnd = [L:L:L^2 L*(L-1)+1:L^2-1];

% add sand
for n = 1:N
	if bnd		% only add to boundary
		iadd = choose(bnd,1);
		p(iadd) = p(iadd) + 1;
	else		% free add
		iadd = max(round(rand*L),1);
		jadd = max(round(rand*L),1);
		p(iadd,jadd) = p(iadd,jadd) + 1;
	end
	[p,s] = relax(p,Zc,L,vis);
	S(n) = s;
	H(n) = max(max(p));
	if mod(n,N/10)==0
		P = visprog(n,N,60,1);
		clc,disp('Progress:'); disp(P);
	end
end

% ----------------------------------------------------------------------- %
function [p,s] = relax(p,Zc,L,vis)

% find locations z(i,j) > Zc
% for each location {
% 	z(i,j) = z(i,j) - Zc
%	distribute Zc to surrounding cells (Moore 'hoods?)
% 	}

fall = [2 2 0
    	2 0 0
		0 0 0];

p0 = zeros(L+2,L+2);
p0(2:L+1,2:L+1) = p;
z0 = gradient(p0,L);
mask = zeros(L+2,L+2);
mask(2:L+1,2:L+1) = ones(L,L);
[I,J] = find(z0 > Zc);
s = 0;

while I
	gain = zeros(L+2,L+2);
	loss = zeros(L+2,L+2);
	for k = 1:length(I)
		i = I(k);
		j = J(k);
		gain(i-1:i+1,j-1:j+1) = gain(i-1:i+1,j-1:j+1) + fall;
		loss(i,j) = loss(i,j) - Zc;
	end
	p0 = max(p0 + gain + loss,0);
	s = s + sum(sum(p0 .* (1-mask)));
	z0 = gradient(p0,L);
	z0 = z0 .* mask;
	[I,J] = find(z0 > Zc);
	if vis
		bar3(flipud(p0));
		axis([0 L+2 0 L+2 0 L*Zc]),axis square,colormap(white),axis off;
		drawnow
	end
end
p = p0(2:L+1,2:L+1);

% ----------------------------------------------------------------------- %
function [z0] = gradient(p0,L)

x(:,:,1) = p0(2:L+1,1:L) - p0(2:L+1,2:L+1);
x(:,:,2) = p0(1:L,1:L) - p0(2:L+1,2:L+1);
x(:,:,3) = p0(1:L,2:L+1) - p0(2:L+1,2:L+1);
z = mean(x,3);
z0 = zeros(L+2,L+2);
z0(2:L+1,2:L+1) = -z;





