Skip to content

Instantly share code, notes, and snippets.

@robinhouston

robinhouston/doyle.js

Last active December 27, 2023 01:09
Show Gist options
  • Star 5 You must be signed in to star a gist
  • Fork 2 You must be signed in to fork a gist
  • Save robinhouston/6096562 to your computer and use it in GitHub Desktop.
Save robinhouston/6096562 to your computer and use it in GitHub Desktop.
Download ZIP
Doyle spiral circle packing
Raw
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};
};
})();
Raw
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>
Raw
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);
};
}());
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment