Skip to content

Instantly share code, notes, and snippets.

@robinhouston

robinhouston/doyle.js

Last active December 23, 2015 09:16
Show Gist options
  • Star 1 You must be signed in to star a gist
  • Fork 1 You must be signed in to fork a gist
  • Save robinhouston/6096950 to your computer and use it in GitHub Desktop.
Save robinhouston/6096950 to your computer and use it in GitHub Desktop.
Download ZIP
Doyle spiral explorer
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 spiral explorer</title>
<script src="rAF.js" charset="utf-8"></script>
<script src="doyle.js" charset="utf-8"></script>
<style>
body { font-family: sans-serif; }
#container { position: relative; width: 960px; height: 500px; }
form { position: absolute; left: 804px; top: 0; padding: 8px; width: 140px; }
form div { margin-bottom: 8px; font-size: 18px; }
form input { width: 50px; font-size: 18px; }
</style>
<div id="container">
<canvas width=804 height=500></canvas>
<form>
<div>Adjust the parameters p and q to change the number of spiral arms:</div>
<div>p: <input type="number" name="p" value="8" min="2" max="23"></div>
<div>q: <input type="number" name="q" value="23" min="8"></div>
</form>
</div>
<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.i + 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;
}
}
// Interaction
var p = 8, q = 23;
function redisplay() {
var root = doyle(p, q);
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 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: 1, max_d: canvas.width * 3
});
start = cmul(start, root.b);
}
}
context.strokeStyle = "rgba(0,0,0,0)";
var params = {};
(function (query, re, match) {
while (match = re.exec(query)) {
params[decodeURIComponent(match[1])] = decodeURIComponent(match[2]);
}
})(window.location.search.substring(1).replace(/\+/g, "%20"), /([^&=]+)=?([^&]*)/g);
var form = document.getElementsByTagName("form")[0],
p_control = form.elements[0],
q_control = form.elements[1];
p_control.addEventListener("change", function(e) {
p = this.value;
q_control.setAttribute("min", Math.max(3, p));
redisplay();
}, false);
q_control.addEventListener("change", function(e) {
q = this.value;
p_control.setAttribute("max", q);
redisplay();
}, false);
if (params.p && params.q) {
p = p_control.value = parseInt(params.p);
q = q_control.value = parseInt(params.q);
p_control.setAttribute("max", q);
q_control.setAttribute("min", Math.max(3, p));
}
redisplay();
</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