**TL;DR** Yes, this can be done!

If you read the article "Hexagonal circle packings and Doyle spirals" by Leys, you will see that for a choice `p`

and `q`

, we need to find the complex values `A`

, `B`

and `r`

. For that purpose, we can steal this part from the demonstration you linked:

```
doyle[pi_, qi_] := Module[{p = pi, q = qi, s, t, r},
r[s_, t_, p_,
q_] := (s^2 + s^(2 p/q) -
2 s^((p + q)/q) Cos[(2 \[Pi] + p t - q t)/q])/(s + s^(p/q))^2;
{s, t} = {s, t} /. FindRoot[
{r[s, t, 0, 1] - r[s, t, p, q] == 0,
r[s, t, 0, 1] - r[s^(p/q), (p*t + 2. Pi)/q, 0, 1] == 0},
{{s, 1.5}, {t, 0}}];
{s*Exp[I*t], s^(p/q)*Exp[I (p*t + 2*Pi)/q], Sqrt[r[s, t, 0, 1]]}
]
{a, b, r} = doyle[2, 12]
```

Now we have the centers for both additional complex circles. Knowing a bit of complex analysis, one understands that for creating all packed circles, we only need to iteratively multiply by e.g. `a`

. So we could write down a function that does this multiplication. I'm keeping complex values until the final visualization, where we use the real part for x and the imaginary part for y:

```
iterate[a_, b_, j_, n_] := Module[{start = b^j},
Table[a^i*start, {i, Range[-n, n]}]
]
Graphics[Circle[ReIm[#], r*Abs[#]] & /@ iterate[a, b, 0, 3]]
```

This shows the $0th$ spiral of packed circles, $3$ circles *inward* to our base circle and $3$ circles outward.
To create the complete plane, we have to create $12$ columns, since `q`

was $12$:

```
toCirle[z_, r_] := Disk[{Re[z], Im[z]}, Abs[z]*r];
pack = Table[iterate[a, b, j, 5], {j, 12}];
gr = Graphics[{EdgeForm[Black],
Map[{RandomColor[], toCirle[#, r] & /@ #} &, pack]},
PlotRange -> {{-10, 10}, {-10, 10}}]
```

Unfortunately, that is not enough, because the artist chose to use the logarithmic spirals through the circle centers as guide to divide each circle into different parts. In order to do this, we need to go further.

Let us make a cut here and divide the following into small sections where we look at the details. These details will be important for the overall approach

## Connection between `a`

, `b`

, `r`

and the circles

As pointed out in the article, the complex numbers `a`

and `b`

are the *generators* of the circles. This means, all circle centers can be obtained by repeated multiplication. The base circle with the center `{1,0}`

is given by $a^0\cdot b^0$ which is 1 (meaning re=1 and im=0).

Now each multiplication by `a`

or by `b`

gives the center of the circle that is next to the base circle. So `1*a*a=a^`

gives the second circle in the direction of `a`

. `a^2*b`

shifts this last circle in the direct of `b`

.

Note, that even `a^(-3)`

is perfectly OK and gives the 3rd circle in the oposite direction. These are the small circles the fill the center. OK, one `Manipulate`

says more than a thousand words. Let us create a dynamic table of all circles in a range for `a^i*b^j`

. Note that, as pointed out in the article, the correct radius for each circle is `Abs[a^i*b^j]*r`

where `r`

is the radius we got from the solution of `doyle`

.

```
Manipulate[
Graphics[{EdgeForm[Black], FaceForm[Gray], Table[
toCirle[a^i*b^j, r],
{i, i1, i2}, {j, j1, j2}]}, PlotRange -> {{-10, 10}, {-10, 10}}],
{{i1, 0}, -5, 3, 1},
{{i2, 0}, i1, 5, 1},
{{j1, 0}, -5, 3, 1},
{{j2, 0}, j1, 5, 1}
]
```

## Circles and spirals

We have seen that we can go from one circle to any neighbour by increasing (or decreasing) either `i`

or `j`

by 1. But what if we don't make *jump* from say `a^3`

to `a^4`

but a smooth transition? Well, the function for such a thing is easy because `a^0=1`

and `a^1=a`

, so we can make a function `a^3*a^t`

and let `t`

run from 0 to 1.

```
Show[
Graphics[{
FaceForm[Gray],
toCirle[#, r] & /@ {a^3, a^4}}
],
ParametricPlot[ReIm[a^3*a^t], {t, 0, 1}, PlotStyle -> White]
]
```

This looks very much like the spirals that were used to divide the circles in the original art. So it seems if we pick out the center any circle next to our base circle, we can create spiral functions that go through the circles. Note that the approach of shifting a spiral to its neighbouring spiral is similar to shifting circles. Here is an example:

```
Show[
gr,
ParametricPlot[Table[ReIm[b^i a^t], {i, 12}], {t, -10, 10},
PlotRange -> {{-10, 10}, {-10, 10}}, PlotStyle -> White],
ParametricPlot[Table[ReIm[a^j b^t], {j, 5}], {t, -10, 10},
PlotRange -> {{-10, 10}, {-10, 10}}, PlotStyle -> White]
]
```

## Spiral functions inside circles

For our later approach, I want to be able to draw the spiral only *inside* a circle. As we have seen, going from `t=0`

to `t=1`

will connect the centers of the circles. This is not what we want. We want values for `t`

that start and end with the circle. Let's make the plot we did earlier again, but use values for `t`

between -1/3 and 1/3

OK, that looks promising. Remember, we know the center of this circle with `a^3`

and we know its radius with `Abs[a^3]*r`

. What are the bounds where our spiral is exactly on the radius? Let us ask `FindRoot`

:

```
tb = t /. FindRoot[Abs[1 - a^t] - r, {t, #}] & /@ {-1/3, 1/3}
(* {-0.565183, 0.433533} *)
```

But wait! I haven't used `a^3`

at all! Correct. The good thing is that the bounds for the circles apply to each circle of the same spiral. Therefore I'm using the next neighbour of the base circle which is `a`

for `FindRoot`

. Look here:

```
Show[
Graphics[{
FaceForm[Gray],
toCirle[#, r] & /@ {a^3, a^4}}
],
ParametricPlot[ReIm[a^3*a^t], {t, tb[[1]], tb[[2]]},
PlotStyle -> White]
]
```

## What spirals did the artist use?

As it turns out she used the spirals of the following direct neighboring circles of the base circle:

```
spoints = {a*b^-1, a, b}
(* {1.46301 - 0.54185 I, 1.67036 + 0.343254 I, 0.927594 + 0.578172 I} *)
```

Let's make a small function that calculates their bounds and returns them with a spiral function. The spiral function will directly incorporate the `i`

and `j`

so that we can easily draw it on every circle we like

```
spiral[pt_] := Module[{t1, t2},
{t1, t2} =
Block[{t},
t /. FindRoot[Abs[1 - pt^(t)] - r, {t, #}] & /@ {-1/3, 1/3}];
{t1, t2, Function[{i, j, t}, a^i*b^j*pt^t]}
]
```

Now let's plot these 3 spirals inside our base circle `{1,0}`

```
Show[{
Graphics[Circle[{1, 0}, r]],
ParametricPlot[ReIm@#3[0, 0, t], {t, #1, #2}] & @@@ spiral /@ spoints
}]
```

Now, we can calculate the points of the spirals inside each circle, we have the radius of each circle and through the spiral's start and end points, we have 6 points on each circle.

## Creating polygons points for the parts of a circle

For each cake-part of a circle, we can now create a polygon by

- starting in the center
- creating points along a spiral outwards to the circle boundary
- going counterclockwise along the circle to the endpoint of the next spiral
- create points along this next spiral from the outer point to the center

However, one tiny point is missing. How do we create points along the circle from when we go from one spiral point to the next. That is not as hard as it sounds.

Assume you have two (complex) points that lie on a circle around a center. Then you can subdivide them and create arbitrarily many points between them that all lie on the circle.

```
circle[z1_, z2_, cent_] :=
Module[{zz1 = z1 - cent, zz2 = z2 - cent, r},
r = Abs[Mean[{zz1, zz2}]];
# + cent & /@
Nest[Riffle[#,
Function[zz, With[{m = Mean[zz]}, m/Abs[m]*Abs[zz1]]] /@
Partition[#, 2, 1]] &, {zz1, zz2}, 5]
]
```

Having this, we can create the points for all cake-parts of circle `i`

, `j`

defined by the provided `spirals`

that divide the circle:

```
createCircleParts[spirals_, i_, j_] :=
Module[{center, outward, inward},
outward = Table[#3[i, j, t], {t, 0, #2, #2/10.}] & @@@ spirals;
inward = Table[#3[i, j, t], {t, 0, #1, #1/10.}] & @@@ spirals;
center = a^i*b^j;
{i, j,
Join[#1[[;; -2]], circle[#1[[-1]], #2[[-1]], center],
Reverse[#2[[;; -2]]]] & @@@
Partition[Join[outward, inward, {First[outward]}], 2, 1]}
]
```

The function returns `{i,j, {part1, part2, ...}}`

and we will use `i`

and `j`

later for the coloring as it gives us information about the position of the circle.

To test this function, let us see what happens with the circle `i=1`

and `j=2`

:

```
Graphics[{RandomColor[], Polygon[ReIm[#]]} & /@
Last@createCircleParts[spiral /@ spoints, 1, 2]
]
```

## Coloring of circles

For one circle we have the information `i`

, `j`

which encodes the global position and of course we have `n`

cake-parts. An easy way would be to provide a list of colors and select a color depending on the information we have.

I could not really find a pattern in the coloring of the artists image, so lets keep it simple but let us use equivalent colors:

```
cols = {Black, RGBColor[0.078, 0.71, 0.964], Orange, Red, Darker@Green, Purple};
colorCircleParts[{i_, j_, parts_}, col_List] :=
Table[{col[[Mod[i + j + n, Length[col]] + 1]],
Polygon[ReIm@parts[[n]]]}, {n, Length[parts]}]
```

## Putting everything together

The last thing we need to do is to create a table containing the circles and their parts for a range of `i`

and `j`

values. Then we color the circle parts and display them:

```
all = Table[colorCircleParts[createCircleParts[spiral /@ spoints, i, j],
cols], {i, -5, 6}, {j, 0, 12}];
Graphics[all, PlotRange -> {{-20, 20}, {-20, 20}}]
```

## Aftermath: Getting something close to the artist's work

The webpage of the artist suggests that

The central part of the picture shows the Doyle spiral circle packing P = 2 Q = 12.

That is not true. The values of P and Q define how many circles you need to close *one loop*. Additionally, the rotation of the circles in the artist's work is clockwise while in mathematics, we usually prefer to go counter-clockwise.

Lucky for us, this is no big deal because to go clockwise we just need to conjugate our complex values `a`

and `b`

. After printing the painting and counting the circles (and paying absolutely no attention to Wjx's comment who already found out that the values are off), I discovered that the painting uses P=3 and Q=8.

Let me show you what that means:

```
pqPlot[p_, q_] := Module[{a, b, r, c1, c2},
{a, b, r} = doyle[p, q];
{a, b} = Conjugate /@ {a, b};
c1 = toCirle[#, r] & /@ NestList[a*# &, 1, p - 1];
c2 = toCirle[#, r] & /@ NestList[b*# &, 1, q - 1];
Graphics[{EdgeForm[Black], FaceForm[LightYellow], c2,
FaceForm[LightBlue], c1, FaceForm[LightGreen], EdgeForm[Thick],
toCirle[1, r], toCirle[a^p, r]}]
]
pqPlot[3, 8]
```

If you include the inner base circle in your counting, you have 3 circles in the first and 8 circles in the other direction until you reach the outer end circle.

Taking this into account and including some of the colors in the original painting, we can come up with a very close optical copy of what the artist did. I played around with the plot-range to make it fit.

```
{a, b, r} = doyle[3, 8];
{a, b} = Conjugate /@ {a, b};
spoints = {a*b^-1, a, b};
cols = {GrayLevel[0.1], RGBColor[0.078, 0.71, 0.964],
RGBColor[0.95, 0.36, 0.09], RGBColor[0.77, 0.17, 0.12],
RGBColor[0.07, 0.6, 0.25], RGBColor[.32, .24, .55]};
range = 5.585;
Graphics[
Table[colorCircleParts[createCircleParts[spiral /@ spoints, i, j],
cols], {i, -5, 2}, {j, 0, 7}],
PlotRange -> {{-range, range}, {-range, range}}
]
```

6It would be good if you could add sources. A link to where you got the image, a link to support the claim that "the painting is actually related to Doyle spirals", and a link to the "Wolfram demo" that you mention. Maybe also a link to a page about Doyle spirals if you have one you think is good, because why not. And the usual question, to which you are not new: What are your thoughts on this problem? Have you tried anything? – C. E. – 2017-05-18T23:01:53.523

5A small discovery: the case in the painting corresponds to parameter

`P=4, Q=9`

in the second Wolfram demo. – Wjx – 2017-05-27T15:17:20.1171

@Wjx I really should pay attention to comments. After seeing that my spirals were off, I printed out the original and counted the circles only to discover what you already have found out. Did you see that she explicitly says that she used P=2 and Q=12 on her webpage? Very nice finding :)

– halirutan – 2017-06-02T23:12:12.670