hammerRetroazimuthal.js

import {geoProjectionMutator as projectionMutator, geoCircle} from "d3-geo";
import {abs, acos, asin, atan2, cos, degrees, epsilon, halfPi, radians, sqrt, sin} from "./math.js";

export function hammerRetroazimuthalRaw(phi0) {
  var sinPhi0 = sin(phi0),
      cosPhi0 = cos(phi0),
      rotate = hammerRetroazimuthalRotation(phi0);

  rotate.invert = hammerRetroazimuthalRotation(-phi0);

  function forward(lambda, phi) {
    var p = rotate(lambda, phi);
    lambda = p[0], phi = p[1];
    var sinPhi = sin(phi),
        cosPhi = cos(phi),
        cosLambda = cos(lambda),
        z = acos(sinPhi0 * sinPhi + cosPhi0 * cosPhi * cosLambda),
        sinz = sin(z),
        K = abs(sinz) > epsilon ? z / sinz : 1;
    return [
      K * cosPhi0 * sin(lambda),
      (abs(lambda) > halfPi ? K : -K) // rotate for back hemisphere
        * (sinPhi0 * cosPhi - cosPhi0 * sinPhi * cosLambda)
    ];
  }

  forward.invert = function(x, y) {
    var rho = sqrt(x * x + y * y),
        sinz = -sin(rho),
        cosz = cos(rho),
        a = rho * cosz,
        b = -y * sinz,
        c = rho * sinPhi0,
        d = sqrt(a * a + b * b - c * c),
        phi = atan2(a * c + b * d, b * c - a * d),
        lambda = (rho > halfPi ? -1 : 1) * atan2(x * sinz, rho * cos(phi) * cosz + y * sin(phi) * sinz);
    return rotate.invert(lambda, phi);
  };

  return forward;
}

// Latitudinal rotation by phi0.
// Temporary hack until D3 supports arbitrary small-circle clipping origins.
function hammerRetroazimuthalRotation(phi0) {
  var sinPhi0 = sin(phi0),
      cosPhi0 = cos(phi0);

  return function(lambda, phi) {
    var cosPhi = cos(phi),
        x = cos(lambda) * cosPhi,
        y = sin(lambda) * cosPhi,
        z = sin(phi);
    return [
      atan2(y, x * cosPhi0 - z * sinPhi0),
      asin(z * cosPhi0 + x * sinPhi0)
    ];
  };
}

export default function() {
  var phi0 = 0,
      m = projectionMutator(hammerRetroazimuthalRaw),
      p = m(phi0),
      rotate_ = p.rotate,
      stream_ = p.stream,
      circle = geoCircle();

  p.parallel = function(_) {
    if (!arguments.length) return phi0 * degrees;
    var r = p.rotate();
    return m(phi0 = _ * radians).rotate(r);
  };

  // Temporary hack; see hammerRetroazimuthalRotation.
  p.rotate = function(_) {
    if (!arguments.length) return (_ = rotate_.call(p), _[1] += phi0 * degrees, _);
    rotate_.call(p, [_[0], _[1] - phi0 * degrees]);
    circle.center([-_[0], -_[1]]);
    return p;
  };

  p.stream = function(stream) {
    stream = stream_(stream);
    stream.sphere = function() {
      stream.polygonStart();
      var epsilon = 1e-2,
          ring = circle.radius(90 - epsilon)().coordinates[0],
          n = ring.length - 1,
          i = -1,
          p;
      stream.lineStart();
      while (++i < n) stream.point((p = ring[i])[0], p[1]);
      stream.lineEnd();
      ring = circle.radius(90 + epsilon)().coordinates[0];
      n = ring.length - 1;
      stream.lineStart();
      while (--i >= 0) stream.point((p = ring[i])[0], p[1]);
      stream.lineEnd();
      stream.polygonEnd();
    };
    return stream;
  };

  return p
      .scale(79.4187)
      .parallel(45)
      .clipAngle(180 - 1e-3);
}