vanDerGrinten.js

import {geoProjection as projection} from "d3-geo";
import {abs, acos, asin, atan, cos, epsilon, halfPi, pi, sign, sin, sqrt, tan} from "./math.js";

export function vanDerGrintenRaw(lambda, phi) {
  if (abs(phi) < epsilon) return [lambda, 0];
  var sinTheta = abs(phi / halfPi),
      theta = asin(sinTheta);
  if (abs(lambda) < epsilon || abs(abs(phi) - halfPi) < epsilon) return [0, sign(phi) * pi * tan(theta / 2)];
  var cosTheta = cos(theta),
      A = abs(pi / lambda - lambda / pi) / 2,
      A2 = A * A,
      G = cosTheta / (sinTheta + cosTheta - 1),
      P = G * (2 / sinTheta - 1),
      P2 = P * P,
      P2_A2 = P2 + A2,
      G_P2 = G - P2,
      Q = A2 + G;
  return [
    sign(lambda) * pi * (A * G_P2 + sqrt(A2 * G_P2 * G_P2 - P2_A2 * (G * G - P2))) / P2_A2,
    sign(phi) * pi * (P * Q - A * sqrt((A2 + 1) * P2_A2 - Q * Q)) / P2_A2
  ];
}

vanDerGrintenRaw.invert = function(x, y) {
  if (abs(y) < epsilon) return [x, 0];
  if (abs(x) < epsilon) return [0, halfPi * sin(2 * atan(y / pi))];
  var x2 = (x /= pi) * x,
      y2 = (y /= pi) * y,
      x2_y2 = x2 + y2,
      z = x2_y2 * x2_y2,
      c1 = -abs(y) * (1 + x2_y2),
      c2 = c1 - 2 * y2 + x2,
      c3 = -2 * c1 + 1 + 2 * y2 + z,
      d = y2 / c3 + (2 * c2 * c2 * c2 / (c3 * c3 * c3) - 9 * c1 * c2 / (c3 * c3)) / 27,
      a1 = (c1 - c2 * c2 / (3 * c3)) / c3,
      m1 = 2 * sqrt(-a1 / 3),
      theta1 = acos(3 * d / (a1 * m1)) / 3;
  return [
    pi * (x2_y2 - 1 + sqrt(1 + 2 * (x2 - y2) + z)) / (2 * x),
    sign(y) * pi * (-m1 * cos(theta1 + pi / 3) - c2 / (3 * c3))
  ];
};

export default function() {
  return projection(vanDerGrintenRaw)
      .scale(79.4183);
}