# Doyle Spiral Circle Packing

Greetings,

I am trying to achieve what is known as Doyle Spiral packing (space filling)

the mathematical equation to map the center points of the circle is complex, on researching i came across the mathematicians original paper with the equation.

i am trying to recreate this on grasshopper, but i am facing some problems

there are a few more links on this topic, which actually share the source code but unfortunately its written in javascript, if someone could take a look at it, would be really helpful. i have been trying it since long.

bridges2008-131 (1).pdf (965.5 KB)

Key features to note about the doyle packing -
Doyle spiral is a pattern of non-crossing circles in the plane, each [tangent] to six others. The sequences of circles linked to each other through opposite points of tangency lie on [logarithmic spirals]

summary on the mathematicianâ€™s work (preferable to read from this point after reading the equation)
Controls/Parameters

at last, thank you to @DanielPiker , @laurent_delrieu , @seghierkhaled for all the input to develop this research

3 Likes

Re hello
with this script in Javascript you could manage to make it in C# or whatever. Yes there are some complex number but you can treat them as double tab = new double{ x,y};

Here is the script. There is a component that allows to output the spiral arms (2 or 3).
Doyle Spiral Legacy.gh (27.6 KB) *

18 Likes

Based on this

doyle_spec.gh (9.9 KB)

3 Likes

Thanks @laurent_delrieu , excellent work
Is the script calculate the rotation angle and the proportion between circles radius?
I notice that the proportion and the angle always the same and changed based on p number.

2 Likes

I canâ€™t say. I just translate from JavaScript to C#. I didnâ€™t try to have a better understanding on what is done. It worked!

2 Likes

I think the values in your script maybe two of them are the angle and the proportion (serie step) but i donâ€™t know how to setup the outputs of all values to find these two numbers

I try to convert js file to python without success

Javascript

``````doyle.js#
/*   Numerics for Doyle spirals.
*   Robin Houston, 2013
*/

(function() {
var pow = Math.pow,
sin = Math.sin,
cos = Math.cos,
pi = Math.PI;

function _d(z,t, p,q) {
// The square of the distance between z*e^(it) and z*e^(it)^(p/q).
var w = pow(z, p/q),
s =(p*t + 2*pi)/q;
return (
pow( z*cos(t) - w*cos(s), 2 )
+ pow( z*sin(t) - w*sin(s), 2 )
);
}

function ddz_d(z,t, p,q) {
// The partial derivative of _d with respect to z.
var w = pow(z, p/q),
s = (p*t + 2*pi)/q,
ddz_w = (p/q)*pow(z, (p-q)/q);
return (
2*(w*cos(s) - z*cos(t))*(ddz_w*cos(s) - cos(t))
+ 2*(w*sin(s) - z*sin(t))*(ddz_w*sin(s) - sin(t))
);
}

function ddt_d(z,t, p,q) {
// The partial derivative of _d with respect to t.
var w = pow(z, p/q),
s = (p*t + 2*pi)/q,
dds_t = (p/q);
return (
2*( z*cos(t) - w*cos(s) )*( -z*sin(t) + w*sin(s)*dds_t )
+ 2*( z*sin(t) - w*sin(s) )*( z*cos(t) - w*cos(s)*dds_t )
);
}

function _s(z,t, p,q) {
// The square of the sum of the origin-distance of z*e^(it) and
// the origin-distance of z*e^(it)^(p/q).
return pow(z + pow(z, p/q), 2);
}

function ddz_s(z,t, p,q) {
// The partial derivative of _s with respect to z.
var w = pow(z, p/q),
ddz_w = (p/q)*pow(z, (p-q)/q);
return 2*(w+z)*(ddz_w+1);
}

/*
function ddt_s(z,t, p,q) {
// The partial derivative of _s with respect to t.
return 0;
}
*/

function _r(z,t, p,q) {
// The square of the radius-ratio implied by having touching circles
// centred at z*e^(it) and z*e^(it)^(p/q).
return _d(z,t,p,q) / _s(z,t,p,q);
}

function ddz_r(z,t, p,q) {
// The partial derivative of _r with respect to z.
return (
ddz_d(z,t,p,q) * _s(z,t,p,q)
- _d(z,t,p,q) * ddz_s(z,t,p,q)
) / pow( _s(z,t,p,q), 2 );
}

function ddt_r(z,t, p,q) {
// The partial derivative of _r with respect to t.
return (
ddt_d(z,t,p,q) * _s(z,t,p,q)
/* - _d(z,t,p,q) * ddt_s(z,t,p,q) */  // omitted because ddt_s is constant at zero
) / pow( _s(z,t,p,q), 2 );
}

var epsilon = 1e-10;
window.doyle = function(p, q) {
// We want to find (z, t) such that:
//    _r(z,t,0,1) = _r(z,t,p,q) = _r(pow(z, p/q), (p*t + 2*pi)/q, 0,1)
//
// so we define functions _f and _g to be zero when these equalities hold,
// and use 2d Newton-Raphson to find a joint root of _f and _g.

function _f(z, t) {
return _r(z,t,0,1) - _r(z,t,p,q);
}

function ddz_f(z, t) {
return ddz_r(z,t,0,1) - ddz_r(z,t,p,q);
}

function ddt_f(z, t) {
return ddt_r(z,t,0,1) - ddt_r(z,t,p,q);
}

function _g(z, t) {
return _r(z,t,0,1) - _r(pow(z, p/q), (p*t + 2*pi)/q, 0,1);
}

function ddz_g(z, t) {
return ddz_r(z,t,0,1) - ddz_r(pow(z, p/q), (p*t + 2*pi)/q, 0,1) * (p/q)*pow(z, (p-q)/q);
}

function ddt_g(z, t) {
return ddt_r(z,t,0,1) - ddt_r(pow(z, p/q), (p*t + 2*pi)/q, 0,1) * (p/q);
}

function find_root(z, t) {
for(;;) {
var v_f = _f(z, t),
v_g = _g(z, t);
if (-epsilon < v_f && v_f < epsilon && -epsilon < v_g && v_g < epsilon)
return {ok: true, z: z, t: t, r: Math.sqrt(_r(z,t,0,1))};

var a = ddz_f(z,t), b = ddt_f(z,t), c = ddz_g(z,t), d = ddt_g(z,t);
var det = a*d-b*c;
if (-epsilon < det && det < epsilon)
return {ok: false};

z -= (d*v_f - b*v_g)/det;
t -= (a*v_g - c*v_f)/det;

if (z < epsilon)
return {ok: false};
}
}

var root = find_root(2, 0);
if (!root.ok) throw "Failed to find root for p=" + p + ", q=" + q;

var a = [root.z * cos(root.t), root.z * sin(root.t) ],
coroot = {z: pow(root.z, p/q), t: (p*root.t+2*pi)/q},
b = [coroot.z * cos(coroot.t), coroot.z * sin(coroot.t) ];
return {a: a, b: b, r: root.r, mod_a: root.z, arg_a: root.t};
};
})();
``````

Python

``````import math

pow = math.pow
sin = math.sin
cos = math.cos
pi = math.pi

def _d(z,t,p,q):
w = pow(z,p/q)
s = (p*t+2*pi)/q
v = pow(z*cos(t)-w*cos(s),2)+pow(z*sin(t)-w*sin(s),2)
return v

def ddz_d(z,t,p,q):
w = pow(z,p/q)
s = (p*t+2*pi)/q
ddz_w = (p/q)*pow(z,(p-q)/q)
v = 2*(w*cos(s) - z*cos(t))*(ddz_w*cos(s) - cos(t))+ 2*(w*sin(s) - z*sin(t))*(ddz_w*sin(s) - sin(t))
return v

def ddt_d(z,t, p,q):
w = pow(z, p/q)
s = (p*t + 2*pi)/q
dds_t = (p/q)
v = 2*( z*cos(t) - w*cos(s) )*( -z*sin(t) + w*sin(s)*dds_t )+ 2*( z*sin(t) - w*sin(s) )*( z*cos(t) - w*cos(s)*dds_t )
return v

def _s(z,t, p,q):
v = pow(z + pow(z, p/q), 2)
return v

def ddz_s(z,t, p,q):
w = pow(z, p/q)
ddz_w = (p/q)*pow(z, (p-q)/q)
v = 2*(w+z)*(ddz_w+1)
return v

def _r(z,t, p,q):
v = _d(z,t,p,q) / _s(z,t,p,q)
return v

def ddz_r(z,t, p,q):
v = (ddz_d(z,t,p,q) * _s(z,t,p,q) - _d(z,t,p,q) * ddz_s(z,t,p,q)) / pow( _s(z,t,p,q), 2 )
return v

def ddt_r(z,t, p,q):
v = (ddt_d(z,t,p,q) * _s(z,t,p,q)) / pow( _s(z,t,p,q), 2 )
return v

epsilon = 1e-10

def _f(z, t):
return _r(z,t,0,1) - _r(z,t,p,q)

def ddz_f(z, t):
return ddz_r(z,t,0,1) - ddz_r(z,t,p,q)

def ddt_f(z, t):
return ddt_r(z,t,0,1) - ddt_r(z,t,p,q)

def _g(z, t):
return _r(z,t,0,1) - _r(pow(z, p/q), (p*t + 2*pi)/q, 0,1)

def ddz_g(z, t):
return ddz_r(z,t,0,1) - ddz_r(pow(z, p/q), (p*t + 2*pi)/q, 0,1) * (p/q)*pow(z, (p-q)/q)

def ddt_g(z, t):
return ddt_r(z,t,0,1) - ddt_r(pow(z, p/q), (p*t + 2*pi)/q, 0,1) * (p/q)

def find_root(z, t):
ok = False
z = 0.0
t = 0.0
r = 0.0
for i in range(100):
v_f = _f(z, t)
v_g = _g(z, t)
if -epsilon < v_f and v_f < epsilon and -epsilon < v_g and v_g < epsilon:
ok = True; z= z; t= t; r= Math.sqrt(_r(z,t,0,1))
return ok,z,t,r
break

a = ddz_f(z,t)
b = ddt_f(z,t)
c = ddz_g(z,t)
d = ddt_g(z,t)
det = a*d-b*c

if -epsilon < det and det < epsilon:
ok = False
return ok
break

z -= (d*v_f - b*v_g)/det
t -= (a*v_g - c*v_f)/det

if z < epsilon:
ok = False
return ok
break

root = find_root(2,0)

if not root.ok:
a = [root.z * cos(root.t), root.z * sin(root.t) ]
coroot.z = pow(root.z, p/q)
coroot.t = (p*root.t+2*pi)/q
coroot.ok = True
coroot.r = 0.0
b = [coroot.z * cos(coroot.t), coroot.z * sin(coroot.t) ]
a= a; b= b; r=root.r; mod_a= root.z; arg_a= root.t
``````

doyle_python_try.gh (3.5 KB)

@seghierkhaled, you can enclose your code blocks in triple backticks and declare the language right after the first triplet in markdown. Makes things much more legible!

``````if __name__ == "__main__":
print "Hello World!"
``````

``````include <iostream>

int main() {
std::cout << "Hello World!\n";
}
``````
1 Like

Thanks but i donâ€™t understand what you mean

this is the try to create circles, i donâ€™t know where is the error

doyle_python_try.gh (8.9 KB)

1 Like

With python and a little changes in the code

doyle_pyt.gh (6.8 KB)

10 Likes

awesome. good thing you kept at it

1 Like

I try to create Double Doyle Spiral but the code miss something

doyle_python.gh (19.1 KB)

3 Likes

You could get the loxodrome arrangement by using a MĂ¶bius transformation on the output

doubledoyle.gh (12.5 KB)

17 Likes
10 Likes

This is one of the cases where i would like the thread to be tagged with the â€śGalleryâ€ť tabâ€¦

Wonderful stuff guys!

3 Likes

Whatâ€™s the second to last component please? The one with the rainbow colour circle and inputs C & S, output C.

Can this pattern be applied to a curved mesh?

This component from Fennec addon, you can create one with python.
I think you can apply the pattern to a curved mesh with transformation components

11 Likes
10 Likes

OMG, it looks like a tie dye t-shirt now.

1 Like