How to pack spheres in python?
<p>I am trying to model random closed packing spheres of uniform size in a square using python. <strong>And the spheres should not overlap</strong>  but I do not know how to do this</p>

<p>I have so far:
<a href="http://i.stack.imgur.com/fwScD.png" rel="nofollow"><img src="http://i.stack.imgur.com/fwScD.png" alt="enter image description here"></a></p>

<p>Code:</p>

<pre><code>import random, math, pylab

def show_conf(L, sigma, title, fname):
    pylab.axes()
    for [x, y] in L:
        for ix in range(-1, 2):
            for iy in range(-1, 2):
                cir = pylab.Circle((x + ix, y + iy), radius=sigma,  fc='r')
                pylab.gca().add_patch(cir)
    pylab.axis('scaled')
    pylab.xlabel('eixo x')
    pylab.ylabel('eixo y')
    pylab.title(title)
    pylab.axis([0.0, 1.0, 0.0, 1.0])
    pylab.savefig(fname)
    pylab.close()

L = []
N = 8 ** 2

for i in range(N):
    posx = float(random.uniform(0, 1))
    posy = float(random.uniform(0, 1))
    L.append([posx, posy])

print L

N = 8 ** 2
eta = 0.3
sigma = math.sqrt(eta / (N * math.pi))
Q = 20
ltilde = 5*sigma

N_sqrt = int(math.sqrt(N) + 0.5)


titulo1 = '$N=$'+str(N)+', $\eta =$'+str(eta)
nome1 = 'inicial'+'_N_'+str(N) + '_eta_'+str(eta) + '.png'
show_conf(L, sigma, titulo1, nome1)
</code></pre>

<p>This is a very hard problem (and probably <strong>np-hard</strong>). There should be a lot of ressources available.</p>

<p>Before i present some more general approach, check out <a href="https://en.wikipedia.org/wiki/Circle_packing_in_a_square" rel="nofollow">this wikipedia-site</a> for an overview of the currently best known packing-patterns for some N (N circles in a square).</p>

<p>You are lucky that there is an existing circle-packing implementation in python (<strong>heuristic!</strong>) which is heavily based on modern optimization-theory (<a href="https://ocw.mit.edu/courses/sloan-school-of-management/15-097-prediction-machine-learning-and-statistics-spring-2012/projects/MIT15_097S12_proj5.pdf" rel="nofollow">difference of convex-functions</a> + <a href="http://papers.nips.cc/paper/2125-the-concave-convex-procedure-cccp.pdf" rel="nofollow">Concave-convex-procedure</a>).</p>

<ul>
<li>The method used is described <a href="https://stanford.edu/~boyd/papers/dccp.html" rel="nofollow">here</a> (academic paper &amp; link to software; 2016!)</li>
<li>The software package used is <a href="https://github.com/cvxgrp/dccp" rel="nofollow">here</a>

<ul>
<li>There is an example directory with <code>circle_packing.py</code> (which is posted below together with the output)</li>
</ul></li>
<li>The following example also works for circles of different shapes</li>
</ul>

<h3>Example taken from the above software-package (example by Xinyue Shen)</h3>

<pre><code>__author__ = 'Xinyue'
from cvxpy import *
import numpy as np
import matplotlib.pyplot as plt
import dccp

n = 10
r = np.linspace(1,5,n)

c = Variable(n,2)
constr = []
for i in range(n-1):
    for j in range(i+1,n):
        constr.append(norm(c[i,:]-c[j,:])&gt;=r[i]+r[j])
prob = Problem(Minimize(max_entries(max_entries(abs(c),axis=1)+r)), constr)
#prob = Problem(Minimize(max_entries(normInf(c,axis=1)+r)), constr)
prob.solve(method = 'dccp', ccp_times = 1)

l = max_entries(max_entries(abs(c),axis=1)+r).value*2
pi = np.pi
ratio = pi*sum_entries(square(r)).value/square(l).value
print "ratio =", ratio
print prob.status

# plot
plt.figure(figsize=(5,5))
circ = np.linspace(0,2*pi)
x_border = [-l/2, l/2, l/2, -l/2, -l/2]
y_border = [-l/2, -l/2, l/2, l/2, -l/2]
for i in xrange(n):
    plt.plot(c[i,0].value+r[i]*np.cos(circ),c[i,1].value+r[i]*np.sin(circ),'b')
plt.plot(x_border,y_border,'g')
plt.axes().set_aspect('equal')
plt.xlim([-l/2,l/2])
plt.ylim([-l/2,l/2])
plt.show()
</code></pre>

<h3>Output</h3>

<p><a href="http://i.stack.imgur.com/Ryupg.png" rel="nofollow"><img src="http://i.stack.imgur.com/Ryupg.png" alt="enter image description here"></a></p>

<h3>Modification for your task: equal-sized circles</h3>

<p>Just replace:</p>

<pre><code>r = np.linspace(1,5,n)
</code></pre>

<p>With:</p>

<pre><code>r = [1 for i in range(n)]
</code></pre>

<h3>Output</h3>

<p><a href="http://i.stack.imgur.com/x3cBf.png" rel="nofollow"><img src="http://i.stack.imgur.com/x3cBf.png" alt="enter image description here"></a></p>

<h3>Fun example with 64 circles (this will take some time!)</h3>

<p><a href="http://i.stack.imgur.com/Th3zK.png" rel="nofollow"><img src="http://i.stack.imgur.com/Th3zK.png" alt="enter image description here"></a></p>

python optimization mathematical-optimization packing