DiffusionLimitedAggregation.jl

using Plots
gr()

const N=60
const C=convert(Int,N/2)
grid=zeros(N,N)
grid[C,C]=1.0

const motion=[[1,0],[0,1],[-1,0],[0,-1]]

b(a)=mod((a-1),N)+1

function neighbor(a::Array{Int64,1})
    down=b.(a+[1,0])
    up=b.(a+[-1,0])
    right=b.(a+[0,1])
    left=b.(a+[0,-1])
    
    sum=0
    sum+=grid[down[1],down[2]]
    sum+=grid[up[1],up[2]]
    sum+=grid[right[1],right[2]]
    sum+=grid[left[1],left[2]]
    return sum
end

@time for t in 1:300
    initialPoint=rand([[1,rand(1:N)],[N,rand(1:N)],[rand(1:N),1],[rand(1:N),N]])
    p=initialPoint

    while neighbor(p)==0
        p.=b.(p+rand(motion))        
    end

    grid[p[1],p[2]]=1.0
end

heatmap(grid,
        aspect_ratio=:equal)