doyle = function (p, q) {
var epsilon = 1e-10;
var pow = Math.pow,
sin = Math.sin,
cos = Math.cos,
pi = Math.PI;
function _d(z, t, 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) {
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) {
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) {
return pow(z + pow(z, p / q), 2);
}
function ddz_s(z, t, p, q) {
var w = pow(z, p / q),
ddz_w = (p / q) * pow(z, (p - q) / q);
return 2 * (w + z) * (ddz_w + 1);
}
function _r(z, t, p, q) {
return _d(z, t, p, q) / _s(z, t, p, q);
}
function ddz_r(z, t, p, q) {
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) {
return (
(ddt_d(z, t, p, q) * _s(z, t, p, q)) /
pow(_s(z, t, p, q), 2)
);
}
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)];
console.log(a, b);
return { a: a, b: b, r: root.r, mod_a: root.z, arg_a: root.t };
}