function hexadecant(lambda, phi) {
if (phi === Math.PI / 2) return [0, 0];
if (lambda > Math.PI / 4 - 1e-12) lambda = Math.PI / 4 - 1e-12;
const phi0 = (3 * Math.PI) / 8;
const cos_phi0 = Math.cos(phi0);
const sin_phi0 = Math.sin(phi0);
const hprime =
(12 / Math.PI) *
(Math.asin(1 / Math.sqrt(2 - cos_phi0 * cos_phi0)) +
Math.asin((2 * sin_phi0) / Math.sqrt(3 - Math.cos(2 * phi0))) -
Math.PI / 2);
const xiprime = Math.atan(
(Math.PI * (hprime - 3) * (hprime - 3)) /
(Math.sqrt(3) *
(Math.PI * (hprime * hprime - 2 * hprime + 45) -
96 * Math.asin(1 / Math.sqrt(2 - cos_phi0 * cos_phi0)) -
48 * Math.asin((2 * sin_phi0) / Math.sqrt(3 - Math.cos(2 * phi0)))))
);
const lambda0 =
(Math.floor((lambda + Math.PI / 4) / (Math.PI / 2)) * Math.PI) / 2;
const theta = Math.atan2(
Math.cos(phi) * Math.sin(lambda - lambda0),
sin_phi0 * Math.cos(phi) * Math.cos(lambda - lambda0) -
cos_phi0 * Math.sin(phi)
);
const r = Math.acos(
sin_phi0 * Math.sin(phi) +
cos_phi0 * Math.cos(phi) * Math.cos(lambda - lambda0)
);
const psi0 = Math.asin(1 / Math.sqrt(2 - cos_phi0 * cos_phi0));
const psi1 = Math.PI - 2 * psi0;
let beta;
if (theta <= psi0) beta = psi0 - theta;
else if (theta <= psi0 + psi1) beta = theta - psi0;
else beta = Math.PI - theta;
const c =
theta <= psi0 + psi1
? Math.acos(cos_phi0 / Math.sqrt(2))
: Math.PI / 2 - phi0;
const G = Math.abs(Math.PI / 2 - (theta % Math.PI)) < psi1 / 2 ? psi1 : psi0;
let Gprime;
let F;
if (theta <= psi0) {
Gprime = Math.PI / 2 - Math.atan(hprime / Math.sqrt(3));
F = Math.asin((2 * sin_phi0) / Math.sqrt(3 - Math.cos(2 * phi0)));
} else if (theta <= psi0 + psi1) {
Gprime = (2 * Math.PI) / 3 + Math.atan(hprime / Math.sqrt(3)) - xiprime;
F =
Math.PI / 2 -
Math.asin((2 * sin_phi0) / Math.sqrt(3 - Math.cos(2 * phi0)));
} else {
Gprime = xiprime - Math.PI / 6;
F = Math.PI / 4;
}
const aprime =
theta <= psi0
? hprime
: (Math.sqrt(hprime * hprime + 3) *
Math.sin(Math.PI / 3 - Math.atan(hprime / Math.sqrt(3)))) /
Math.sin(xiprime);
const cprime =
theta <= psi0 + psi1 ? Math.sqrt(hprime * hprime + 3) : 3 - hprime;
const x = Math.acos(
Math.cos(r) * Math.cos(c) + Math.sin(r) * Math.sin(c) * Math.cos(beta)
);
const gamma = Math.asin((Math.sin(beta) * Math.sin(r)) / Math.sin(x));
const epsilon = Math.acos(
Math.sin(G) * Math.sin(gamma) * Math.cos(c) - Math.cos(G) * Math.cos(gamma)
);
const upupvp = (gamma + G + epsilon - Math.PI) / (F + G - Math.PI / 2);
const cos_xy = Math.sqrt(
1 - Math.pow((Math.sin(G) * Math.sin(c)) / Math.sin(epsilon), 2)
);
const xpxpyp = Math.sqrt((1 - Math.cos(x)) / (1 - cos_xy));
const uprime = aprime * upupvp;
const xpyp = Math.sqrt(
uprime * uprime + cprime * cprime - 2 * uprime * cprime * Math.cos(Gprime)
);
const cos_gammaprime = Math.sqrt(
1 - Math.pow((uprime * Math.sin(Gprime)) / xpyp, 2)
);
const xprime = xpxpyp * xpyp;
const yprime = xpyp - xprime;
let rprime = Math.sqrt(
xprime * xprime + cprime * cprime - 2 * xprime * cprime * cos_gammaprime
);
const alphaprime = Math.acos(
(yprime * yprime - uprime * uprime - rprime * rprime) /
(-2 * uprime * rprime)
);
let thetaprime;
if (theta <= psi0) thetaprime = alphaprime;
else if (theta <= psi0 + psi1) thetaprime = Math.PI - (xiprime - Math.PI / 6) - alphaprime;
else thetaprime = Math.PI - (xiprime - Math.PI / 6) + alphaprime;
if (lambda < 1e-12) thetaprime = phi < phi0 ? 0 : Math.PI;
if (x === 0) thetaprime = Gprime, rprime = cprime;
const xm = rprime * Math.sin(thetaprime);
const ym = -rprime * Math.cos(thetaprime) - (3 - hprime);
let xy0 = xm;
let xy1 = ym / Math.sqrt(3);
xy0 *= Math.sqrt(3) / 3;
xy1 *= Math.sqrt(3) / 3;
return [xy0, xy1];
}