type SandPile
    M::Array{Int,2}
    n::Tuple{Int,Int}
    SandPile(a::Array{Int,2})=new(a,size(a))
end    

function avalanche(A::Array{Int,2},limit=3)
    while any(A.>limit) 
        R = CartesianRange(size(A))
        B=zeros(A)
        for i in R
            if A[i]>limit
                B[i]-=(limit+1)
                modPts=map(c->i+CartesianIndex(c),vcat([(z,0) for z=(-1,1)],[(0,z) for z=(-1,1)]))
                for j in modPts
                    (j in R) && (B[j]+=1)
                end
            end
        end
        A=A+B
    end
    A
end

import Base.+
function +(a::SandPile,b::SandPile) 
    a.n==b.n || error("Sandpiles must have the same dimenstion")
    c=a.M+b.M
    SandPile(avalanche(c))
end

sp_zero=SandPile([2 1 2
                  1 0 1
                  2 1 2])

sp_zero+sp_zero

SandPile([2 2 2
          2 2 2
          2 2 2])+sp_zero

s3=SandPile(ones(Int,3,3)*3)
ns3=ones(Int,3,3)*3
ns3[2,2]=1
s̄3=SandPile(ns3)
s3+s̄3

@show sp_rnd=rand(0:3,3,3)
new_sp=SandPile(avalanche(ones(Int,3,3)*3+sp_rnd))

new_sp+sp_zero