I am currently trying to solve and understand a problem about packing the maximum of squares into a circle. What I am trying to do, is to calculate the maximum number of variable sized squares (with fixed side) into a circle of variable diameter. For example the result you get at WolframAlpha typing 'pack squares of side 50mm into a circle of diameter 300mm':

Densest known packing:

enter image description here

Furthermore, I would like to be able to add a variable for horizontal and vertical spacing between each square. I have found this to be exactly what I am trying to understand and build.


1 Answer 1

(* i, j, -n ...0 ...n*)
(* calculate a square center with vertical and horizontal spacings*)
center[i_, j_, s_, spx_, spy_] := {i (s + spx), j (s + spy)}
(* calculate the vartices of a square given the center*)
verts[{cx_, cy_}, s_] := {cx, cy} + # s/2 & /@ Tuples[{-1, 1}, 2]
(* calculate all squares to consider, given a circle center*)
(* note: it's overestimating, but it doesn't matter *)
allSquares[cirCenterX_, cirCenterY_, circR_, s_, spx_, spy_] := 
                         N@center[i, j, s, spx, spy] /. Solve[(#.#) &@({cirCenterX, cirCenterY} - 
                                  center[i, j, s, spx, spy]) < (circR + s + spx + spy)^2, {i, j}, 
(* calculate all verttices to consider, given a circle center*) 
allVerts[cirCenterX_, cirCenterY_, circR_, s_, spx_, spy_] := 
                     N@verts[#, s] & /@ allSquares[cirCenterX, cirCenterY, circR, s, spx, spy]
(* test if a square is inside a given circle *)
testVertInsideCirc[cirCenterX_, cirCenterY_, circR_, vert_] := 
                     Norm[{cirCenterX, cirCenterY} - vert] <= circR

(* set a problem*)
s = 10; (* square side *)
r = 50; (* circle radius *)
spx = 2/10; (*horiz spacing*)
spy = 3/10; (*vert spacing*)
(* calculate all possible involved squares for a given problem*)
allvs = allVerts[0, 0, r, s, spx, spy];
(* Solve*)
nm = NMaximize[{Tr[ Boole /@  And @@@ Map[testVertInsideCirc[ccx, ccy, r, #] &, N@allvs, {2}]],
               {-s/2 <= ccx <= s/2, -s/2 <= ccy <= s/2}}, {ccx, ccy},
               Method -> "DifferentialEvolution"]

{59., {ccx -> -3.33338, ccy -> -0.0896481}}

Equivalent problem solved with your linked reference:

Mathematica graphics


Drawing our solution:

With[{cx = ccx /. nm[[2]], cy = ccy /. nm[[2]]}, 
     Graphics[{EdgeForm[Directive[Thick, Black]], 
              Flatten@Transpose@{And @@@ 
                     Map[testVertInsideCirc[cx, cy, r, #] &, N@allvs, {2}] /. 
                                                      {True -> Blue, False -> Yellow}, 
                                    Rectangle @@@ (allvs[[All, {1, 4}]])}, 
              Red, Thick, Circle[{cx, cy}, r], PointSize[Large], Point[{cx, cy}]}]]

Mathematica graphics

  • $\begingroup$ This is awesome! But I have tried playing with the spacing settings a little. As I go smaller, calculations fail. $\endgroup$
    – Alex
    Aug 27, 2013 at 18:49
  • $\begingroup$ @Alex Could you please try it again? I corrected a minor bug. If you still have problems please post the example so I can look for the problem $\endgroup$ Sep 7, 2013 at 19:14
  • $\begingroup$ I seem to get calculation errors, when I use s=1,r=150,spx=spy=1/1000 for example. I am not quite sure, but it seems, that the s variable is causing this? $\endgroup$
    – Alex
    Sep 19, 2013 at 17:07

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.