-
-
Save davo/22f15ef4b034d4bee5655508611185a9 to your computer and use it in GitHub Desktop.
Doyle spiral circle packing
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
/* 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}; | |
}; | |
})(); |
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
<!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> |
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
// 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