%sand?pile
clear?all
clf
nx=52;?%must?be?divisible?by?4
ny=100;
Pbridge?=?.05;
z=zeros（nxny）;
o=ones（nxny）;
sand?=?z?;
sandNew?=?z;
gnd?=?z?;
gnd（1:nxny-3）=1?;?%?the?ground?line
gnd（nx/4:nx/2+4ny-15）=1;?%the?hole?line
gnd（nx/2+6:nxny-15）=1;?%the?hole?line
gnd（nx/4?ny-15:ny）?=?1;?%side?line
gnd（3*nx/4?1:ny）?=?1?;
imh?=?image（cat（3z‘sand‘gnd‘））;
set（imh?‘erasemode‘?‘none‘）
axis?equal
axis?tight
for?i=1:1000
p=mod（i2）;?%margolis?neighborhood
sand（nx/2ny/2）?=?1;?%add?a?grain?at?the?top
%upper?left?cell?update
xind?=?[1+p:2:nx-2+p];
yind?=?[1+p:2:ny-2+p];
vary?=?rand（nxny）<.95?;
vary1?=?1-vary;
sandNew（xindyind）?=?...
gnd（xindyind）.*sand（xindyind）?+?...
（1-gnd（xindyind））.*sand（xindyind）.*sand（xindyind+1）?.*?...
（sand（xind+1yind+1）+（1-sand（xin