sand.jl

using ProgressMeter
using Plots
plotly(size=(500,500))

function updatesandpile!(d, i, j)
    n = get(d, (i,j), 0)
    if n >= 4
        d[(i-1,j)] = get(d, (i-1,j), 0) + 1
        d[(i+1,j)] = get(d, (i+1,j), 0) + 1
        d[(i,j-1)] = get(d, (i,j-1), 0) + 1
        d[(i,j+1)] = get(d, (i,j+1), 0) + 1
        d[(i,j)]   = get(d, (i,j), 0) - 4
        return true
    end
    return false
end

N = 100000

d = Dict([((0,0), N)])

changed = true
@showprogress for k in 1:2N
    for (i,j) in keys(d)
        changed = updatesandpile!(d, i, j)
    end
    if !changed
        break
    end
end

xs = [i for (i,j) in keys(d)]
ys = [j for (i,j) in keys(d)]
xmin = minimum(xs)
xmax = maximum(xs)
ymin = minimum(ys)
ymax = maximum(ys)

mat = zeros(xmax - xmin + 1, ymax - ymin + 1)

for ((i,j), v) in pairs(d)
    mat[i - xmin + 1, j - ymin + 1] = v
end

display(heatmap(mat, aspect_ratio=:equal))

Source: https://www.ams.org/journals/notices/201903/201903FullIssue.pdf p. 411-413