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.

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

https://demonstrations.wolfram.com/DoyleSpiralsAndMoebiusTransformations/

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

10 Likes

Based on this

doyle_spec.gh (9.9 KB)

2 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.

1 Like

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

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)

@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

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)

With python and a little changes in the code

doyle_pyt.gh (6.8 KB)

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

2 Likes

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


doubledoyle.gh (12.5 KB)
loxodrome_000

12 Likes

Thank you , very nice

6 Likes

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

Wonderful stuff guys!

2 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

doyle_spiral.gh (11.5 KB)

8 Likes

7 Likes

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

1 Like