davo · December 27, 2023 01:18 · Dec 27, 2023 · Dec 27, 2023 · Jul 28, 2013 · Jul 27, 2013
diff --git a/index.html b/index.html
@@ -1,50 +1,21 @@
 <!DOCTYPE html>
-<meta charset="utf-8">
-<title>Doyle spiral + Möbius transformation = &#x1F632;</title>
-<script src="http://bl.ocks.org/robinhouston/raw/6096562/rAF.js" charset="utf-8"></script>
-<script src="http://bl.ocks.org/robinhouston/raw/6096562/doyle.js" charset="utf-8"></script>
+<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 and transformation
 
-    // Möbius transformation that maps (0, 1, ∞) to (-1, 0, 1)
-    function transform_point(x, y) {
-        var denom = (x+1)*(x+1) + y*y;
-        return [ (x*x - 1 + y*y)/denom, 2*y/denom ];
-    }
-    // The image of a circle under the Möbius transformation
-    function transform_circle(x, y, r) {
-        var a = transform_point(x-r, y),
-            b = transform_point(x+r, y),
-            c = transform_point(x, y+r);
-        return circle_through_points(a,b,c);
-    }
-    // The unique circle passing through three non-collinear points
-    function circle_through_points(a, b, c) {
-        var na = a[0]*a[0] + a[1]*a[1],
-            nb = b[0]*b[0] + b[1]*b[1],
-            nc = c[0]*c[0] + c[1]*c[1];
-        var y = (
-            (a[0]-b[0])*(nb-nc) - (b[0]-c[0])*(na-nb)
-        ) / (
-            2*(b[1]-a[1])*(b[0]-c[0]) - 2*(a[0]-b[0])*(c[1]-b[1])
-        ),
-            x = (na-nb + 2*(b[1]-a[1])*y) / (2*(a[0]-b[0])),
-            r = Math.sqrt( (x-a[0])*(x-a[0]) + (y-a[1])*(y-a[1]) );
-        return [x, y, r];
-    }
-    // Draw a circle!
+    // Circle drawing
     function circle(x, y, r) {
-        var c = transform_circle(x, y, r);
-        if (c[2] > 10) return;
         context.beginPath();
-        context.arc(c[0], c[1], c[2], 0, 7, false);
+        context.arc(x, y, r, 0, 7, false);
         context.fill();
+        context.lineWidth = r/10;
+        context.stroke();
     }
 
 
@@ -55,9 +26,6 @@
             w[0]*z[1] + w[1]*z[0]
         ];
     }
-    function rotate(z, theta) {
-        return cmul(z, [Math.cos(theta), Math.sin(theta)]);
-    }
     function modulus(p) {
         return Math.sqrt(p[0]*p[0] + p[1]*p[1]);
     }
@@ -104,20 +72,23 @@
     var p = 9, q = 24;
     var root = doyle(p, q);
 
-    var ms_per_repeat = 10000;
+    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));
-        context.scale(200, 200);
+        var scale = Math.pow(root.mod_a, t);
+        context.scale(scale, scale);
+        context.rotate(root.arg_a * t);
 
-        var start = rotate(root.a, Math.PI*2*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: 1e-3, max_d: 1e3
+                min_d: min_d, max_d: max_d
             });
             start = cmul(start, root.b);
         }
@@ -126,9 +97,10 @@
     var first_timestamp;
     function loop(timestamp) {
         if (!first_timestamp) first_timestamp = timestamp;
-        frame(((timestamp - first_timestamp) % ms_per_repeat) / ms_per_repeat);
+        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>
diff --git a/index.html b/index.html
@@ -1,21 +1,50 @@
 <!DOCTYPE html>
-<title>Doyle spirals</title>
-<script src="rAF.js" charset="utf-8"></script>
-<script src="doyle.js" charset="utf-8"></script>
+<meta charset="utf-8">
+<title>Doyle spiral + Möbius transformation = &#x1F632;</title>
+<script src="http://bl.ocks.org/robinhouston/raw/6096562/rAF.js" charset="utf-8"></script>
+<script src="http://bl.ocks.org/robinhouston/raw/6096562/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 and transformation
 
-    // Circle drawing
+    // Möbius transformation that maps (0, 1, ∞) to (-1, 0, 1)
+    function transform_point(x, y) {
+        var denom = (x+1)*(x+1) + y*y;
+        return [ (x*x - 1 + y*y)/denom, 2*y/denom ];
+    }
+    // The image of a circle under the Möbius transformation
+    function transform_circle(x, y, r) {
+        var a = transform_point(x-r, y),
+            b = transform_point(x+r, y),
+            c = transform_point(x, y+r);
+        return circle_through_points(a,b,c);
+    }
+    // The unique circle passing through three non-collinear points
+    function circle_through_points(a, b, c) {
+        var na = a[0]*a[0] + a[1]*a[1],
+            nb = b[0]*b[0] + b[1]*b[1],
+            nc = c[0]*c[0] + c[1]*c[1];
+        var y = (
+            (a[0]-b[0])*(nb-nc) - (b[0]-c[0])*(na-nb)
+        ) / (
+            2*(b[1]-a[1])*(b[0]-c[0]) - 2*(a[0]-b[0])*(c[1]-b[1])
+        ),
+            x = (na-nb + 2*(b[1]-a[1])*y) / (2*(a[0]-b[0])),
+            r = Math.sqrt( (x-a[0])*(x-a[0]) + (y-a[1])*(y-a[1]) );
+        return [x, y, r];
+    }
+    // Draw a circle!
     function circle(x, y, r) {
+        var c = transform_circle(x, y, r);
+        if (c[2] > 10) return;
         context.beginPath();
-        context.arc(x, y, r, 0, 7, false);
+        context.arc(c[0], c[1], c[2], 0, 7, false);
         context.fill();
-        context.lineWidth = r/10;
-        context.stroke();
     }
 
 
@@ -26,6 +55,9 @@
             w[0]*z[1] + w[1]*z[0]
         ];
     }
+    function rotate(z, theta) {
+        return cmul(z, [Math.cos(theta), Math.sin(theta)]);
+    }
     function modulus(p) {
         return Math.sqrt(p[0]*p[0] + p[1]*p[1]);
     }
@@ -72,23 +104,20 @@
     var p = 9, q = 24;
     var root = doyle(p, q);
 
-    var ms_per_repeat = 300;
+    var ms_per_repeat = 10000;
 
     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);
+        context.scale(200, 200);
 
-        var min_d = 1/scale, max_d = canvas.width * 2;
-        var start = root.a;
+        var start = rotate(root.a, Math.PI*2*t);
         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
+                min_d: 1e-3, max_d: 1e3
             });
             start = cmul(start, root.b);
         }
@@ -97,10 +126,9 @@
     var first_timestamp;
     function loop(timestamp) {
         if (!first_timestamp) first_timestamp = timestamp;
-        frame(((timestamp - first_timestamp) % (ms_per_repeat*3)) / ms_per_repeat);
+        frame(((timestamp - first_timestamp) % ms_per_repeat) / ms_per_repeat);
         requestAnimationFrame(loop);
     }
 
-    context.strokeStyle = "rgba(0,0,0,0)"; //"white"; //"hsl(228, 10%, 35%)";
     requestAnimationFrame(loop);
-</script>
+</script>
diff --git a/index.html b/index.html
@@ -69,7 +69,7 @@
 
 
     // Animation
-    var p = 8, q = 23;
+    var p = 9, q = 24;
     var root = doyle(p, q);
 
     var ms_per_repeat = 300;
@@ -83,11 +83,11 @@
         context.scale(scale, scale);
         context.rotate(root.arg_a * t);
 
-        var min_d = 1/scale, max_d = canvas.width * 0.71 / scale;
+        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,
+                fill: ["#49B", "#483352", "#486078"], i: (2*i)%3,
                 min_d: min_d, max_d: max_d
             });
             start = cmul(start, root.b);

diff --git a/index.html b/index.html
@@ -75,7 +75,6 @@
     var ms_per_repeat = 300;
 
     function frame(t) {
-        context.save();
         context.setTransform(1, 0, 0, 1, 0, 0);
         context.clearRect(0, 0, canvas.width, canvas.height);
 
@@ -93,7 +92,6 @@
             });
             start = cmul(start, root.b);
         }
-        context.restore();
     }
 
     var first_timestamp;
@@ -103,6 +101,6 @@
         requestAnimationFrame(loop);
     }
 
-    context.strokeStyle = "none"; //"white"; //"hsl(228, 10%, 35%)";
+    context.strokeStyle = "rgba(0,0,0,0)"; //"white"; //"hsl(228, 10%, 35%)";
     requestAnimationFrame(loop);
 </script>
diff --git a/doyle.js b/doyle.js
@@ -0,0 +1,146 @@
+/*   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};
+    };
+})();
diff --git a/rAF.js b/rAF.js
@@ -0,0 +1,31 @@
+// 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);
+        };
+}());
diff --git a/index.html b/index.html
@@ -0,0 +1,108 @@
+<!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 = 8, q = 23;
+    var root = doyle(p, q);
+
+    var ms_per_repeat = 300;
+
+    function frame(t) {
+        context.save();
+        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 * 0.71 / scale;
+        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);
+        }
+        context.restore();
+    }
+
+    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 = "none"; //"white"; //"hsl(228, 10%, 35%)";
+    requestAnimationFrame(loop);
+</script>