Julia code using Divide and Concur algorithm to pack circles in a circle

Circle packing

function dm(Pa, Pb, y, β, n)
    for i in 1:n
        if i%100 == 0
            print("Iteration $i/$n\n")
        end
        fay = (1-1/β)*Pa(y)+(1/β)*y
        fby = (1+1/β)*Pb(y)-(1/β)*y
        y = y+β*(Pa(fby)-Pb(fay))
#        @show(y)
    end

    fby = (1+1/β)*Pb(y)-(1/β)*y
    return Pa(fby)[1,:]
end

# N rows of width K
# Each row is a replica
function Pc(x)
    pc = repmat(mean(x,1),size(x)[1])
    return pc
end

function Pd(P, x)
#    @show("Pd.1")
#    @show(x)
    pd = x
    for i in 1:size(x)[1]
        pd[i, :] = P[i](x[i, :])
    end
    return pd
end

function apart(i, w, r, N, x)
#    @show("apart.1")
    d = 2*r
    z = x
#    @show(z, i)
    p = z[2*i-1:2*i] # centre of repulsion
#    @show("apart.2")
    for j in 1:N
        if j != i
            q = z[2*j-1:2*j]
            d_actual = norm(q-p)
            if d_actual < d
                unit = (d-d_actual)/d_actual*(q-p)
                z[2*j-1:2*j] += unit
            end
        end
    end
    return z
end

function in_circle(w, r, N, x)
#    @show("in_circle.1")
    z = x
    t = 0.5*w-r
    for j in 1:N
#    @show("in_circle.2")
        p = [0.5, 0.5]
        q = z[2*j-1:2*j]
#    @show("in_circle.3")
        d_actual = norm(q-p)
        if d_actual > t
#    @show("in_circle.4")
            unit = t/d_actual*(q-p)
            z[2*j-1:2*j] = p+unit
        end
    end
    return z
end

# For N disks we need
# 1. projection of everything back into box
# 2. N projections where all other disks get pushed away
# That's N+1 projections
# So N+1 rows, 2N columns
N = 14
r = 0.116 # disk radius
w = 1 # box size
P = [x -> apart(i, w, r, N, x) for i in 1:N]
push!(P, x -> in_circle(w, r, N, x))
y = repmat(rand((1, 2*N)), N+1)
@show(y)
y = dm(Pc, x -> Pd(P, x), y, 0.3, 100000)
#@show(w)
#@show(y)
#@show(2*r)
#@show(norm(y[1:2]-y[3:4]))
#@show(y[1]-r, y[3]-r, w-r-y[1], w-r-y[3])
#@show(y[2]-r, y[4]-r, w-r-y[2], w-r-y[4])

f = open("pack_circle.ps", "w")
write(f, "%!PS-Adobe-1.0\n")
write(f, "0 setlinewidth\n")
write(f, "100 200 translate\n")
write(f, "400 400 scale\n")
for i in 1:N
    write(f, "$(y[2*i-1]) $(y[2*i]) $r 0 360 arc closepath stroke\n")
end
write(f, "0.5 0.5 0.5 0 360 arc closepath stroke\n")
#write(f, "$(y[3]) $(y[4]) $r 0 360 arc closepath stroke\n")
write(f, "closepath\n")
write(f, "stroke\n")
close(f)