mpetroff/index.html

## index.html
<!DOCTYPE html>
<html lang="en">
<meta charset="utf-8">
<title>New Projection</title>
<style>
.stroke {
  fill: none;
  stroke: #000;
  stroke-width: 2px;
}

.fill {
  fill: #fff;
}

.graticule {
  fill: none;
  stroke: #666;
  stroke-width: 2px;
  mix-blend-mode: exclusion;
}

.land {
  fill: #222;
}
</style>
<script src="https://cdn.jsdelivr.net/npm/d3@5.16.0/dist/d3.min.js" integrity="sha256-Xb6SSzhH3wEPC4Vy3W70Lqh9Y3Du/3KxPqI2JHQSpTw=" crossorigin="anonymous"></script>
<script src="https://cdn.jsdelivr.net/npm/d3-geo-projection@2.9.0/dist/d3-geo-projection.min.js" integrity="sha256-LY1qVgxfyQnPl8CVKGRKGIfi5BUjq9aMVgxEWC26gT0=" crossorigin="anonymous"></script>
<script src="https://cdn.jsdelivr.net/npm/topojson@3.0.2/dist/topojson.min.js" integrity="sha256-tHoAPGoNdhIR28YHl9DWLzeRfdwigkH7OCBXMrHXhoM=" crossorigin="anonymous"></script>

<script src="newprojection.js"></script>

<svg width="678" height="678" id="small"></svg>
<script>
function projectionRaw(lambda, phi) {
    const sLambda = Math.sign(lambda),
        sPhi = Math.sign(phi),
        cosPhi = Math.cos(phi),
        x = Math.cos(lambda) * cosPhi,
        y = Math.sin(lambda) * cosPhi;
    let z = Math.sin(sPhi * phi);
    lambda = Math.abs(Math.atan2(y, z));
    phi = Math.asin(x);
    if (Math.abs(lambda - Math.PI / 2) > 1e-6) lambda %= Math.PI / 2;
    const point = hexadecant(lambda > Math.PI / 4 ? Math.PI / 2 - lambda : lambda, phi);
    if (lambda > Math.PI / 4) z = point[0], point[0] = -point[1], point[1] = -z;
    return (point[0] *= sLambda, point[1] *= -sPhi, point);
};

function projectionRawInvert(x, y) {
    if (Math.abs(x) > 1)
        x = Math.sign(x) * 2 - x;
    if (Math.abs(y) > 1)
        y = Math.sign(y) * 2 - y;
    const sx = Math.sign(x),
        sy = Math.sign(y),
        x0 = -sx * x,
        y0 = -sy * y,
        t = y0 / x0 < 1,
        p = hexadecantInvert(t ? y0 : x0, t ? x0 : y0),
        lambda = t ? -Math.PI / 2 - p[0] : p[0],
        phi = p[1],
        cosPhi = Math.cos(phi);
    return [sx * (Math.atan2(Math.sin(lambda) * cosPhi, -Math.sin(phi)) + Math.PI), sy * Math.asin(Math.cos(lambda) * cosPhi)];
};

function projectionRaw2(lambda, phi) {
    // Go through inverse to make sure it works
    var tmp = projectionRaw(lambda, phi);
    tmp = projectionRawInvert(tmp[0], tmp[1]);
    return projectionRaw(tmp[0], tmp[1]);
}

projection = d3.geoQuincuncial(projectionRaw2).scale(176.423);

var svg = d3.select("#small"),
    width = +svg.attr("width"),
    height = +svg.attr("height");

var projection = projection
    .scale(239)
    .translate([width / 2, height / 2])
    .precision(0.1);

var path = d3.geoPath()
    .projection(projection);

var graticule = d3.geoGraticule();

svg.append("defs").append("path")
    .datum({type: "Sphere"})
    .attr("id", "sphere")
    .attr("d", path);

svg.append("use")
    .attr("class", "fill")
    .attr("xlink:href", "#sphere");

svg.append("path")
    .datum(graticule)
    .attr("class", "graticule")
    .attr("d", path);

svg.append("use")
    .attr("class", "stroke")
    .attr("xlink:href", "#sphere");

d3.json("https://cdn.jsdelivr.net/npm/world-atlas@2.0.2/land-110m.json").then(world => {
  svg.insert("path", ".graticule")
      .datum(topojson.feature(world, world.objects.land))
      .attr("class", "land")
      .attr("d", path);
});
</script>

## newprojection.js
/*
 * Matthew Petroff, October 2021
 *
 * Implementation of "A Square Equal-Area Map Projection with Low Angular
 * Distortion, Minimal Cusps, and Closed-Form Solutions":
 * https://doi.org/10.1145/3460521
 *
 * This work is placed into the public domain via the CC0 1.0 public domain
 * dedication: https://creativecommons.org/publicdomain/zero/1.0/
 *
 */

function hexadecant(lambda, phi) {
    // Handle edge cases
    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 psi0 = Math.asin(1 / Math.sqrt(2 - cos_phi0 * cos_phi0));
    const psi1 = Math.PI - 2 * psi0;
    const rho = Math.asin(2 * sin_phi0 / Math.sqrt(3 - Math.cos(2 * phi0)));

    const hprime = 12 / Math.PI * (psi0 + rho - Math.PI / 2);
    const xiprime = Math.atan((Math.PI * (hprime - 3) * (hprime - 3)) / (Math.sqrt(3) *
        (Math.PI * (hprime * hprime - 2 * hprime + 45) - 96 * psi0 - 48 * rho)));

    const psi0prime = Math.atan(Math.sqrt(3) / hprime);
    const psi1prime = 7 * Math.PI / 6 - psi0prime - xiprime;
    const psi2prime = xiprime - Math.PI / 6;
    const rhoprime = Math.atan(hprime / Math.sqrt(3));

    const lambda0 = Math.floor((lambda + Math.PI / 4) / (Math.PI / 2)) * Math.PI / 2;

    const theta = Math.abs(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));

    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.SQRT2) : Math.PI / 2 - phi0;

    let G, Gprime, F;
    if (theta <= psi0) {
        G = psi0;
        Gprime = psi0prime;
        F = rho;
    } else if (theta <= psi0 + psi1) {
        G = psi1;
        Gprime = psi1prime;
        F = Math.PI / 2 - rho;
    } else {
        G = psi0;
        Gprime = psi2prime;
        F = Math.PI / 4;
    }

    const aprime = theta <= psi0 ? hprime : Math.sqrt(hprime * hprime + 3) *
        Math.sin(Math.PI / 3 - rhoprime) / 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 = xpyp * xpxpyp;
    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));

    // Put theta back in the correct section
    let thetaprime;
    if (theta <= psi0)
        thetaprime = alphaprime;
    else if (theta <= psi0 + psi1)
        thetaprime = 7 * Math.PI / 6 - xiprime - alphaprime;
    else
        thetaprime = 7 * Math.PI / 6 - xiprime + alphaprime;

    // Handle edge cases
    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);

    const xy0 = xm * Math.sqrt(3) / 3;
    const xy1 = ym / 3;

    return [xy0, xy1];
}

function hexadecantInvert(x_m, y_m) {
    const phi0 = 3 * Math.PI / 8,
        cos_phi0 = Math.cos(phi0),
        sin_phi0 = Math.sin(phi0);

    const x_c = 3 * x_m / Math.sqrt(3),
        y_h = 3 * y_m;

    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 rprime = Math.sqrt(x_c * x_c + (hprime - 3 - y_h) * (hprime - 3 - y_h));
    const thetaprime = Math.abs(Math.atan2(x_c, hprime - 3 - y_h));

    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 psi0prime = Math.atan(Math.sqrt(3) / hprime);
    const psi1prime = 7 * Math.PI / 6 - psi0prime - xiprime;
    const psi2prime = xiprime - Math.PI / 6;

    let alphaprime;
    if (thetaprime <= psi0prime)
        alphaprime = thetaprime;
    else if (thetaprime <= psi0prime + psi1prime)
        alphaprime = Math.PI - psi2prime - thetaprime;
    else
        alphaprime = thetaprime + psi2prime - Math.PI;

    let c = thetaprime <= psi0prime + psi1prime ? Math.acos(cos_phi0 / Math.SQRT2) : Math.PI / 2 - phi0;

    const psi0 = Math.asin(1 / Math.sqrt(2 - cos_phi0 * cos_phi0));
    const psi1 = Math.PI - 2 * psi0;
    const rho = Math.asin(2 * sin_phi0 / Math.sqrt(3 - Math.cos(2 * phi0)));

    // The original publication neglected to mention that the non-prime to
    // prime replacement should only be done in the condition statements, so G
    // is still set to psi0 / psi1.
    let G, Gprime, F;
    if (thetaprime <= psi0prime) {
        G = psi0;
        Gprime = psi0prime;
        F = rho;
    } else if (thetaprime <= psi0prime + psi1prime) {
        G = psi1;
        Gprime = psi1prime;
        F = Math.PI / 2 - rho;
    } else {
        G = psi0;
        Gprime = psi2prime;
        F = Math.PI / 4;
    }

    const aprime = thetaprime <= psi0prime ? hprime : Math.sqrt(hprime * hprime + 3) *
        Math.sin(Math.PI / 3 - Math.atan(hprime / Math.sqrt(3))) / Math.sin(xiprime);
    const cprime = thetaprime <= psi0prime + psi1prime ? Math.sqrt(hprime * hprime + 3) : 3 - hprime;

    let b;
    if (thetaprime <= psi0prime)
        b = Math.PI / 4;
    else if (thetaprime <= psi0prime + psi1prime)
        b = Math.atan(Math.SQRT2 * Math.tan(phi0));
    else
        b = Math.PI / 2 - Math.atan(Math.SQRT2 * Math.tan(phi0));

    const betaprime = Gprime - alphaprime;
    const xprime = Math.sqrt(rprime * rprime + cprime * cprime - 2 * rprime *
        cprime * Math.cos(betaprime));
    const gammaprime = xprime > 0 ? Math.asin(rprime * Math.sin(betaprime) / xprime) : 0;  // Avoid divide by zero
    const epsilonprime = Math.PI - Gprime - gammaprime;
    const yprime = epsilonprime > 0 ? rprime * Math.sin(alphaprime) / Math.sin(epsilonprime) : 0;  // Avoid divide by zero
    const uprime = Math.sqrt(Math.abs(cprime * cprime +
        (xprime + yprime) * (xprime + yprime) -
        2 * cprime * (xprime + yprime) * Math.cos(gammaprime)));  // Take absolute value to avoid negative numbers due numerical precision issues
    const vprime = aprime - uprime;

    const delta = Math.atan(
        (-Math.sin(vprime * (F + G - Math.PI / 2) / aprime))
        / (Math.cos(b) - Math.cos(vprime * (F + G - Math.PI / 2) / aprime))
    );

    const gamma = F - delta;
    const cos_xy = 1 / Math.sqrt(1 + (Math.tan(b) / Math.cos(delta)) * (Math.tan(b) / Math.cos(delta)));

    const x = Math.acos(1 - (xprime / (xprime + yprime)) * (xprime / (xprime + yprime)) * (1 - cos_xy));

    const r = Math.acos(Math.cos(x) * Math.cos(c) + Math.sin(x) * Math.sin(c) * Math.cos(gamma));
    const beta = Math.asin(Math.sin(x) * Math.sin(gamma) / Math.sin(r));

    let alpha;
    if (thetaprime <= psi0prime)
        alpha = psi0 - beta;
    else if (thetaprime <= psi0prime + psi1prime)
        alpha = beta + psi0;  // There was a sign error here in the original publication
    else
        alpha = Math.PI - beta;

    const phi = Math.asin(sin_phi0 * Math.cos(r) - cos_phi0 * Math.sin(r) * Math.cos(alpha));
    let lambda = Math.sign(x_c) * Math.atan(Math.sin(alpha) *
        Math.sin(r) * cos_phi0 / (Math.cos(r) - sin_phi0 * Math.sin(phi)));
    lambda *= Math.abs(phi) != Math.PI / 2  // Handle edge case

    return [lambda, phi];
}
	<!DOCTYPE html>
	<html lang="en">
	<meta charset="utf-8">
	<title>New Projection</title>
	<style>
	.stroke {
	fill: none;
	stroke: #000;
	stroke-width: 2px;
	}

	.fill {
	fill: #fff;
	}

	.graticule {
	fill: none;
	stroke: #666;
	stroke-width: 2px;
	mix-blend-mode: exclusion;
	}

	.land {
	fill: #222;
	}
	</style>
	<script src="https://cdn.jsdelivr.net/npm/d3@5.16.0/dist/d3.min.js" integrity="sha256-Xb6SSzhH3wEPC4Vy3W70Lqh9Y3Du/3KxPqI2JHQSpTw=" crossorigin="anonymous"></script>
	<script src="https://cdn.jsdelivr.net/npm/d3-geo-projection@2.9.0/dist/d3-geo-projection.min.js" integrity="sha256-LY1qVgxfyQnPl8CVKGRKGIfi5BUjq9aMVgxEWC26gT0=" crossorigin="anonymous"></script>
	<script src="https://cdn.jsdelivr.net/npm/topojson@3.0.2/dist/topojson.min.js" integrity="sha256-tHoAPGoNdhIR28YHl9DWLzeRfdwigkH7OCBXMrHXhoM=" crossorigin="anonymous"></script>

	<script src="newprojection.js"></script>

	<svg width="678" height="678" id="small"></svg>
	<script>
	function projectionRaw(lambda, phi) {
	const sLambda = Math.sign(lambda),
	sPhi = Math.sign(phi),
	cosPhi = Math.cos(phi),
	x = Math.cos(lambda) * cosPhi,
	y = Math.sin(lambda) * cosPhi;
	let z = Math.sin(sPhi * phi);
	lambda = Math.abs(Math.atan2(y, z));
	phi = Math.asin(x);
	if (Math.abs(lambda - Math.PI / 2) > 1e-6) lambda %= Math.PI / 2;
	const point = hexadecant(lambda > Math.PI / 4 ? Math.PI / 2 - lambda : lambda, phi);
	if (lambda > Math.PI / 4) z = point[0], point[0] = -point[1], point[1] = -z;
	return (point[0] = sLambda, point[1] = -sPhi, point);
	};

	function projectionRawInvert(x, y) {
	if (Math.abs(x) > 1)
	x = Math.sign(x) * 2 - x;
	if (Math.abs(y) > 1)
	y = Math.sign(y) * 2 - y;
	const sx = Math.sign(x),
	sy = Math.sign(y),
	x0 = -sx * x,
	y0 = -sy * y,
	t = y0 / x0 < 1,
	p = hexadecantInvert(t ? y0 : x0, t ? x0 : y0),
	lambda = t ? -Math.PI / 2 - p[0] : p[0],
	phi = p[1],
	cosPhi = Math.cos(phi);
	return [sx * (Math.atan2(Math.sin(lambda) * cosPhi, -Math.sin(phi)) + Math.PI), sy * Math.asin(Math.cos(lambda) * cosPhi)];
	};

	function projectionRaw2(lambda, phi) {
	// Go through inverse to make sure it works
	var tmp = projectionRaw(lambda, phi);
	tmp = projectionRawInvert(tmp[0], tmp[1]);
	return projectionRaw(tmp[0], tmp[1]);
	}

	projection = d3.geoQuincuncial(projectionRaw2).scale(176.423);

	var svg = d3.select("#small"),
	width = +svg.attr("width"),
	height = +svg.attr("height");

	var projection = projection
	.scale(239)
	.translate([width / 2, height / 2])
	.precision(0.1);

	var path = d3.geoPath()
	.projection(projection);

	var graticule = d3.geoGraticule();

	svg.append("defs").append("path")
	.datum({type: "Sphere"})
	.attr("id", "sphere")
	.attr("d", path);

	svg.append("use")
	.attr("class", "fill")
	.attr("xlink:href", "#sphere");

	svg.append("path")
	.datum(graticule)
	.attr("class", "graticule")
	.attr("d", path);

	svg.append("use")
	.attr("class", "stroke")
	.attr("xlink:href", "#sphere");

	d3.json("https://cdn.jsdelivr.net/npm/world-atlas@2.0.2/land-110m.json").then(world => {
	svg.insert("path", ".graticule")
	.datum(topojson.feature(world, world.objects.land))
	.attr("class", "land")
	.attr("d", path);
	});
	</script>
	/*
	* Matthew Petroff, October 2021
	*
	* Implementation of "A Square Equal-Area Map Projection with Low Angular
	* Distortion, Minimal Cusps, and Closed-Form Solutions":
	* https://doi.org/10.1145/3460521
	*
	* This work is placed into the public domain via the CC0 1.0 public domain
	* dedication: https://creativecommons.org/publicdomain/zero/1.0/
	*
	*/

	function hexadecant(lambda, phi) {
	// Handle edge cases
	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 psi0 = Math.asin(1 / Math.sqrt(2 - cos_phi0 * cos_phi0));
	const psi1 = Math.PI - 2 * psi0;
	const rho = Math.asin(2 * sin_phi0 / Math.sqrt(3 - Math.cos(2 * phi0)));

	const hprime = 12 / Math.PI * (psi0 + rho - Math.PI / 2);
	const xiprime = Math.atan((Math.PI * (hprime - 3) * (hprime - 3)) / (Math.sqrt(3) *
	(Math.PI * (hprime * hprime - 2 * hprime + 45) - 96 * psi0 - 48 * rho)));

	const psi0prime = Math.atan(Math.sqrt(3) / hprime);
	const psi1prime = 7 * Math.PI / 6 - psi0prime - xiprime;
	const psi2prime = xiprime - Math.PI / 6;
	const rhoprime = Math.atan(hprime / Math.sqrt(3));

	const lambda0 = Math.floor((lambda + Math.PI / 4) / (Math.PI / 2)) * Math.PI / 2;

	const theta = Math.abs(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));

	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.SQRT2) : Math.PI / 2 - phi0;

	let G, Gprime, F;
	if (theta <= psi0) {
	G = psi0;
	Gprime = psi0prime;
	F = rho;
	} else if (theta <= psi0 + psi1) {
	G = psi1;
	Gprime = psi1prime;
	F = Math.PI / 2 - rho;
	} else {
	G = psi0;
	Gprime = psi2prime;
	F = Math.PI / 4;
	}

	const aprime = theta <= psi0 ? hprime : Math.sqrt(hprime * hprime + 3) *
	Math.sin(Math.PI / 3 - rhoprime) / 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 = xpyp * xpxpyp;
	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));

	// Put theta back in the correct section
	let thetaprime;
	if (theta <= psi0)
	thetaprime = alphaprime;
	else if (theta <= psi0 + psi1)
	thetaprime = 7 * Math.PI / 6 - xiprime - alphaprime;
	else
	thetaprime = 7 * Math.PI / 6 - xiprime + alphaprime;

	// Handle edge cases
	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);

	const xy0 = xm * Math.sqrt(3) / 3;
	const xy1 = ym / 3;

	return [xy0, xy1];
	}

	function hexadecantInvert(x_m, y_m) {
	const phi0 = 3 * Math.PI / 8,
	cos_phi0 = Math.cos(phi0),
	sin_phi0 = Math.sin(phi0);

	const x_c = 3 * x_m / Math.sqrt(3),
	y_h = 3 * y_m;

	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 rprime = Math.sqrt(x_c * x_c + (hprime - 3 - y_h) * (hprime - 3 - y_h));
	const thetaprime = Math.abs(Math.atan2(x_c, hprime - 3 - y_h));

	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 psi0prime = Math.atan(Math.sqrt(3) / hprime);
	const psi1prime = 7 * Math.PI / 6 - psi0prime - xiprime;
	const psi2prime = xiprime - Math.PI / 6;

	let alphaprime;
	if (thetaprime <= psi0prime)
	alphaprime = thetaprime;
	else if (thetaprime <= psi0prime + psi1prime)
	alphaprime = Math.PI - psi2prime - thetaprime;
	else
	alphaprime = thetaprime + psi2prime - Math.PI;

	let c = thetaprime <= psi0prime + psi1prime ? Math.acos(cos_phi0 / Math.SQRT2) : Math.PI / 2 - phi0;

	const psi0 = Math.asin(1 / Math.sqrt(2 - cos_phi0 * cos_phi0));
	const psi1 = Math.PI - 2 * psi0;
	const rho = Math.asin(2 * sin_phi0 / Math.sqrt(3 - Math.cos(2 * phi0)));

	// The original publication neglected to mention that the non-prime to
	// prime replacement should only be done in the condition statements, so G
	// is still set to psi0 / psi1.
	let G, Gprime, F;
	if (thetaprime <= psi0prime) {
	G = psi0;
	Gprime = psi0prime;
	F = rho;
	} else if (thetaprime <= psi0prime + psi1prime) {
	G = psi1;
	Gprime = psi1prime;
	F = Math.PI / 2 - rho;
	} else {
	G = psi0;
	Gprime = psi2prime;
	F = Math.PI / 4;
	}

	const aprime = thetaprime <= psi0prime ? hprime : Math.sqrt(hprime * hprime + 3) *
	Math.sin(Math.PI / 3 - Math.atan(hprime / Math.sqrt(3))) / Math.sin(xiprime);
	const cprime = thetaprime <= psi0prime + psi1prime ? Math.sqrt(hprime * hprime + 3) : 3 - hprime;

	let b;
	if (thetaprime <= psi0prime)
	b = Math.PI / 4;
	else if (thetaprime <= psi0prime + psi1prime)
	b = Math.atan(Math.SQRT2 * Math.tan(phi0));
	else
	b = Math.PI / 2 - Math.atan(Math.SQRT2 * Math.tan(phi0));

	const betaprime = Gprime - alphaprime;
	const xprime = Math.sqrt(rprime * rprime + cprime * cprime - 2 * rprime *
	cprime * Math.cos(betaprime));
	const gammaprime = xprime > 0 ? Math.asin(rprime * Math.sin(betaprime) / xprime) : 0; // Avoid divide by zero
	const epsilonprime = Math.PI - Gprime - gammaprime;
	const yprime = epsilonprime > 0 ? rprime * Math.sin(alphaprime) / Math.sin(epsilonprime) : 0; // Avoid divide by zero
	const uprime = Math.sqrt(Math.abs(cprime * cprime +
	(xprime + yprime) * (xprime + yprime) -
	2 * cprime * (xprime + yprime) * Math.cos(gammaprime))); // Take absolute value to avoid negative numbers due numerical precision issues
	const vprime = aprime - uprime;

	const delta = Math.atan(
	(-Math.sin(vprime * (F + G - Math.PI / 2) / aprime))
	/ (Math.cos(b) - Math.cos(vprime * (F + G - Math.PI / 2) / aprime))
	);

	const gamma = F - delta;
	const cos_xy = 1 / Math.sqrt(1 + (Math.tan(b) / Math.cos(delta)) * (Math.tan(b) / Math.cos(delta)));

	const x = Math.acos(1 - (xprime / (xprime + yprime)) * (xprime / (xprime + yprime)) * (1 - cos_xy));

	const r = Math.acos(Math.cos(x) * Math.cos(c) + Math.sin(x) * Math.sin(c) * Math.cos(gamma));
	const beta = Math.asin(Math.sin(x) * Math.sin(gamma) / Math.sin(r));

	let alpha;
	if (thetaprime <= psi0prime)
	alpha = psi0 - beta;
	else if (thetaprime <= psi0prime + psi1prime)
	alpha = beta + psi0; // There was a sign error here in the original publication
	else
	alpha = Math.PI - beta;

	const phi = Math.asin(sin_phi0 * Math.cos(r) - cos_phi0 * Math.sin(r) * Math.cos(alpha));
	let lambda = Math.sign(x_c) * Math.atan(Math.sin(alpha) *
	Math.sin(r) * cos_phi0 / (Math.cos(r) - sin_phi0 * Math.sin(phi)));
	lambda *= Math.abs(phi) != Math.PI / 2 // Handle edge case

	return [lambda, phi];
	}