Doyle Spiral Circle Packing

Grasshopper
1

Greetings,

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

ref
ref828Ă—795 74.3 KB

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.

pd002
pd002685Ă—837 60.3 KB
pd 001
pd 001692Ă—586 39.5 KB

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.

Doyle Spiral

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
2

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};

image
image1073Ă—727 118 KB

image
image911Ă—728 238 KB

image
image902Ă—648 151 KB

image
image911Ă—728 167 KB

image
image1118Ă—893 187 KB

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) *

14 Likes
3

Based on this

image
image1641Ă—375 73.1 KB

image
image1629Ă—628 58.8 KB

image
image1604Ă—576 66.1 KB

doyle_spec.gh (9.9 KB)

3 Likes
4

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.

image
image1509Ă—662 91.1 KB

image
image1684Ă—705 94.2 KB

2 Likes
5

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
6

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

7

I try to convert js file to python without success

image

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)

8

@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!

Screenshot 2021-05-02 at 08.16.10

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

Screenshot 2021-05-02 at 08.16.35

include <iostream>

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

Thanks but i don’t understand what you mean

this is the try to create circles, i don’t know where is the error

image
image1299Ă—690 43.6 KB

doyle_python_try.gh (8.9 KB)

1 Like
10

With python and a little changes in the code

doyle_pyt.gh (6.8 KB)

image
image1087Ă—661 131 KB

10 Likes
11

awesome. good thing you kept at it

1 Like
12

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

image
image1596Ă—597 58.1 KB

image
image960Ă—500 123 KB

doyle_python.gh (19.1 KB)

3 Likes
13

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

image
image1038Ă—575 74.6 KB

doubledoyle.gh (12.5 KB)
loxodrome_000

15 Likes
14

Thank you , very nice

image
image1677Ă—716 167 KB

9 Likes
15

This is one of the cases where i would like the thread to be tagged with the “Gallery” tab…

Wonderful stuff guys!

3 Likes
16

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?

17

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

18

doyle_spiral.gh (11.5 KB)

10 Likes
19

image
image808Ă—694 38.7 KB

image
image1000Ă—705 55.2 KB

9 Likes
20

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

1 Like