Doyle spiral circle packing

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

index.html

<!DOCTYPE html>
<title>Doyle spirals</title>
<script src="rAF.js" charset="utf-8"></script>
<script src="doyle.js" charset="utf-8"></script>
<canvas width=960 height=500></canvas>
<script>
    // Initialisation
    var canvas = document.getElementsByTagName("canvas")[0],
        context = canvas.getContext("2d");
    
    
    // Circle drawing
    function circle(x, y, r) {
        context.beginPath();
        context.arc(x, y, r, 0, 7, false);
        context.fill();
        context.lineWidth = r/10;
        context.stroke();
    }
    
    
    // Complex arithmetic
    function cmul(w, z) {
        return [
            w[0]*z[0] - w[1]*z[1],
            w[0]*z[1] + w[1]*z[0]
        ];
    }
    function modulus(p) {
        return Math.sqrt(p[0]*p[0] + p[1]*p[1]);
    }
    function crecip(z) {
        var d = z[0]*z[0] + z[1]*z[1];
        return [z[0]/d, -z[1]/d];
    }
    
    
    // Doyle spiral drawing
    function spiral(r, start_point, delta, options) {
        var recip_delta = crecip(delta),
            mod_delta = modulus(delta),
            mod_recip_delta = 1/mod_delta,
            color_index = options.i,
            colors = options.fill,
            min_d = options.min_d,
            max_d = options.max_d;
        
        // Spiral outwards
        for (var q = start_point, mod_q = modulus(q);
            mod_q < max_d;
            q = cmul(q, delta), mod_q *= mod_delta
        ) {
            context.fillStyle = colors[color_index];
            circle(q[0], q[1], mod_q*r);
            color_index = (color_index + 1) % colors.length;
        }
        
        // Spiral inwards
        color_index = ((options ? options.i : 0) + colors.length - 1) % colors.length;
        for (var q = cmul(start_point, recip_delta), mod_q = modulus(q);
            mod_q > min_d;
            q = cmul(q, recip_delta), mod_q *= mod_recip_delta
        ) {
            context.fillStyle = colors[color_index];
            circle(q[0], q[1], mod_q*r);
            color_index = (color_index + colors.length - 1) % colors.length;
        }
    }
    
    
    // Animation
    var p = 9, q = 24;
    var root = doyle(p, q);
    
    var ms_per_repeat = 300;
    
    function frame(t) {
        context.setTransform(1, 0, 0, 1, 0, 0);
        context.clearRect(0, 0, canvas.width, canvas.height);
        
        context.translate(Math.round(canvas.width/2), cy = Math.round(canvas.height/2));
        var scale = Math.pow(root.mod_a, t);
        context.scale(scale, scale);
        context.rotate(root.arg_a * t);
        
        var min_d = 1/scale, max_d = canvas.width * 2;
        var start = root.a;
        for (var i=0; i<q; i++) {
            spiral(root.r, start, root.a, {
                fill: ["#49B", "#483352", "#486078"], i: (2*i)%3,
                min_d: min_d, max_d: max_d
            });
            start = cmul(start, root.b);
        }
    }
    
    var first_timestamp;
    function loop(timestamp) {
        if (!first_timestamp) first_timestamp = timestamp;
        frame(((timestamp - first_timestamp) % (ms_per_repeat*3)) / ms_per_repeat);
        requestAnimationFrame(loop);
    }
    
    context.strokeStyle = "rgba(0,0,0,0)"; //"white"; //"hsl(228, 10%, 35%)";
    requestAnimationFrame(loop);
</script>

rAF.js

// http://paulirish.com/2011/requestanimationframe-for-smart-animating/
// http://my.opera.com/emoller/blog/2011/12/20/requestanimationframe-for-smart-er-animating

// requestAnimationFrame polyfill by Erik Möller. fixes from Paul Irish and Tino Zijdel

// MIT license

(function() {
    var lastTime = 0;
    var vendors = ['ms', 'moz', 'webkit', 'o'];
    for(var x = 0; x < vendors.length && !window.requestAnimationFrame; ++x) {
        window.requestAnimationFrame = window[vendors[x]+'RequestAnimationFrame'];
        window.cancelAnimationFrame = window[vendors[x]+'CancelAnimationFrame'] 
                                   || window[vendors[x]+'CancelRequestAnimationFrame'];
    }
 
    if (!window.requestAnimationFrame)
        window.requestAnimationFrame = function(callback, element) {
            var currTime = new Date().getTime();
            var timeToCall = Math.max(0, 16 - (currTime - lastTime));
            var id = window.setTimeout(function() { callback(currTime + timeToCall); }, 
              timeToCall);
            lastTime = currTime + timeToCall;
            return id;
        };
 
    if (!window.cancelAnimationFrame)
        window.cancelAnimationFrame = function(id) {
            clearTimeout(id);
        };
}());

Thanks for the code. I ported your Doyle Spiral code to PovRay and made modifications (of course) so it displays a static picture. I checked the project into GitHub.