Skip to content

Instantly share code, notes, and snippets.

@davo
Forked from robinhouston/index.html
Last active December 27, 2023 01:21
Show Gist options
  • Save davo/f78fe7d576d660daf70548942b71bd05 to your computer and use it in GitHub Desktop.
Save davo/f78fe7d576d660daf70548942b71bd05 to your computer and use it in GitHub Desktop.
Stumbling blocks
/* 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};
};
})();
<!DOCTYPE html>
<title>Stumbling blocks</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");
// 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, gamma, 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,
delta_gamma = cmul(delta, gamma),
delta_over_gamma = cmul(delta, crecip(gamma)),
gamma_over_delta = cmul(gamma, crecip(delta));
// Spiral inwards
for (var q = start_point, mod_q = modulus(q);
mod_q > min_d;
q = cmul(q, recip_delta), mod_q *= mod_recip_delta
) {
color_index = (color_index + colors.length - 1) % colors.length;
}
// Spiral outwards
for (;mod_q < max_d; q = cmul(q, delta), mod_q *= mod_delta) {
context.fillStyle = colors[color_index];
context.beginPath();
context.moveTo.apply(context, q);
if (color_index == 0)
context.lineTo.apply(context, cmul(q, delta_over_gamma));
context.lineTo.apply(context, cmul(q, delta));
if (color_index == 2)
context.lineTo.apply(context, cmul(q, delta_gamma));
context.lineTo.apply(context, cmul(q, gamma));
if (color_index == 1)
context.lineTo.apply(context, cmul(q, gamma_over_delta));
context.closePath();
context.fill();
color_index = (color_index + 1) % colors.length;
}
}
// Animation
var p = 9, q = 24;
var root = doyle(p, q);
var ms_per_repeat = 1000;
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, -3*t);
context.scale(scale, scale);
context.rotate(-3*root.arg_a * t);
var min_d = 1/scale, max_d = canvas.width * 0.7 / scale;
var start = root.a;
for (var i=0; i<q; i++) {
spiral(root.r, start, root.a, root.b, {
fill: ["#D2D2D2", "#676767", "#050505"], 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);
}
requestAnimationFrame(loop);
</script>
// 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