const

{atan2, cos, exp, log, sin} = Math,

rotate1 = (array) => [...array.slice(1), array[0]], // returns [A[1], A[2], ... A[n-1], A[0]]

even = (array) => array.filter((x, i) => !(i % 2)),

odd = (array) => array.filter((x, i) => i % 2),

evenodd = (array) => [even(array), odd(array)],

displacement = (array0, array1) => array0.map((x, i) => array1[i] - x),

data = scpie_diskmap_data,

center = data.center,

n = data.vertices.length >> 1,

[nodes0, nodes1] = evenodd(data.nodes),

[values0, values1] = evenodd(data.values),

[weights0, weights1] = evenodd(data.weights),

[vertices0x, vertices0y] = evenodd(data.vertices),

vertices1x = rotate1(vertices0x),

vertices1y = rotate1(vertices0y),

sidesx = displacement(vertices0x, vertices1x),

sidesy = displacement(vertices0y, vertices1y),

[prevertices0x, prevertices0y] = evenodd(data.prevertices),

prevertices1x = rotate1(prevertices0x),

prevertices1y = rotate1(prevertices0y),

[piewidthsx, piewidthsy] = evenodd(data.piewidths),

exponents0 = data.exponents,

exponents1 = rotate1(exponents0);

const scpie_eval = function scpie_eval (x, y) {

if (!isFinite(x + y)) return NaN;

// figure out which sector the point is in

let k;

for (k = 0;; k++) {

if (k === n) return center.slice();

const z_left_z0 = (prevertices0x[k] * y - prevertices0y[k] * x) >= 0;

const z_right_z1 = (prevertices1x[k] * y - prevertices1y[k] * x) < 0;

if (z_left_z0 & z_right_z1) break;

// special case when the sector between k & k+1 is wider than half the disk:

if (piewidthsx[k] <= 0) { if (z_left_z0 | z_right_z1) break; }

}

const

d0x = x - prevertices0x[k],

d0y = y - prevertices0y[k],

d1x = x - prevertices1x[k],

d1y = y - prevertices1y[k],

hx = piewidthsx[k],

hy = piewidthsy[k];

// half-sector closest to prevertex k

if ((d0x*d0x + d0y*d0y) <= (d1x*d1x + d1y*d1y)) {

const // q = denominator of moebius map

qx = d0x - (hx * d1x + hy * d1y),

qy = d0y - (hx * d1y - hy * d1x);

x = (d0x * qx + d0y * qy) * (1 / (qx * qx + qy * qy));

y = (d0y * qx - d0x * qy) * (1 / (qx * qx + qy * qy));

const // complex power

r = exp(0.5 * log(x*x + y*y) * exponents0[k]),

a = atan2(y, x) * exponents0[k];

x = r * cos(a);

y = r * sin(a);

// rational approximation

let out = reval(nodes0[k], values0[k], weights0[k], x, y);

x = out[0]; y = out[1];

// transform to final position; reuse out array

out[0] = vertices0x[k] + (sidesx[k] * x - sidesy[k] * y);

out[1] = vertices0y[k] + (sidesy[k] * x + sidesx[k] * y);

return out;

}

// half-sector closest to prevertex k+1

else {

const // q = denominator of moebius map

qx = d1x - (hx * d0x - hy * d0y),

qy = d1y - (hx * d0y + hy * d0x);

x = (d1x * qx + d1y * qy) * (1 / (qx * qx + qy * qy));

y = (d1y * qx - d1x * qy) * (1 / (qx * qx + qy * qy));

const // complex power

r = exp(0.5 * log(x*x + y*y) * exponents1[k]),

a = atan2(y, x) * exponents1[k];

x = r * cos(a);

y = r * sin(a);

// rational approximation

let out = reval(nodes1[k], values1[k], weights1[k], x, y);

x = out[0]; y = out[1];

// transform to final position; reuse out array

out[0] = vertices1x[k] - (sidesx[k] * x - sidesy[k] * y);

out[1] = vertices1y[k] - (sidesy[k] * x + sidesx[k] * y);

return out;

}

}

Object.assign(scpie_eval, scpie_diskmap_data);

return scpie_eval;

}

const n = nodes.length;

// running total for numerator & denominator

let px = 0, py = 0, qx = 0, qy = 0;

for (let j = 0; j < n; j += 2) {

// q += w_j / (z - z_j)

// p += f_j * w_j / (z - z_j)

const

xj = x - nodes[j], yj = nodes[j+1] - y,

wxj = weights[j], wyj = weights[j+1],

fxj = values[j], fyj = values[j+1],

dj = 1 / (xj*xj + yj*yj),

qxj = dj * (wxj*xj - wyj*yj),

qyj = dj * (wxj*yj + wyj*xj);

px += fxj * qxj - fyj * qyj;

py += fxj * qyj + fyj * qxj;

qx += qxj, qy += qyj;

}

let

d = 1 / (qx * qx + qy * qy),

fx = (px * qx + py * qy) * d,

fy = (py * qx - px * qy) * d;

// Edge case where x + iy is one of the nodes; directly use value

if (fx + fy !== fx + fy) { // true if x or y is NaN

for (let j = 0; j < n; j += 2) {

const

xj = x - nodes[j],

yj = nodes[j+1] - y;

if ((xj === 0) & (yj === 0)) {

fx = values[j]; fy = values[j+1];

break;

}

}

}

return [fx, fy];

const p1 = 0.7548776662466927; // inverse of plastic number

const p2 = 0.5698402909980532;

return (n) => {

let squarepts = [...Array(2*n).keys()].map(i => [

2 * ((p1 * i) % 1) - 1,

2 * ((p2 * i) % 1) - 1]);

return squarepts.slice(1)

.filter(([x,y]) => x*x + y*y <= 1)

.slice(0, n);

}

JCh_in_gamut = (jch) => {

const [J, C, h] = jch;

let cm = C / 2;

return xyz_to_srgb(JzCzhz_inverse(jch).map(X => .5*X))

// TODO: implement this