# Randomly packing spheres of fixed radius within a cube

How can I have Mathematica randomly place spheres in a cube so they won't overlap? The cube is $20 \times 20 \times 20$, and the spheres have a radius of $0.7$.

-
What have you tried doing? –  rm -rf Oct 12 '12 at 1:30
Interesting problem. Could you explain the context of the problem? Just curiosity? –  belisarius Oct 12 '12 at 1:37
How many spheres? –  Mark McClure Oct 12 '12 at 1:38
I'm trying to solve a packing density problem with mathematica. Eventually, I want to fill the cube to a volume fraction of 1/3 of the spheres. Right now I'm struggling to plot the cube and hard spheres. –  Jen Oct 12 '12 at 1:43
So you are just trying to plot right now? You can do cube = {Opacity[.3], Cuboid[{0, 0, 0}, {20, 20, 20}]}; spheres = Table[{Black, Sphere[RandomReal[{.35, 20 - .35}, 3], .7]}, {j, 1, 10}]; Graphics3D[{cube, spheres}] –  David Slater Oct 12 '12 at 1:54

Here's a brute force and soporifically slow method.

spheres = {};
Dynamic[Length[spheres]]


Then create a new random sphere, and if it's not too close to any existing sphere add it to the list:

While[Length[spheres] < 1856, s = RandomReal[{0.7, 20 - 0.7}, 3];
If[And @@ (Norm[# - s] > 1.4 & /@ spheres), AppendTo[spheres, s]]];


According to my calculations, you need 1856 spheres of radius 0.7 to fill 1/3 of the volume of a 20 x 20 x 20 cube.

To display:

cube = {Opacity[0.3], Cuboid[{0, 0, 0}, {20, 20, 20}]};
Graphics3D[{cube, Sphere[#, 0.7] & /@ spheres}, Boxed -> False]


As I stated at the beginning, this is very slow and inefficient. The above image is from when I gave up after 1500 spheres. You can Alt. and resume at your leisure. It's entirely possible, what with the randomness, that you will never fit the 1856 spheres in.

-
It's a good way to get the OP started. An extension would be to allow for wiggling the locations. One could also treat it as like an electromagnetics problem, for instance — think of the spheres as all having unit charge and reposition the spheres just as they would when another charged sphere is added to the existing bunch, subject to the constraint that none leave the box. –  rm -rf Oct 12 '12 at 3:36
The alternate type approach is to randomly place spheres so that each exactly touches three others. Much tougher to implement but it will get you ultimately to higher density (50+%). –  george2079 Oct 12 '12 at 4:01
A relatively simple scheme could be to select a random sphere in the list, then creating a new sphere that touches it in a random direction. This got me to 1856 spheres in a few minutes. –  wxffles Oct 12 '12 at 4:54
Why is the standing assumption that the metric space is Euclidean (L2-norm)? –  alancalvitti Oct 12 '12 at 6:22

I will present another brute-force method for densely packing spheres that works by randomly depositing the spheres in successive "layers". The first layer is built by solving the two-dimensional analog of the problem (that is, packing disks in a square). Having constructed the first layer, successive layers are built up by dropping new spheres, checking that the incoming spheres do not intersect with spheres in previous layers (as well as with other spheres in the same layer). For the second layer, the new spheres are compared against spheres in the first layer, while for the third and succeeding layers, the new sphere is tested for intersections with the two previous layers. All this is done up until we hit the top of the cube.

Here is a Mathematica implementation of the idea presented above:

sphrPak[s_, r_, its_Integer: 2000] :=
Module[{j, k, lap, lap2, layer, layers, p0, pc, sphrs, z},

(* pack disks into square *)

(* pick first disk *)
While[
p0 = RandomReal[{r, s - r}, 2];
And @@ Thread[Map[Min[#, s - #] &, p0] < r]];

(* successively pack disks into square *)
layer = {p0}; k = 0;
While[k < its,
pc = RandomReal[{r, s - r}, 2]; (* candidate disk *)
If[And @@ Thread[Map[EuclideanDistance[pc, #] &, layer] >= 2 r],
k = 0; AppendTo[layer, pc];,
++k];
];

(* build and store first layer *)
layers = {PadRight[layer, {Length[layer], 3}, r]};

(* build successive layers *)
j = 2;
While[True,

(* previous layers *)
lap = layers[[j - 1]]; If[j > 2, lap2 = layers[[j - 2]]];

(* find "seed" sphere for layer *)
p0 = {}; k = 0;
While[k < its,
z = RandomReal[{j, j + 1} r];
pc = Append[RandomReal[{r, s - r}, 2], z]; ++k; (* candidate seed *)
If[z <= s - r &&
(And @@ Thread[Map[EuclideanDistance[pc, #] &, lap] >= 2 r]),
If[j > 2,
If[And @@ Thread[Map[EuclideanDistance[pc, #] &, lap2] >= 2 r],
p0 = pc; Break[]],
p0 = pc; Break[]]]];

If[p0 === {}, Break[]]; (* no seed found *)

(* start building new layer *)
layer = {p0}; k = 0;
While[k < its,
z = RandomReal[{j, j + 1} r];
pc = Append[RandomReal[{r, s - r}, 2], z]; (* candidate point *)
If[z <= s - r &&
(And @@ Thread[Map[EuclideanDistance[pc, #] &, layer] >= 2 r]) &&
(And @@ Thread[Map[EuclideanDistance[pc, #] &, lap] >= 2 r]),
If[j > 2,
If[And @@ Thread[Map[EuclideanDistance[pc, #] &, lap2] >= 2 r],
k = 0; AppendTo[layer, pc], ++k],
k = 0; AppendTo[layer, pc]], ++k];
];
AppendTo[layers, layer];
++j];
Print[StringForm["1 spheres generated.", Length[sphrs = Join @@ layers]]];
Graphics3D[Sphere[sphrs, r], Axes -> Automatic, PlotRange -> {{0, s}, {0, s}, {0, s}}]]


Try it out:

BlockRandom[SeedRandom[2012, Method -> "MersenneTwister"]; (* for reproducibility *)
sphrPak[20, 0.7]]
1996 spheres generated.


BlockRandom[SeedRandom[1024, Method -> "MersenneTwister"];
sphrPak[20, 0.7]]
2041 spheres generated.


-
being nit picky, if i understand this scheme correctly it seems you are introducing a certain non-randomness in that you have a layered structure with constant z planes. That may be important depending on what you want to do with the result. Perhaps throw in a small random z within each layer. Might need to check two layers below if you do that though. –  george2079 Oct 15 '12 at 12:25
No @george, the heights at each layer (except the first) are not constant; that is, unless z = RandomReal[{j, j + 1} r] is doing something peculiar that you know and I don't. Yes, I check two layers below as well; that's what the If[j > 2,(* stuff *)] sections of the seeding and layer-building parts are for. –  Ｊ. Ｍ. Oct 15 '12 at 12:30
sorry..in my head that was outside the loop. Now i'm not sure that seeding the layer is necessary, you could find the mean z for each layer and move up a fixed amount. In the end i expect that seeding is a small part of the time though. –  george2079 Oct 15 '12 at 13:08

Here's a brute-force method (based on a bounded random walk) that attempts to generate n nonintersecting spheres of radius r within a cube of side length s, as well as reporting the actual number of spheres it managed to generate after lim unsuccessful iterations:

sphrPak[s_, r_, n_Integer?Positive, lim_Integer?Positive] :=
Module[{sh = s/2, k, p0, pc, sphList},

(* generate initial sphere within box *)
While[
p0 = RandomReal[{-sh, sh}, 3];
And @@ Thread[Map[Min[sh - #, sh + #] &, p0] < r]];
sphList = {p0};

k = 0;
While[Length[sphList] < n && k < lim,
(* center for new sphere chosen in random direction *)
pc = p0 + 2 r Normalize[RandomVariate[NormalDistribution[], 3]];
If[(* is sphere within the cube? *)
PolyhedronData["Cube", "RegionFunction"] @@
ScalingTransform[ConstantArray[1/(s - 2 r), 3]][pc] &&
(* does the sphere not overlap with other spheres? *)
(And @@ Thread[Map[EuclideanDistance[pc, #] &, sphList] >= 2 r]),
k = 0; AppendTo[sphList, pc]; p0 = pc,
(* else *) ++k];
];

If[k == lim, Print[StringForm["Only 1 spheres were generated.", Length[sphList]]]];

Graphics3D[Sphere[sphList, r], Axes -> Automatic,
PlotRange -> {{-sh, sh}, {-sh, sh}, {-sh, sh}}]];

sphrPak[s_, r_, n_Integer] := sphrPak[s, r, n, Quotient[2 n, 3]]


Try it out:

BlockRandom[SeedRandom[1024, Method -> "MersenneTwister"]; (* for reproducibility *)
sphrPak[20, 0.7, 1500]]
Only 725 spheres were generated.


Of course, a different attempt might yield a result with more spheres within the cube.

More attempts:

BlockRandom[SeedRandom[4092, Method -> "MersenneTwister"];
sphrPak[20, 0.7, 1500]]
Only 805 spheres were generated.


BlockRandom[SeedRandom[2012, Method -> "MersenneTwister"];
sphrPak[20, 0.7, 1500]]
Only 932 spheres were generated.


Here's a variation that attempts to do more than one bounded random walk, where the starting point of the $n$th random walk is a randomly chosen point from the $(n-1)$-th random walk:

sphrPak2[s_, r_, n_Integer?Positive, walks_Integer?Positive, its_Integer?Positive] :=
Module[{sh = s/2, j, k, p0, pc, sphList},
While[
p0 = RandomReal[{-sh, sh}, 3];
And @@ Thread[Map[Min[sh - #, sh + #] &, p0] < r]];
sphList = {p0}; j = 0;
While[Length[sphList] < n && j < walks,
k = 0;
While[Length[sphList] < n && k < its,
pc = p0 + 2 r Normalize[RandomVariate[NormalDistribution[], 3]];
If[(PolyhedronData["Cube", "RegionFunction"] @@
ScalingTransform[ConstantArray[1/(s - 2 r), 3]][pc]) &&
(And @@ Thread[Map[EuclideanDistance[pc, #] &, sphList] >= 2 r]),
k = 0; AppendTo[sphList, pc]; p0 = pc,
++k];
];
If[Length[sphList] < n,
p0 = RandomChoice[Map[EuclideanDistance[#, sphList[[1]]] &,
sphList] -> sphList]];
++j;
];
If[j == walks, Print[StringForm["Only 1 spheres were generated.",
Length[sphList]]]];
Graphics3D[Sphere[sphList, r], Axes -> Automatic,
PlotRange -> {{-sh, sh}, {-sh, sh}, {-sh, sh}}]]


Try it out:

BlockRandom[SeedRandom[85, Method -> "Rule50025CA"];
sphrPak2[20, 0.7, 1500, 1, 1000]]
Only 662 spheres were generated.


BlockRandom[SeedRandom[85, Method -> "Rule50025CA"];
sphrPak2[20, 0.7, 1500, 20, 1000]]
Only 1335 spheres were generated.


-
Why do you need to check if the center AND the hull are inside the cube? –  belisarius Oct 12 '12 at 12:14
@bel: I'm exploiting the fact that And[] does short-circuit evaluation; if it fails the easy test, then the other tests are not touched. On the other hand, I suppose I can merge those first two tests; wait... –  Ｊ. Ｍ. Oct 12 '12 at 12:19