function [Z,S,LS,TE,LT] = sand_btw(L,N,varargin)

% SAND_BTW
%   Fractal implementation of sand-pile model. In this case, a simple
%   gradient matrix is considered (no physical 'height' as such) following
%   more closely the initial model. Sand is added one grain at a time
%   to the very centre of the 'pile' and relaxation applied symmetrically
%   in the Moore neighbourhood. This model yields interesting Total Energy
%   dissapation curves, having a power-law component (alpha ~ -1.0) giving
%   way to an exponential decay component.
% AUTHOR
%   (c) S.Angus, January 2007 (UNSW)
% INPUTS
%   (L,N) the side length L and number of grains to add N.
% OPTIONS
%   (...,VIS) VIS=1 will show changing structure with each update using
%     PCOLOR (default VIS=0).
%   (...,Z0) allows an initial state for the gradient matrix to be
%     used -- useful when 'warming-up' the system (watch S with time)
%     so that data can be taken thereafter.
% OUTPUTS
%   [Z,S,LS,TE,LT] The final gradient matrix Z, the average Z score S,
%     Critical Response statistics:
%       LS - linear size of the avalanche
%       TE - total energy (sum of sites that had to be relaxed)
%       LT - lifetime total number of relaxation events needed
% Depends: lib-misc
% See also: PCOLOR

addlib({'lib-misc'});

% DEFAULTS ---------------------------------------------------------------- %
Zc = 4;		% critical gradient
vis = 0;	% show pile
bnd = 1;	% add grains only to boundary
% ------------------------------------------------------------------------- %
% initialize
Z = zeros(L,L);
S = zeros(N,1);
LS = zeros(N,1);
TE = zeros(N,1);
LT = zeros(N,1);

% variables
if nargin > 2, vis = varargin{1}; end
if nargin > 3, Z = varargin{2}; end

% add sand
for n = 1:N
	iadd = round(L/2);
	Z(iadd,iadd) = Z(iadd,iadd) + 1;
	[Z,ls,te,lt] = relax(Z,Zc,L,vis);
	S(n) = sum(sum(Z));
	LS(n) = ls;
	TE(n) = te;
	LT(n) = lt;
	if mod(n,N/10)==0
		P = visprog(n,N,60,1);
		clc,disp('Progress:'); disp(P);
	end
end

% ----------------------------------------------------------------------- %
function [Z,ls,te,lt] = relax(Z,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 = [0  1  0
		1 -4  1
		0  1  0];

z0 = zeros(L+2,L+2);
z0(2:L+1,2:L+1) = Z;
mask = zeros(L+2,L+2);
mask(2:L+1,2:L+1) = ones(L,L);
[I,J] = find(z0 > Zc);
Is = []; Js = [];
s = 0;
c1 = 0;

while I
	c1 = c1 + 1;
	Is = [Is;I];
	Js = [Js;J];
	for k = 1:length(I)
		i = I(k);
		j = J(k);
		z0(i-1:i+1,j-1:j+1) = z0(i-1:i+1,j-1:j+1) + fall;
	end
	z0 = z0 .* mask;
	[I,J] = find(z0 > Zc);
	if vis
		pcolor(flipud(z0));
		axis square,colormap(hsv),axis off;
		drawnow
	end
end
[Th,R] = cart2pol(Is-mean(Is),Js-mean(Js));
ls = mean(R);
te = length(Is);
lt = c1;
Z = z0(2:L+1,2:L+1);




