|
// https://d3js.org/d3-geo-projection/ Version 2.1.0. Copyright 2017 Mike Bostock. |
|
(function (global, factory) { |
|
typeof exports === 'object' && typeof module !== 'undefined' ? factory(exports, require('d3-geo'), require('d3-array')) : |
|
typeof define === 'function' && define.amd ? define(['exports', 'd3-geo', 'd3-array'], factory) : |
|
(factory((global.d3 = global.d3 || {}),global.d3,global.d3)); |
|
}(this, (function (exports,d3Geo,d3Array) { 'use strict'; |
|
|
|
var abs = Math.abs; |
|
var atan = Math.atan; |
|
var atan2 = Math.atan2; |
|
|
|
var cos = Math.cos; |
|
var exp = Math.exp; |
|
var floor = Math.floor; |
|
var log = Math.log; |
|
var max = Math.max; |
|
var min = Math.min; |
|
var pow = Math.pow; |
|
var round = Math.round; |
|
var sign = Math.sign || function(x) { return x > 0 ? 1 : x < 0 ? -1 : 0; }; |
|
var sin = Math.sin; |
|
var tan = Math.tan; |
|
|
|
var epsilon = 1e-6; |
|
var epsilon2 = 1e-12; |
|
var pi = Math.PI; |
|
var halfPi = pi / 2; |
|
var quarterPi = pi / 4; |
|
var sqrt1_2 = Math.SQRT1_2; |
|
var sqrt2 = sqrt(2); |
|
var sqrtPi = sqrt(pi); |
|
var tau = pi * 2; |
|
var degrees = 180 / pi; |
|
var radians = pi / 180; |
|
|
|
function sinci(x) { |
|
return x ? x / Math.sin(x) : 1; |
|
} |
|
|
|
function asin(x) { |
|
return x > 1 ? halfPi : x < -1 ? -halfPi : Math.asin(x); |
|
} |
|
|
|
function acos(x) { |
|
return x > 1 ? 0 : x < -1 ? pi : Math.acos(x); |
|
} |
|
|
|
function sqrt(x) { |
|
return x > 0 ? Math.sqrt(x) : 0; |
|
} |
|
|
|
function tanh(x) { |
|
x = exp(2 * x); |
|
return (x - 1) / (x + 1); |
|
} |
|
|
|
function sinh(x) { |
|
return (exp(x) - exp(-x)) / 2; |
|
} |
|
|
|
function cosh(x) { |
|
return (exp(x) + exp(-x)) / 2; |
|
} |
|
|
|
function arsinh(x) { |
|
return log(x + sqrt(x * x + 1)); |
|
} |
|
|
|
function arcosh(x) { |
|
return log(x + sqrt(x * x - 1)); |
|
} |
|
|
|
function airyRaw(beta) { |
|
var tanBeta_2 = tan(beta / 2), |
|
b = 2 * log(cos(beta / 2)) / (tanBeta_2 * tanBeta_2); |
|
|
|
function forward(x, y) { |
|
var cosx = cos(x), |
|
cosy = cos(y), |
|
siny = sin(y), |
|
cosz = cosy * cosx, |
|
k = -((1 - cosz ? log((1 + cosz) / 2) / (1 - cosz) : -0.5) + b / (1 + cosz)); |
|
return [k * cosy * sin(x), k * siny]; |
|
} |
|
|
|
forward.invert = function(x, y) { |
|
var r = sqrt(x * x + y * y), |
|
z = -beta / 2, |
|
i = 50, delta; |
|
if (!r) return [0, 0]; |
|
do { |
|
var z_2 = z / 2, |
|
cosz_2 = cos(z_2), |
|
sinz_2 = sin(z_2), |
|
tanz_2 = tan(z_2), |
|
lnsecz_2 = log(1 / cosz_2); |
|
z -= delta = (2 / tanz_2 * lnsecz_2 - b * tanz_2 - r) / (-lnsecz_2 / (sinz_2 * sinz_2) + 1 - b / (2 * cosz_2 * cosz_2)); |
|
} while (abs(delta) > epsilon && --i > 0); |
|
var sinz = sin(z); |
|
return [atan2(x * sinz, r * cos(z)), asin(y * sinz / r)]; |
|
}; |
|
|
|
return forward; |
|
} |
|
|
|
var airy = function() { |
|
var beta = halfPi, |
|
m = d3Geo.geoProjectionMutator(airyRaw), |
|
p = m(beta); |
|
|
|
p.radius = function(_) { |
|
return arguments.length ? m(beta = _ * radians) : beta * degrees; |
|
}; |
|
|
|
return p |
|
.scale(179.976) |
|
.clipAngle(147); |
|
}; |
|
|
|
function aitoffRaw(x, y) { |
|
var cosy = cos(y), sincia = sinci(acos(cosy * cos(x /= 2))); |
|
return [2 * cosy * sin(x) * sincia, sin(y) * sincia]; |
|
} |
|
|
|
// Abort if [x, y] is not within an ellipse centered at [0, 0] with |
|
// semi-major axis pi and semi-minor axis pi/2. |
|
aitoffRaw.invert = function(x, y) { |
|
if (x * x + 4 * y * y > pi * pi + epsilon) return; |
|
var x1 = x, y1 = y, i = 25; |
|
do { |
|
var sinx = sin(x1), |
|
sinx_2 = sin(x1 / 2), |
|
cosx_2 = cos(x1 / 2), |
|
siny = sin(y1), |
|
cosy = cos(y1), |
|
sin_2y = sin(2 * y1), |
|
sin2y = siny * siny, |
|
cos2y = cosy * cosy, |
|
sin2x_2 = sinx_2 * sinx_2, |
|
c = 1 - cos2y * cosx_2 * cosx_2, |
|
e = c ? acos(cosy * cosx_2) * sqrt(f = 1 / c) : f = 0, |
|
f, |
|
fx = 2 * e * cosy * sinx_2 - x, |
|
fy = e * siny - y, |
|
dxdx = f * (cos2y * sin2x_2 + e * cosy * cosx_2 * sin2y), |
|
dxdy = f * (0.5 * sinx * sin_2y - e * 2 * siny * sinx_2), |
|
dydx = f * 0.25 * (sin_2y * sinx_2 - e * siny * cos2y * sinx), |
|
dydy = f * (sin2y * cosx_2 + e * sin2x_2 * cosy), |
|
z = dxdy * dydx - dydy * dxdx; |
|
if (!z) break; |
|
var dx = (fy * dxdy - fx * dydy) / z, |
|
dy = (fx * dydx - fy * dxdx) / z; |
|
x1 -= dx, y1 -= dy; |
|
} while ((abs(dx) > epsilon || abs(dy) > epsilon) && --i > 0); |
|
return [x1, y1]; |
|
}; |
|
|
|
var aitoff = function() { |
|
return d3Geo.geoProjection(aitoffRaw) |
|
.scale(152.63); |
|
}; |
|
|
|
function armadilloRaw(phi0) { |
|
var sinPhi0 = sin(phi0), |
|
cosPhi0 = cos(phi0), |
|
sPhi0 = phi0 >= 0 ? 1 : -1, |
|
tanPhi0 = tan(sPhi0 * phi0), |
|
k = (1 + sinPhi0 - cosPhi0) / 2; |
|
|
|
function forward(lambda, phi) { |
|
var cosPhi = cos(phi), |
|
cosLambda = cos(lambda /= 2); |
|
return [ |
|
(1 + cosPhi) * sin(lambda), |
|
(sPhi0 * phi > -atan2(cosLambda, tanPhi0) - 1e-3 ? 0 : -sPhi0 * 10) + k + sin(phi) * cosPhi0 - (1 + cosPhi) * sinPhi0 * cosLambda // TODO D3 core should allow null or [NaN, NaN] to be returned. |
|
]; |
|
} |
|
|
|
forward.invert = function(x, y) { |
|
var lambda = 0, |
|
phi = 0, |
|
i = 50; |
|
do { |
|
var cosLambda = cos(lambda), |
|
sinLambda = sin(lambda), |
|
cosPhi = cos(phi), |
|
sinPhi = sin(phi), |
|
A = 1 + cosPhi, |
|
fx = A * sinLambda - x, |
|
fy = k + sinPhi * cosPhi0 - A * sinPhi0 * cosLambda - y, |
|
dxdLambda = A * cosLambda / 2, |
|
dxdPhi = -sinLambda * sinPhi, |
|
dydLambda = sinPhi0 * A * sinLambda / 2, |
|
dydPhi = cosPhi0 * cosPhi + sinPhi0 * cosLambda * sinPhi, |
|
denominator = dxdPhi * dydLambda - dydPhi * dxdLambda, |
|
dLambda = (fy * dxdPhi - fx * dydPhi) / denominator / 2, |
|
dPhi = (fx * dydLambda - fy * dxdLambda) / denominator; |
|
lambda -= dLambda, phi -= dPhi; |
|
} while ((abs(dLambda) > epsilon || abs(dPhi) > epsilon) && --i > 0); |
|
return sPhi0 * phi > -atan2(cos(lambda), tanPhi0) - 1e-3 ? [lambda * 2, phi] : null; |
|
}; |
|
|
|
return forward; |
|
} |
|
|
|
var armadillo = function() { |
|
var phi0 = 20 * radians, |
|
sPhi0 = phi0 >= 0 ? 1 : -1, |
|
tanPhi0 = tan(sPhi0 * phi0), |
|
m = d3Geo.geoProjectionMutator(armadilloRaw), |
|
p = m(phi0), |
|
stream_ = p.stream; |
|
|
|
p.parallel = function(_) { |
|
if (!arguments.length) return phi0 * degrees; |
|
tanPhi0 = tan((sPhi0 = (phi0 = _ * radians) >= 0 ? 1 : -1) * phi0); |
|
return m(phi0); |
|
}; |
|
|
|
p.stream = function(stream) { |
|
var rotate = p.rotate(), |
|
rotateStream = stream_(stream), |
|
sphereStream = (p.rotate([0, 0]), stream_(stream)); |
|
p.rotate(rotate); |
|
rotateStream.sphere = function() { |
|
sphereStream.polygonStart(), sphereStream.lineStart(); |
|
for (var lambda = sPhi0 * -180; sPhi0 * lambda < 180; lambda += sPhi0 * 90) sphereStream.point(lambda, sPhi0 * 90); |
|
while (sPhi0 * (lambda -= phi0) >= -180) { // TODO precision? |
|
sphereStream.point(lambda, sPhi0 * -atan2(cos(lambda * radians / 2), tanPhi0) * degrees); |
|
} |
|
sphereStream.lineEnd(), sphereStream.polygonEnd(); |
|
}; |
|
return rotateStream; |
|
}; |
|
|
|
return p |
|
.scale(218.695) |
|
.center([0, 28.0974]); |
|
}; |
|
|
|
function augustRaw(lambda, phi) { |
|
var tanPhi = tan(phi / 2), |
|
k = sqrt(1 - tanPhi * tanPhi), |
|
c = 1 + k * cos(lambda /= 2), |
|
x = sin(lambda) * k / c, |
|
y = tanPhi / c, |
|
x2 = x * x, |
|
y2 = y * y; |
|
return [ |
|
4 / 3 * x * (3 + x2 - 3 * y2), |
|
4 / 3 * y * (3 + 3 * x2 - y2) |
|
]; |
|
} |
|
|
|
augustRaw.invert = function(x, y) { |
|
x *= 3 / 8, y *= 3 / 8; |
|
if (!x && abs(y) > 1) return null; |
|
var x2 = x * x, |
|
y2 = y * y, |
|
s = 1 + x2 + y2, |
|
sin3Eta = sqrt((s - sqrt(s * s - 4 * y * y)) / 2), |
|
eta = asin(sin3Eta) / 3, |
|
xi = sin3Eta ? arcosh(abs(y / sin3Eta)) / 3 : arsinh(abs(x)) / 3, |
|
cosEta = cos(eta), |
|
coshXi = cosh(xi), |
|
d = coshXi * coshXi - cosEta * cosEta; |
|
return [ |
|
sign(x) * 2 * atan2(sinh(xi) * cosEta, 0.25 - d), |
|
sign(y) * 2 * atan2(coshXi * sin(eta), 0.25 + d) |
|
]; |
|
}; |
|
|
|
var august = function() { |
|
return d3Geo.geoProjection(augustRaw) |
|
.scale(66.1603); |
|
}; |
|
|
|
var sqrt8 = sqrt(8); |
|
var phi0 = log(1 + sqrt2); |
|
|
|
function bakerRaw(lambda, phi) { |
|
var phi0 = abs(phi); |
|
return phi0 < quarterPi |
|
? [lambda, log(tan(quarterPi + phi / 2))] |
|
: [lambda * cos(phi0) * (2 * sqrt2 - 1 / sin(phi0)), sign(phi) * (2 * sqrt2 * (phi0 - quarterPi) - log(tan(phi0 / 2)))]; |
|
} |
|
|
|
bakerRaw.invert = function(x, y) { |
|
if ((y0 = abs(y)) < phi0) return [x, 2 * atan(exp(y)) - halfPi]; |
|
var phi = quarterPi, i = 25, delta, y0; |
|
do { |
|
var cosPhi_2 = cos(phi / 2), tanPhi_2 = tan(phi / 2); |
|
phi -= delta = (sqrt8 * (phi - quarterPi) - log(tanPhi_2) - y0) / (sqrt8 - cosPhi_2 * cosPhi_2 / (2 * tanPhi_2)); |
|
} while (abs(delta) > epsilon2 && --i > 0); |
|
return [x / (cos(phi) * (sqrt8 - 1 / sin(phi))), sign(y) * phi]; |
|
}; |
|
|
|
var baker = function() { |
|
return d3Geo.geoProjection(bakerRaw) |
|
.scale(112.314); |
|
}; |
|
|
|
function berghausRaw(lobes) { |
|
var k = 2 * pi / lobes; |
|
|
|
function forward(lambda, phi) { |
|
var p = d3Geo.geoAzimuthalEquidistantRaw(lambda, phi); |
|
if (abs(lambda) > halfPi) { // back hemisphere |
|
var theta = atan2(p[1], p[0]), |
|
r = sqrt(p[0] * p[0] + p[1] * p[1]), |
|
theta0 = k * round((theta - halfPi) / k) + halfPi, |
|
alpha = atan2(sin(theta -= theta0), 2 - cos(theta)); // angle relative to lobe end |
|
theta = theta0 + asin(pi / r * sin(alpha)) - alpha; |
|
p[0] = r * cos(theta); |
|
p[1] = r * sin(theta); |
|
} |
|
return p; |
|
} |
|
|
|
forward.invert = function(x, y) { |
|
var r = sqrt(x * x + y * y); |
|
if (r > halfPi) { |
|
var theta = atan2(y, x), |
|
theta0 = k * round((theta - halfPi) / k) + halfPi, |
|
s = theta > theta0 ? -1 : 1, |
|
A = r * cos(theta0 - theta), |
|
cotAlpha = 1 / tan(s * acos((A - pi) / sqrt(pi * (pi - 2 * A) + r * r))); |
|
theta = theta0 + 2 * atan((cotAlpha + s * sqrt(cotAlpha * cotAlpha - 3)) / 3); |
|
x = r * cos(theta), y = r * sin(theta); |
|
} |
|
return d3Geo.geoAzimuthalEquidistantRaw.invert(x, y); |
|
}; |
|
|
|
return forward; |
|
} |
|
|
|
var berghaus = function() { |
|
var lobes = 5, |
|
m = d3Geo.geoProjectionMutator(berghausRaw), |
|
p = m(lobes), |
|
projectionStream = p.stream, |
|
epsilon$$1 = 1e-2, |
|
cr = -cos(epsilon$$1 * radians), |
|
sr = sin(epsilon$$1 * radians); |
|
|
|
p.lobes = function(_) { |
|
return arguments.length ? m(lobes = +_) : lobes; |
|
}; |
|
|
|
p.stream = function(stream) { |
|
var rotate = p.rotate(), |
|
rotateStream = projectionStream(stream), |
|
sphereStream = (p.rotate([0, 0]), projectionStream(stream)); |
|
p.rotate(rotate); |
|
rotateStream.sphere = function() { |
|
sphereStream.polygonStart(), sphereStream.lineStart(); |
|
for (var i = 0, delta = 360 / lobes, delta0 = 2 * pi / lobes, phi = 90 - 180 / lobes, phi0 = halfPi; i < lobes; ++i, phi -= delta, phi0 -= delta0) { |
|
sphereStream.point(atan2(sr * cos(phi0), cr) * degrees, asin(sr * sin(phi0)) * degrees); |
|
if (phi < -90) { |
|
sphereStream.point(-90, -180 - phi - epsilon$$1); |
|
sphereStream.point(-90, -180 - phi + epsilon$$1); |
|
} else { |
|
sphereStream.point(90, phi + epsilon$$1); |
|
sphereStream.point(90, phi - epsilon$$1); |
|
} |
|
} |
|
sphereStream.lineEnd(), sphereStream.polygonEnd(); |
|
}; |
|
return rotateStream; |
|
}; |
|
|
|
return p |
|
.scale(87.8076) |
|
.center([0, 17.1875]) |
|
.clipAngle(180 - 1e-3); |
|
}; |
|
|
|
function mollweideBromleyTheta(cp, phi) { |
|
var cpsinPhi = cp * sin(phi), i = 30, delta; |
|
do phi -= delta = (phi + sin(phi) - cpsinPhi) / (1 + cos(phi)); |
|
while (abs(delta) > epsilon && --i > 0); |
|
return phi / 2; |
|
} |
|
|
|
function mollweideBromleyRaw(cx, cy, cp) { |
|
|
|
function forward(lambda, phi) { |
|
return [cx * lambda * cos(phi = mollweideBromleyTheta(cp, phi)), cy * sin(phi)]; |
|
} |
|
|
|
forward.invert = function(x, y) { |
|
return y = asin(y / cy), [x / (cx * cos(y)), asin((2 * y + sin(2 * y)) / cp)]; |
|
}; |
|
|
|
return forward; |
|
} |
|
|
|
var mollweideRaw = mollweideBromleyRaw(sqrt2 / halfPi, sqrt2, pi); |
|
|
|
var mollweide = function() { |
|
return d3Geo.geoProjection(mollweideRaw) |
|
.scale(169.529); |
|
}; |
|
|
|
var k = 2.00276; |
|
var w = 1.11072; |
|
|
|
function boggsRaw(lambda, phi) { |
|
var theta = mollweideBromleyTheta(pi, phi); |
|
return [k * lambda / (1 / cos(phi) + w / cos(theta)), (phi + sqrt2 * sin(theta)) / k]; |
|
} |
|
|
|
boggsRaw.invert = function(x, y) { |
|
var ky = k * y, theta = y < 0 ? -quarterPi : quarterPi, i = 25, delta, phi; |
|
do { |
|
phi = ky - sqrt2 * sin(theta); |
|
theta -= delta = (sin(2 * theta) + 2 * theta - pi * sin(phi)) / (2 * cos(2 * theta) + 2 + pi * cos(phi) * sqrt2 * cos(theta)); |
|
} while (abs(delta) > epsilon && --i > 0); |
|
phi = ky - sqrt2 * sin(theta); |
|
return [x * (1 / cos(phi) + w / cos(theta)) / k, phi]; |
|
}; |
|
|
|
var boggs = function() { |
|
return d3Geo.geoProjection(boggsRaw) |
|
.scale(160.857); |
|
}; |
|
|
|
var parallel1 = function(projectAt) { |
|
var phi0 = 0, |
|
m = d3Geo.geoProjectionMutator(projectAt), |
|
p = m(phi0); |
|
|
|
p.parallel = function(_) { |
|
return arguments.length ? m(phi0 = _ * radians) : phi0 * degrees; |
|
}; |
|
|
|
return p; |
|
}; |
|
|
|
function sinusoidalRaw(lambda, phi) { |
|
return [lambda * cos(phi), phi]; |
|
} |
|
|
|
sinusoidalRaw.invert = function(x, y) { |
|
return [x / cos(y), y]; |
|
}; |
|
|
|
var sinusoidal = function() { |
|
return d3Geo.geoProjection(sinusoidalRaw) |
|
.scale(152.63); |
|
}; |
|
|
|
function bonneRaw(phi0) { |
|
if (!phi0) return sinusoidalRaw; |
|
var cotPhi0 = 1 / tan(phi0); |
|
|
|
function forward(lambda, phi) { |
|
var rho = cotPhi0 + phi0 - phi, |
|
e = rho ? lambda * cos(phi) / rho : rho; |
|
return [rho * sin(e), cotPhi0 - rho * cos(e)]; |
|
} |
|
|
|
forward.invert = function(x, y) { |
|
var rho = sqrt(x * x + (y = cotPhi0 - y) * y), |
|
phi = cotPhi0 + phi0 - rho; |
|
return [rho / cos(phi) * atan2(x, y), phi]; |
|
}; |
|
|
|
return forward; |
|
} |
|
|
|
var bonne = function() { |
|
return parallel1(bonneRaw) |
|
.scale(123.082) |
|
.center([0, 26.1441]) |
|
.parallel(45); |
|
}; |
|
|
|
function bottomleyRaw(sinPsi) { |
|
|
|
function forward(lambda, phi) { |
|
var rho = halfPi - phi, |
|
eta = rho ? lambda * sinPsi * sin(rho) / rho : rho; |
|
return [rho * sin(eta) / sinPsi, halfPi - rho * cos(eta)]; |
|
} |
|
|
|
forward.invert = function(x, y) { |
|
var x1 = x * sinPsi, |
|
y1 = halfPi - y, |
|
rho = sqrt(x1 * x1 + y1 * y1), |
|
eta = atan2(x1, y1); |
|
return [(rho ? rho / sin(rho) : 1) * eta / sinPsi, halfPi - rho]; |
|
}; |
|
|
|
return forward; |
|
} |
|
|
|
var bottomley = function() { |
|
var sinPsi = 0.5, |
|
m = d3Geo.geoProjectionMutator(bottomleyRaw), |
|
p = m(sinPsi); |
|
|
|
p.fraction = function(_) { |
|
return arguments.length ? m(sinPsi = +_) : sinPsi; |
|
}; |
|
|
|
return p |
|
.scale(158.837); |
|
}; |
|
|
|
var bromleyRaw = mollweideBromleyRaw(1, 4 / pi, pi); |
|
|
|
var bromley = function() { |
|
return d3Geo.geoProjection(bromleyRaw) |
|
.scale(152.63); |
|
}; |
|
|
|
// Azimuthal distance. |
|
function distance(dPhi, c1, s1, c2, s2, dLambda) { |
|
var cosdLambda = cos(dLambda), r; |
|
if (abs(dPhi) > 1 || abs(dLambda) > 1) { |
|
r = acos(s1 * s2 + c1 * c2 * cosdLambda); |
|
} else { |
|
var sindPhi = sin(dPhi / 2), sindLambda = sin(dLambda / 2); |
|
r = 2 * asin(sqrt(sindPhi * sindPhi + c1 * c2 * sindLambda * sindLambda)); |
|
} |
|
return abs(r) > epsilon ? [r, atan2(c2 * sin(dLambda), c1 * s2 - s1 * c2 * cosdLambda)] : [0, 0]; |
|
} |
|
|
|
// Angle opposite a, and contained between sides of lengths b and c. |
|
function angle(b, c, a) { |
|
return acos((b * b + c * c - a * a) / (2 * b * c)); |
|
} |
|
|
|
// Normalize longitude. |
|
function longitude(lambda) { |
|
return lambda - 2 * pi * floor((lambda + pi) / (2 * pi)); |
|
} |
|
|
|
function chamberlinRaw(p0, p1, p2) { |
|
var points = [ |
|
[p0[0], p0[1], sin(p0[1]), cos(p0[1])], |
|
[p1[0], p1[1], sin(p1[1]), cos(p1[1])], |
|
[p2[0], p2[1], sin(p2[1]), cos(p2[1])] |
|
]; |
|
|
|
for (var a = points[2], b, i = 0; i < 3; ++i, a = b) { |
|
b = points[i]; |
|
a.v = distance(b[1] - a[1], a[3], a[2], b[3], b[2], b[0] - a[0]); |
|
a.point = [0, 0]; |
|
} |
|
|
|
var beta0 = angle(points[0].v[0], points[2].v[0], points[1].v[0]), |
|
beta1 = angle(points[0].v[0], points[1].v[0], points[2].v[0]), |
|
beta2 = pi - beta0; |
|
|
|
points[2].point[1] = 0; |
|
points[0].point[0] = -(points[1].point[0] = points[0].v[0] / 2); |
|
|
|
var mean = [ |
|
points[2].point[0] = points[0].point[0] + points[2].v[0] * cos(beta0), |
|
2 * (points[0].point[1] = points[1].point[1] = points[2].v[0] * sin(beta0)) |
|
]; |
|
|
|
function forward(lambda, phi) { |
|
var sinPhi = sin(phi), |
|
cosPhi = cos(phi), |
|
v = new Array(3), i; |
|
|
|
// Compute distance and azimuth from control points. |
|
for (i = 0; i < 3; ++i) { |
|
var p = points[i]; |
|
v[i] = distance(phi - p[1], p[3], p[2], cosPhi, sinPhi, lambda - p[0]); |
|
if (!v[i][0]) return p.point; |
|
v[i][1] = longitude(v[i][1] - p.v[1]); |
|
} |
|
|
|
// Arithmetic mean of interception points. |
|
var point = mean.slice(); |
|
for (i = 0; i < 3; ++i) { |
|
var j = i == 2 ? 0 : i + 1; |
|
var a = angle(points[i].v[0], v[i][0], v[j][0]); |
|
if (v[i][1] < 0) a = -a; |
|
|
|
if (!i) { |
|
point[0] += v[i][0] * cos(a); |
|
point[1] -= v[i][0] * sin(a); |
|
} else if (i == 1) { |
|
a = beta1 - a; |
|
point[0] -= v[i][0] * cos(a); |
|
point[1] -= v[i][0] * sin(a); |
|
} else { |
|
a = beta2 - a; |
|
point[0] += v[i][0] * cos(a); |
|
point[1] += v[i][0] * sin(a); |
|
} |
|
} |
|
|
|
point[0] /= 3, point[1] /= 3; |
|
return point; |
|
} |
|
|
|
return forward; |
|
} |
|
|
|
function pointRadians(p) { |
|
return p[0] *= radians, p[1] *= radians, p; |
|
} |
|
|
|
function chamberlinAfrica() { |
|
return chamberlin([0, 22], [45, 22], [22.5, -22]) |
|
.scale(380) |
|
.center([22.5, 2]); |
|
} |
|
|
|
function chamberlin(p0, p1, p2) { // TODO order matters! |
|
var c = d3Geo.geoCentroid({type: "MultiPoint", coordinates: [p0, p1, p2]}), |
|
R = [-c[0], -c[1]], |
|
r = d3Geo.geoRotation(R), |
|
p = d3Geo.geoProjection(chamberlinRaw(pointRadians(r(p0)), pointRadians(r(p1)), pointRadians(r(p2)))).rotate(R), |
|
center = p.center; |
|
|
|
delete p.rotate; |
|
|
|
p.center = function(_) { |
|
return arguments.length ? center(r(_)) : r.invert(center()); |
|
}; |
|
|
|
return p |
|
.clipAngle(90); |
|
} |
|
|
|
function collignonRaw(lambda, phi) { |
|
var alpha = sqrt(1 - sin(phi)); |
|
return [(2 / sqrtPi) * lambda * alpha, sqrtPi * (1 - alpha)]; |
|
} |
|
|
|
collignonRaw.invert = function(x, y) { |
|
var lambda = (lambda = y / sqrtPi - 1) * lambda; |
|
return [lambda > 0 ? x * sqrt(pi / lambda) / 2 : 0, asin(1 - lambda)]; |
|
}; |
|
|
|
var collignon = function() { |
|
return d3Geo.geoProjection(collignonRaw) |
|
.scale(95.6464) |
|
.center([0, 30]); |
|
}; |
|
|
|
function craigRaw(phi0) { |
|
var tanPhi0 = tan(phi0); |
|
|
|
function forward(lambda, phi) { |
|
return [lambda, (lambda ? lambda / sin(lambda) : 1) * (sin(phi) * cos(lambda) - tanPhi0 * cos(phi))]; |
|
} |
|
|
|
forward.invert = tanPhi0 ? function(x, y) { |
|
if (x) y *= sin(x) / x; |
|
var cosLambda = cos(x); |
|
return [x, 2 * atan2(sqrt(cosLambda * cosLambda + tanPhi0 * tanPhi0 - y * y) - cosLambda, tanPhi0 - y)]; |
|
} : function(x, y) { |
|
return [x, asin(x ? y * tan(x) / x : y)]; |
|
}; |
|
|
|
return forward; |
|
} |
|
|
|
var craig = function() { |
|
return parallel1(craigRaw) |
|
.scale(249.828) |
|
.clipAngle(90); |
|
}; |
|
|
|
var sqrt3 = sqrt(3); |
|
|
|
function crasterRaw(lambda, phi) { |
|
return [sqrt3 * lambda * (2 * cos(2 * phi / 3) - 1) / sqrtPi, sqrt3 * sqrtPi * sin(phi / 3)]; |
|
} |
|
|
|
crasterRaw.invert = function(x, y) { |
|
var phi = 3 * asin(y / (sqrt3 * sqrtPi)); |
|
return [sqrtPi * x / (sqrt3 * (2 * cos(2 * phi / 3) - 1)), phi]; |
|
}; |
|
|
|
var craster = function() { |
|
return d3Geo.geoProjection(crasterRaw) |
|
.scale(156.19); |
|
}; |
|
|
|
function cylindricalEqualAreaRaw(phi0) { |
|
var cosPhi0 = cos(phi0); |
|
|
|
function forward(lambda, phi) { |
|
return [lambda * cosPhi0, sin(phi) / cosPhi0]; |
|
} |
|
|
|
forward.invert = function(x, y) { |
|
return [x / cosPhi0, asin(y * cosPhi0)]; |
|
}; |
|
|
|
return forward; |
|
} |
|
|
|
var cylindricalEqualArea = function() { |
|
return parallel1(cylindricalEqualAreaRaw) |
|
.parallel(38.58) // acos(sqrt(width / height / pi)) * radians |
|
.scale(195.044); // width / (sqrt(width / height / pi) * 2 * pi) |
|
}; |
|
|
|
function cylindricalStereographicRaw(phi0) { |
|
var cosPhi0 = cos(phi0); |
|
|
|
function forward(lambda, phi) { |
|
return [lambda * cosPhi0, (1 + cosPhi0) * tan(phi / 2)]; |
|
} |
|
|
|
forward.invert = function(x, y) { |
|
return [x / cosPhi0, atan(y / (1 + cosPhi0)) * 2]; |
|
}; |
|
|
|
return forward; |
|
} |
|
|
|
var cylindricalStereographic = function() { |
|
return parallel1(cylindricalStereographicRaw) |
|
.scale(124.75); |
|
}; |
|
|
|
function eckert1Raw(lambda, phi) { |
|
var alpha = sqrt(8 / (3 * pi)); |
|
return [ |
|
alpha * lambda * (1 - abs(phi) / pi), |
|
alpha * phi |
|
]; |
|
} |
|
|
|
eckert1Raw.invert = function(x, y) { |
|
var alpha = sqrt(8 / (3 * pi)), |
|
phi = y / alpha; |
|
return [ |
|
x / (alpha * (1 - abs(phi) / pi)), |
|
phi |
|
]; |
|
}; |
|
|
|
var eckert1 = function() { |
|
return d3Geo.geoProjection(eckert1Raw) |
|
.scale(165.664); |
|
}; |
|
|
|
function eckert2Raw(lambda, phi) { |
|
var alpha = sqrt(4 - 3 * sin(abs(phi))); |
|
return [ |
|
2 / sqrt(6 * pi) * lambda * alpha, |
|
sign(phi) * sqrt(2 * pi / 3) * (2 - alpha) |
|
]; |
|
} |
|
|
|
eckert2Raw.invert = function(x, y) { |
|
var alpha = 2 - abs(y) / sqrt(2 * pi / 3); |
|
return [ |
|
x * sqrt(6 * pi) / (2 * alpha), |
|
sign(y) * asin((4 - alpha * alpha) / 3) |
|
]; |
|
}; |
|
|
|
var eckert2 = function() { |
|
return d3Geo.geoProjection(eckert2Raw) |
|
.scale(165.664); |
|
}; |
|
|
|
function eckert3Raw(lambda, phi) { |
|
var k = sqrt(pi * (4 + pi)); |
|
return [ |
|
2 / k * lambda * (1 + sqrt(1 - 4 * phi * phi / (pi * pi))), |
|
4 / k * phi |
|
]; |
|
} |
|
|
|
eckert3Raw.invert = function(x, y) { |
|
var k = sqrt(pi * (4 + pi)) / 2; |
|
return [ |
|
x * k / (1 + sqrt(1 - y * y * (4 + pi) / (4 * pi))), |
|
y * k / 2 |
|
]; |
|
}; |
|
|
|
var eckert3 = function() { |
|
return d3Geo.geoProjection(eckert3Raw) |
|
.scale(180.739); |
|
}; |
|
|
|
function eckert4Raw(lambda, phi) { |
|
var k = (2 + halfPi) * sin(phi); |
|
phi /= 2; |
|
for (var i = 0, delta = Infinity; i < 10 && abs(delta) > epsilon; i++) { |
|
var cosPhi = cos(phi); |
|
phi -= delta = (phi + sin(phi) * (cosPhi + 2) - k) / (2 * cosPhi * (1 + cosPhi)); |
|
} |
|
return [ |
|
2 / sqrt(pi * (4 + pi)) * lambda * (1 + cos(phi)), |
|
2 * sqrt(pi / (4 + pi)) * sin(phi) |
|
]; |
|
} |
|
|
|
eckert4Raw.invert = function(x, y) { |
|
var A = y * sqrt((4 + pi) / pi) / 2, |
|
k = asin(A), |
|
c = cos(k); |
|
return [ |
|
x / (2 / sqrt(pi * (4 + pi)) * (1 + c)), |
|
asin((k + A * (c + 2)) / (2 + halfPi)) |
|
]; |
|
}; |
|
|
|
var eckert4 = function() { |
|
return d3Geo.geoProjection(eckert4Raw) |
|
.scale(180.739); |
|
}; |
|
|
|
function eckert5Raw(lambda, phi) { |
|
return [ |
|
lambda * (1 + cos(phi)) / sqrt(2 + pi), |
|
2 * phi / sqrt(2 + pi) |
|
]; |
|
} |
|
|
|
eckert5Raw.invert = function(x, y) { |
|
var k = sqrt(2 + pi), |
|
phi = y * k / 2; |
|
return [ |
|
k * x / (1 + cos(phi)), |
|
phi |
|
]; |
|
}; |
|
|
|
var eckert5 = function() { |
|
return d3Geo.geoProjection(eckert5Raw) |
|
.scale(173.044); |
|
}; |
|
|
|
function eckert6Raw(lambda, phi) { |
|
var k = (1 + halfPi) * sin(phi); |
|
for (var i = 0, delta = Infinity; i < 10 && abs(delta) > epsilon; i++) { |
|
phi -= delta = (phi + sin(phi) - k) / (1 + cos(phi)); |
|
} |
|
k = sqrt(2 + pi); |
|
return [ |
|
lambda * (1 + cos(phi)) / k, |
|
2 * phi / k |
|
]; |
|
} |
|
|
|
eckert6Raw.invert = function(x, y) { |
|
var j = 1 + halfPi, |
|
k = sqrt(j / 2); |
|
return [ |
|
x * 2 * k / (1 + cos(y *= k)), |
|
asin((y + sin(y)) / j) |
|
]; |
|
}; |
|
|
|
var eckert6 = function() { |
|
return d3Geo.geoProjection(eckert6Raw) |
|
.scale(173.044); |
|
}; |
|
|
|
var eisenlohrK = 3 + 2 * sqrt2; |
|
|
|
function eisenlohrRaw(lambda, phi) { |
|
var s0 = sin(lambda /= 2), |
|
c0 = cos(lambda), |
|
k = sqrt(cos(phi)), |
|
c1 = cos(phi /= 2), |
|
t = sin(phi) / (c1 + sqrt2 * c0 * k), |
|
c = sqrt(2 / (1 + t * t)), |
|
v = sqrt((sqrt2 * c1 + (c0 + s0) * k) / (sqrt2 * c1 + (c0 - s0) * k)); |
|
return [ |
|
eisenlohrK * (c * (v - 1 / v) - 2 * log(v)), |
|
eisenlohrK * (c * t * (v + 1 / v) - 2 * atan(t)) |
|
]; |
|
} |
|
|
|
eisenlohrRaw.invert = function(x, y) { |
|
if (!(p = augustRaw.invert(x / 1.2, y * 1.065))) return null; |
|
var lambda = p[0], phi = p[1], i = 20, p; |
|
x /= eisenlohrK, y /= eisenlohrK; |
|
do { |
|
var _0 = lambda / 2, |
|
_1 = phi / 2, |
|
s0 = sin(_0), |
|
c0 = cos(_0), |
|
s1 = sin(_1), |
|
c1 = cos(_1), |
|
cos1 = cos(phi), |
|
k = sqrt(cos1), |
|
t = s1 / (c1 + sqrt2 * c0 * k), |
|
t2 = t * t, |
|
c = sqrt(2 / (1 + t2)), |
|
v0 = (sqrt2 * c1 + (c0 + s0) * k), |
|
v1 = (sqrt2 * c1 + (c0 - s0) * k), |
|
v2 = v0 / v1, |
|
v = sqrt(v2), |
|
vm1v = v - 1 / v, |
|
vp1v = v + 1 / v, |
|
fx = c * vm1v - 2 * log(v) - x, |
|
fy = c * t * vp1v - 2 * atan(t) - y, |
|
deltatDeltaLambda = s1 && sqrt1_2 * k * s0 * t2 / s1, |
|
deltatDeltaPhi = (sqrt2 * c0 * c1 + k) / (2 * (c1 + sqrt2 * c0 * k) * (c1 + sqrt2 * c0 * k) * k), |
|
deltacDeltat = -0.5 * t * c * c * c, |
|
deltacDeltaLambda = deltacDeltat * deltatDeltaLambda, |
|
deltacDeltaPhi = deltacDeltat * deltatDeltaPhi, |
|
A = (A = 2 * c1 + sqrt2 * k * (c0 - s0)) * A * v, |
|
deltavDeltaLambda = (sqrt2 * c0 * c1 * k + cos1) / A, |
|
deltavDeltaPhi = -(sqrt2 * s0 * s1) / (k * A), |
|
deltaxDeltaLambda = vm1v * deltacDeltaLambda - 2 * deltavDeltaLambda / v + c * (deltavDeltaLambda + deltavDeltaLambda / v2), |
|
deltaxDeltaPhi = vm1v * deltacDeltaPhi - 2 * deltavDeltaPhi / v + c * (deltavDeltaPhi + deltavDeltaPhi / v2), |
|
deltayDeltaLambda = t * vp1v * deltacDeltaLambda - 2 * deltatDeltaLambda / (1 + t2) + c * vp1v * deltatDeltaLambda + c * t * (deltavDeltaLambda - deltavDeltaLambda / v2), |
|
deltayDeltaPhi = t * vp1v * deltacDeltaPhi - 2 * deltatDeltaPhi / (1 + t2) + c * vp1v * deltatDeltaPhi + c * t * (deltavDeltaPhi - deltavDeltaPhi / v2), |
|
denominator = deltaxDeltaPhi * deltayDeltaLambda - deltayDeltaPhi * deltaxDeltaLambda; |
|
if (!denominator) break; |
|
var deltaLambda = (fy * deltaxDeltaPhi - fx * deltayDeltaPhi) / denominator, |
|
deltaPhi = (fx * deltayDeltaLambda - fy * deltaxDeltaLambda) / denominator; |
|
lambda -= deltaLambda; |
|
phi = max(-halfPi, min(halfPi, phi - deltaPhi)); |
|
} while ((abs(deltaLambda) > epsilon || abs(deltaPhi) > epsilon) && --i > 0); |
|
return abs(abs(phi) - halfPi) < epsilon ? [0, phi] : i && [lambda, phi]; |
|
}; |
|
|
|
var eisenlohr = function() { |
|
return d3Geo.geoProjection(eisenlohrRaw) |
|
.scale(62.5271); |
|
}; |
|
|
|
var faheyK = cos(35 * radians); |
|
|
|
function faheyRaw(lambda, phi) { |
|
var t = tan(phi / 2); |
|
return [lambda * faheyK * sqrt(1 - t * t), (1 + faheyK) * t]; |
|
} |
|
|
|
faheyRaw.invert = function(x, y) { |
|
var t = y / (1 + faheyK); |
|
return [x && x / (faheyK * sqrt(1 - t * t)), 2 * atan(t)]; |
|
}; |
|
|
|
var fahey = function() { |
|
return d3Geo.geoProjection(faheyRaw) |
|
.scale(137.152); |
|
}; |
|
|
|
function foucautRaw(lambda, phi) { |
|
var k = phi / 2, cosk = cos(k); |
|
return [ 2 * lambda / sqrtPi * cos(phi) * cosk * cosk, sqrtPi * tan(k)]; |
|
} |
|
|
|
foucautRaw.invert = function(x, y) { |
|
var k = atan(y / sqrtPi), cosk = cos(k), phi = 2 * k; |
|
return [x * sqrtPi / 2 / (cos(phi) * cosk * cosk), phi]; |
|
}; |
|
|
|
var foucaut = function() { |
|
return d3Geo.geoProjection(foucautRaw) |
|
.scale(135.264); |
|
}; |
|
|
|
function gilbertForward(point) { |
|
return [point[0] / 2, asin(tan(point[1] / 2 * radians)) * degrees]; |
|
} |
|
|
|
function gilbertInvert(point) { |
|
return [point[0] * 2, 2 * atan(sin(point[1] * radians)) * degrees]; |
|
} |
|
|
|
var gilbert = function(projectionType) { |
|
if (projectionType == null) projectionType = d3Geo.geoOrthographic; |
|
var projection = projectionType(), |
|
equirectangular = d3Geo.geoEquirectangular().scale(degrees).precision(0).clipAngle(null).translate([0, 0]); // antimeridian cutting |
|
|
|
function gilbert(point) { |
|
return projection(gilbertForward(point)); |
|
} |
|
|
|
if (projection.invert) gilbert.invert = function(point) { |
|
return gilbertInvert(projection.invert(point)); |
|
}; |
|
|
|
gilbert.stream = function(stream) { |
|
var s1 = projection.stream(stream), s0 = equirectangular.stream({ |
|
point: function(lambda, phi) { s1.point(lambda / 2, asin(tan(-phi / 2 * radians)) * degrees); }, |
|
lineStart: function() { s1.lineStart(); }, |
|
lineEnd: function() { s1.lineEnd(); }, |
|
polygonStart: function() { s1.polygonStart(); }, |
|
polygonEnd: function() { s1.polygonEnd(); } |
|
}); |
|
s0.sphere = s1.sphere; |
|
return s0; |
|
}; |
|
|
|
function property(name) { |
|
gilbert[name] = function(_) { |
|
return arguments.length ? (projection[name](_), gilbert) : projection[name](); |
|
}; |
|
} |
|
|
|
gilbert.rotate = function(_) { |
|
return arguments.length ? (equirectangular.rotate(_), gilbert) : equirectangular.rotate(); |
|
}; |
|
|
|
gilbert.center = function(_) { |
|
return arguments.length ? (projection.center(gilbertForward(_)), gilbert) : gilbertInvert(projection.center()); |
|
}; |
|
|
|
property("clipAngle"); |
|
property("clipExtent"); |
|
property("scale"); |
|
property("translate"); |
|
property("precision"); |
|
|
|
return gilbert |
|
.scale(249.5); |
|
}; |
|
|
|
function gingeryRaw(rho, n) { |
|
var k = 2 * pi / n, |
|
rho2 = rho * rho; |
|
|
|
function forward(lambda, phi) { |
|
var p = d3Geo.geoAzimuthalEquidistantRaw(lambda, phi), |
|
x = p[0], |
|
y = p[1], |
|
r2 = x * x + y * y; |
|
|
|
if (r2 > rho2) { |
|
var r = sqrt(r2), |
|
theta = atan2(y, x), |
|
theta0 = k * round(theta / k), |
|
alpha = theta - theta0, |
|
rhoCosAlpha = rho * cos(alpha), |
|
k_ = (rho * sin(alpha) - alpha * sin(rhoCosAlpha)) / (halfPi - rhoCosAlpha), |
|
s_ = gingeryLength(alpha, k_), |
|
e = (pi - rho) / gingeryIntegrate(s_, rhoCosAlpha, pi); |
|
|
|
x = r; |
|
var i = 50, delta; |
|
do { |
|
x -= delta = (rho + gingeryIntegrate(s_, rhoCosAlpha, x) * e - r) / (s_(x) * e); |
|
} while (abs(delta) > epsilon && --i > 0); |
|
|
|
y = alpha * sin(x); |
|
if (x < halfPi) y -= k_ * (x - halfPi); |
|
|
|
var s = sin(theta0), |
|
c = cos(theta0); |
|
p[0] = x * c - y * s; |
|
p[1] = x * s + y * c; |
|
} |
|
return p; |
|
} |
|
|
|
forward.invert = function(x, y) { |
|
var r2 = x * x + y * y; |
|
if (r2 > rho2) { |
|
var r = sqrt(r2), |
|
theta = atan2(y, x), |
|
theta0 = k * round(theta / k), |
|
dTheta = theta - theta0; |
|
|
|
x = r * cos(dTheta); |
|
y = r * sin(dTheta); |
|
|
|
var x_halfPi = x - halfPi, |
|
sinx = sin(x), |
|
alpha = y / sinx, |
|
delta = x < halfPi ? Infinity : 0, |
|
i = 10; |
|
|
|
while (true) { |
|
var rhosinAlpha = rho * sin(alpha), |
|
rhoCosAlpha = rho * cos(alpha), |
|
sinRhoCosAlpha = sin(rhoCosAlpha), |
|
halfPi_RhoCosAlpha = halfPi - rhoCosAlpha, |
|
k_ = (rhosinAlpha - alpha * sinRhoCosAlpha) / halfPi_RhoCosAlpha, |
|
s_ = gingeryLength(alpha, k_); |
|
|
|
if (abs(delta) < epsilon2 || !--i) break; |
|
|
|
alpha -= delta = (alpha * sinx - k_ * x_halfPi - y) / ( |
|
sinx - x_halfPi * 2 * ( |
|
halfPi_RhoCosAlpha * (rhoCosAlpha + alpha * rhosinAlpha * cos(rhoCosAlpha) - sinRhoCosAlpha) - |
|
rhosinAlpha * (rhosinAlpha - alpha * sinRhoCosAlpha) |
|
) / (halfPi_RhoCosAlpha * halfPi_RhoCosAlpha)); |
|
} |
|
r = rho + gingeryIntegrate(s_, rhoCosAlpha, x) * (pi - rho) / gingeryIntegrate(s_, rhoCosAlpha, pi); |
|
theta = theta0 + alpha; |
|
x = r * cos(theta); |
|
y = r * sin(theta); |
|
} |
|
return d3Geo.geoAzimuthalEquidistantRaw.invert(x, y); |
|
}; |
|
|
|
return forward; |
|
} |
|
|
|
function gingeryLength(alpha, k) { |
|
return function(x) { |
|
var y_ = alpha * cos(x); |
|
if (x < halfPi) y_ -= k; |
|
return sqrt(1 + y_ * y_); |
|
}; |
|
} |
|
|
|
// Numerical integration: trapezoidal rule. |
|
function gingeryIntegrate(f, a, b) { |
|
var n = 50, |
|
h = (b - a) / n, |
|
s = f(a) + f(b); |
|
for (var i = 1, x = a; i < n; ++i) s += 2 * f(x += h); |
|
return s * 0.5 * h; |
|
} |
|
|
|
var gingery = function() { |
|
var n = 6, |
|
rho = 30 * radians, |
|
cRho = cos(rho), |
|
sRho = sin(rho), |
|
m = d3Geo.geoProjectionMutator(gingeryRaw), |
|
p = m(rho, n), |
|
stream_ = p.stream, |
|
epsilon$$1 = 1e-2, |
|
cr = -cos(epsilon$$1 * radians), |
|
sr = sin(epsilon$$1 * radians); |
|
|
|
p.radius = function(_) { |
|
if (!arguments.length) return rho * degrees; |
|
cRho = cos(rho = _ * radians); |
|
sRho = sin(rho); |
|
return m(rho, n); |
|
}; |
|
|
|
p.lobes = function(_) { |
|
if (!arguments.length) return n; |
|
return m(rho, n = +_); |
|
}; |
|
|
|
p.stream = function(stream) { |
|
var rotate = p.rotate(), |
|
rotateStream = stream_(stream), |
|
sphereStream = (p.rotate([0, 0]), stream_(stream)); |
|
p.rotate(rotate); |
|
rotateStream.sphere = function() { |
|
sphereStream.polygonStart(), sphereStream.lineStart(); |
|
for (var i = 0, delta = 2 * pi / n, phi = 0; i < n; ++i, phi -= delta) { |
|
sphereStream.point(atan2(sr * cos(phi), cr) * degrees, asin(sr * sin(phi)) * degrees); |
|
sphereStream.point(atan2(sRho * cos(phi - delta / 2), cRho) * degrees, asin(sRho * sin(phi - delta / 2)) * degrees); |
|
} |
|
sphereStream.lineEnd(), sphereStream.polygonEnd(); |
|
}; |
|
return rotateStream; |
|
}; |
|
|
|
return p |
|
.rotate([90, -40]) |
|
.scale(91.7095) |
|
.clipAngle(180 - 1e-3); |
|
}; |
|
|
|
var ginzburgPolyconicRaw = function(a, b, c, d, e, f, g, h) { |
|
if (arguments.length < 8) h = 0; |
|
|
|
function forward(lambda, phi) { |
|
if (!phi) return [a * lambda / pi, 0]; |
|
var phi2 = phi * phi, |
|
xB = a + phi2 * (b + phi2 * (c + phi2 * d)), |
|
yB = phi * (e - 1 + phi2 * (f - h + phi2 * g)), |
|
m = (xB * xB + yB * yB) / (2 * yB), |
|
alpha = lambda * asin(xB / m) / pi; |
|
return [m * sin(alpha), phi * (1 + phi2 * h) + m * (1 - cos(alpha))]; |
|
} |
|
|
|
forward.invert = function(x, y) { |
|
var lambda = pi * x / a, |
|
phi = y, |
|
deltaLambda, deltaPhi, i = 50; |
|
do { |
|
var phi2 = phi * phi, |
|
xB = a + phi2 * (b + phi2 * (c + phi2 * d)), |
|
yB = phi * (e - 1 + phi2 * (f - h + phi2 * g)), |
|
p = xB * xB + yB * yB, |
|
q = 2 * yB, |
|
m = p / q, |
|
m2 = m * m, |
|
dAlphadLambda = asin(xB / m) / pi, |
|
alpha = lambda * dAlphadLambda, |
|
xB2 = xB * xB, |
|
dxBdPhi = (2 * b + phi2 * (4 * c + phi2 * 6 * d)) * phi, |
|
dyBdPhi = e + phi2 * (3 * f + phi2 * 5 * g), |
|
dpdPhi = 2 * (xB * dxBdPhi + yB * (dyBdPhi - 1)), |
|
dqdPhi = 2 * (dyBdPhi - 1), |
|
dmdPhi = (dpdPhi * q - p * dqdPhi) / (q * q), |
|
cosAlpha = cos(alpha), |
|
sinAlpha = sin(alpha), |
|
mcosAlpha = m * cosAlpha, |
|
msinAlpha = m * sinAlpha, |
|
dAlphadPhi = ((lambda / pi) * (1 / sqrt(1 - xB2 / m2)) * (dxBdPhi * m - xB * dmdPhi)) / m2, |
|
fx = msinAlpha - x, |
|
fy = phi * (1 + phi2 * h) + m - mcosAlpha - y, |
|
deltaxDeltaPhi = dmdPhi * sinAlpha + mcosAlpha * dAlphadPhi, |
|
deltaxDeltaLambda = mcosAlpha * dAlphadLambda, |
|
deltayDeltaPhi = 1 + dmdPhi - (dmdPhi * cosAlpha - msinAlpha * dAlphadPhi), |
|
deltayDeltaLambda = msinAlpha * dAlphadLambda, |
|
denominator = deltaxDeltaPhi * deltayDeltaLambda - deltayDeltaPhi * deltaxDeltaLambda; |
|
if (!denominator) break; |
|
lambda -= deltaLambda = (fy * deltaxDeltaPhi - fx * deltayDeltaPhi) / denominator; |
|
phi -= deltaPhi = (fx * deltayDeltaLambda - fy * deltaxDeltaLambda) / denominator; |
|
} while ((abs(deltaLambda) > epsilon || abs(deltaPhi) > epsilon) && --i > 0); |
|
return [lambda, phi]; |
|
}; |
|
|
|
return forward; |
|
}; |
|
|
|
var ginzburg4Raw = ginzburgPolyconicRaw(2.8284, -1.6988, 0.75432, -0.18071, 1.76003, -0.38914, 0.042555); |
|
|
|
var ginzburg4 = function() { |
|
return d3Geo.geoProjection(ginzburg4Raw) |
|
.scale(149.995); |
|
}; |
|
|
|
var ginzburg5Raw = ginzburgPolyconicRaw(2.583819, -0.835827, 0.170354, -0.038094, 1.543313, -0.411435,0.082742); |
|
|
|
var ginzburg5 = function() { |
|
return d3Geo.geoProjection(ginzburg5Raw) |
|
.scale(153.93); |
|
}; |
|
|
|
var ginzburg6Raw = ginzburgPolyconicRaw(5 / 6 * pi, -0.62636, -0.0344, 0, 1.3493, -0.05524, 0, 0.045); |
|
|
|
var ginzburg6 = function() { |
|
return d3Geo.geoProjection(ginzburg6Raw) |
|
.scale(130.945); |
|
}; |
|
|
|
function ginzburg8Raw(lambda, phi) { |
|
var lambda2 = lambda * lambda, |
|
phi2 = phi * phi; |
|
return [ |
|
lambda * (1 - 0.162388 * phi2) * (0.87 - 0.000952426 * lambda2 * lambda2), |
|
phi * (1 + phi2 / 12) |
|
]; |
|
} |
|
|
|
ginzburg8Raw.invert = function(x, y) { |
|
var lambda = x, |
|
phi = y, |
|
i = 50, delta; |
|
do { |
|
var phi2 = phi * phi; |
|
phi -= delta = (phi * (1 + phi2 / 12) - y) / (1 + phi2 / 4); |
|
} while (abs(delta) > epsilon && --i > 0); |
|
i = 50; |
|
x /= 1 -0.162388 * phi2; |
|
do { |
|
var lambda4 = (lambda4 = lambda * lambda) * lambda4; |
|
lambda -= delta = (lambda * (0.87 - 0.000952426 * lambda4) - x) / (0.87 - 0.00476213 * lambda4); |
|
} while (abs(delta) > epsilon && --i > 0); |
|
return [lambda, phi]; |
|
}; |
|
|
|
var ginzburg8 = function() { |
|
return d3Geo.geoProjection(ginzburg8Raw) |
|
.scale(131.747); |
|
}; |
|
|
|
var ginzburg9Raw = ginzburgPolyconicRaw(2.6516, -0.76534, 0.19123, -0.047094, 1.36289, -0.13965,0.031762); |
|
|
|
var ginzburg9 = function() { |
|
return d3Geo.geoProjection(ginzburg9Raw) |
|
.scale(131.087); |
|
}; |
|
|
|
var squareRaw = function(project) { |
|
var dx = project(halfPi, 0)[0] - project(-halfPi, 0)[0]; |
|
|
|
function projectSquare(lambda, phi) { |
|
var s = lambda > 0 ? -0.5 : 0.5, |
|
point = project(lambda + s * pi, phi); |
|
point[0] -= s * dx; |
|
return point; |
|
} |
|
|
|
if (project.invert) projectSquare.invert = function(x, y) { |
|
var s = x > 0 ? -0.5 : 0.5, |
|
location = project.invert(x + s * dx, y), |
|
lambda = location[0] - s * pi; |
|
if (lambda < -pi) lambda += 2 * pi; |
|
else if (lambda > pi) lambda -= 2 * pi; |
|
location[0] = lambda; |
|
return location; |
|
}; |
|
|
|
return projectSquare; |
|
}; |
|
|
|
function gringortenRaw(lambda, phi) { |
|
var sLambda = sign(lambda), |
|
sPhi = sign(phi), |
|
cosPhi = cos(phi), |
|
x = cos(lambda) * cosPhi, |
|
y = sin(lambda) * cosPhi, |
|
z = sin(sPhi * phi); |
|
lambda = abs(atan2(y, z)); |
|
phi = asin(x); |
|
if (abs(lambda - halfPi) > epsilon) lambda %= halfPi; |
|
var point = gringortenHexadecant(lambda > pi / 4 ? halfPi - lambda : lambda, phi); |
|
if (lambda > pi / 4) z = point[0], point[0] = -point[1], point[1] = -z; |
|
return (point[0] *= sLambda, point[1] *= -sPhi, point); |
|
} |
|
|
|
gringortenRaw.invert = function(x, y) { |
|
if (abs(x) > 1) x = sign(x) * 2 - x; |
|
if (abs(y) > 1) y = sign(y) * 2 - y; |
|
var sx = sign(x), |
|
sy = sign(y), |
|
x0 = -sx * x, |
|
y0 = -sy * y, |
|
t = y0 / x0 < 1, |
|
p = gringortenHexadecantInvert(t ? y0 : x0, t ? x0 : y0), |
|
lambda = p[0], |
|
phi = p[1], |
|
cosPhi = cos(phi); |
|
if (t) lambda = -halfPi - lambda; |
|
return [sx * (atan2(sin(lambda) * cosPhi, -sin(phi)) + pi), sy * asin(cos(lambda) * cosPhi)]; |
|
}; |
|
|
|
function gringortenHexadecant(lambda, phi) { |
|
if (phi === halfPi) return [0, 0]; |
|
|
|
var sinPhi = sin(phi), |
|
r = sinPhi * sinPhi, |
|
r2 = r * r, |
|
j = 1 + r2, |
|
k = 1 + 3 * r2, |
|
q = 1 - r2, |
|
z = asin(1 / sqrt(j)), |
|
v = q + r * j * z, |
|
p2 = (1 - sinPhi) / v, |
|
p = sqrt(p2), |
|
a2 = p2 * j, |
|
a = sqrt(a2), |
|
h = p * q, |
|
x, |
|
i; |
|
|
|
if (lambda === 0) return [0, -(h + r * a)]; |
|
|
|
var cosPhi = cos(phi), |
|
secPhi = 1 / cosPhi, |
|
drdPhi = 2 * sinPhi * cosPhi, |
|
dvdPhi = (-3 * r + z * k) * drdPhi, |
|
dp2dPhi = (-v * cosPhi - (1 - sinPhi) * dvdPhi) / (v * v), |
|
dpdPhi = (0.5 * dp2dPhi) / p, |
|
dhdPhi = q * dpdPhi - 2 * r * p * drdPhi, |
|
dra2dPhi = r * j * dp2dPhi + p2 * k * drdPhi, |
|
mu = -secPhi * drdPhi, |
|
nu = -secPhi * dra2dPhi, |
|
zeta = -2 * secPhi * dhdPhi, |
|
lambda1 = 4 * lambda / pi, |
|
delta; |
|
|
|
// Slower but accurate bisection method. |
|
if (lambda > 0.222 * pi || phi < pi / 4 && lambda > 0.175 * pi) { |
|
x = (h + r * sqrt(a2 * (1 + r2) - h * h)) / (1 + r2); |
|
if (lambda > pi / 4) return [x, x]; |
|
var x1 = x, x0 = 0.5 * x; |
|
x = 0.5 * (x0 + x1), i = 50; |
|
do { |
|
var g = sqrt(a2 - x * x), |
|
f = (x * (zeta + mu * g) + nu * asin(x / a)) - lambda1; |
|
if (!f) break; |
|
if (f < 0) x0 = x; |
|
else x1 = x; |
|
x = 0.5 * (x0 + x1); |
|
} while (abs(x1 - x0) > epsilon && --i > 0); |
|
} |
|
|
|
// Newton-Raphson. |
|
else { |
|
x = epsilon, i = 25; |
|
do { |
|
var x2 = x * x, |
|
g2 = sqrt(a2 - x2), |
|
zetaMug = zeta + mu * g2, |
|
f2 = x * zetaMug + nu * asin(x / a) - lambda1, |
|
df = zetaMug + (nu - mu * x2) / g2; |
|
x -= delta = g2 ? f2 / df : 0; |
|
} while (abs(delta) > epsilon && --i > 0); |
|
} |
|
|
|
return [x, -h - r * sqrt(a2 - x * x)]; |
|
} |
|
|
|
function gringortenHexadecantInvert(x, y) { |
|
var x0 = 0, |
|
x1 = 1, |
|
r = 0.5, |
|
i = 50; |
|
|
|
while (true) { |
|
var r2 = r * r, |
|
sinPhi = sqrt(r), |
|
z = asin(1 / sqrt(1 + r2)), |
|
v = (1 - r2) + r * (1 + r2) * z, |
|
p2 = (1 - sinPhi) / v, |
|
p = sqrt(p2), |
|
a2 = p2 * (1 + r2), |
|
h = p * (1 - r2), |
|
g2 = a2 - x * x, |
|
g = sqrt(g2), |
|
y0 = y + h + r * g; |
|
if (abs(x1 - x0) < epsilon2 || --i === 0 || y0 === 0) break; |
|
if (y0 > 0) x0 = r; |
|
else x1 = r; |
|
r = 0.5 * (x0 + x1); |
|
} |
|
|
|
if (!i) return null; |
|
|
|
var phi = asin(sinPhi), |
|
cosPhi = cos(phi), |
|
secPhi = 1 / cosPhi, |
|
drdPhi = 2 * sinPhi * cosPhi, |
|
dvdPhi = (-3 * r + z * (1 + 3 * r2)) * drdPhi, |
|
dp2dPhi = (-v * cosPhi - (1 - sinPhi) * dvdPhi) / (v * v), |
|
dpdPhi = 0.5 * dp2dPhi / p, |
|
dhdPhi = (1 - r2) * dpdPhi - 2 * r * p * drdPhi, |
|
zeta = -2 * secPhi * dhdPhi, |
|
mu = -secPhi * drdPhi, |
|
nu = -secPhi * (r * (1 + r2) * dp2dPhi + p2 * (1 + 3 * r2) * drdPhi); |
|
|
|
return [pi / 4 * (x * (zeta + mu * g) + nu * asin(x / sqrt(a2))), phi]; |
|
} |
|
|
|
var gringorten = function() { |
|
return d3Geo.geoProjection(squareRaw(gringortenRaw)) |
|
.scale(239.75); |
|
}; |
|
|
|
// Returns [sn, cn, dn](u + iv|m). |
|
function ellipticJi(u, v, m) { |
|
var a, b, c; |
|
if (!u) { |
|
b = ellipticJ(v, 1 - m); |
|
return [ |
|
[0, b[0] / b[1]], |
|
[1 / b[1], 0], |
|
[b[2] / b[1], 0] |
|
]; |
|
} |
|
a = ellipticJ(u, m); |
|
if (!v) return [[a[0], 0], [a[1], 0], [a[2], 0]]; |
|
b = ellipticJ(v, 1 - m); |
|
c = b[1] * b[1] + m * a[0] * a[0] * b[0] * b[0]; |
|
return [ |
|
[a[0] * b[2] / c, a[1] * a[2] * b[0] * b[1] / c], |
|
[a[1] * b[1] / c, -a[0] * a[2] * b[0] * b[2] / c], |
|
[a[2] * b[1] * b[2] / c, -m * a[0] * a[1] * b[0] / c] |
|
]; |
|
} |
|
|
|
// Returns [sn, cn, dn, ph](u|m). |
|
function ellipticJ(u, m) { |
|
var ai, b, phi, t, twon; |
|
if (m < epsilon) { |
|
t = sin(u); |
|
b = cos(u); |
|
ai = m * (u - t * b) / 4; |
|
return [ |
|
t - ai * b, |
|
b + ai * t, |
|
1 - m * t * t / 2, |
|
u - ai |
|
]; |
|
} |
|
if (m >= 1 - epsilon) { |
|
ai = (1 - m) / 4; |
|
b = cosh(u); |
|
t = tanh(u); |
|
phi = 1 / b; |
|
twon = b * sinh(u); |
|
return [ |
|
t + ai * (twon - u) / (b * b), |
|
phi - ai * t * phi * (twon - u), |
|
phi + ai * t * phi * (twon + u), |
|
2 * atan(exp(u)) - halfPi + ai * (twon - u) / b |
|
]; |
|
} |
|
|
|
var a = [1, 0, 0, 0, 0, 0, 0, 0, 0], |
|
c = [sqrt(m), 0, 0, 0, 0, 0, 0, 0, 0], |
|
i = 0; |
|
b = sqrt(1 - m); |
|
twon = 1; |
|
|
|
while (abs(c[i] / a[i]) > epsilon && i < 8) { |
|
ai = a[i++]; |
|
c[i] = (ai - b) / 2; |
|
a[i] = (ai + b) / 2; |
|
b = sqrt(ai * b); |
|
twon *= 2; |
|
} |
|
|
|
phi = twon * a[i] * u; |
|
do { |
|
t = c[i] * sin(b = phi) / a[i]; |
|
phi = (asin(t) + phi) / 2; |
|
} while (--i); |
|
|
|
return [sin(phi), t = cos(phi), t / cos(phi - b), phi]; |
|
} |
|
|
|
// Calculate F(phi+iPsi|m). |
|
// See Abramowitz and Stegun, 17.4.11. |
|
function ellipticFi(phi, psi, m) { |
|
var r = abs(phi), |
|
i = abs(psi), |
|
sinhPsi = sinh(i); |
|
if (r) { |
|
var cscPhi = 1 / sin(r), |
|
cotPhi2 = 1 / (tan(r) * tan(r)), |
|
b = -(cotPhi2 + m * (sinhPsi * sinhPsi * cscPhi * cscPhi) - 1 + m), |
|
c = (m - 1) * cotPhi2, |
|
cotLambda2 = (-b + sqrt(b * b - 4 * c)) / 2; |
|
return [ |
|
ellipticF(atan(1 / sqrt(cotLambda2)), m) * sign(phi), |
|
ellipticF(atan(sqrt((cotLambda2 / cotPhi2 - 1) / m)), 1 - m) * sign(psi) |
|
]; |
|
} |
|
return [ |
|
0, |
|
ellipticF(atan(sinhPsi), 1 - m) * sign(psi) |
|
]; |
|
} |
|
|
|
// Calculate F(phi|m) where m = k² = sin²α. |
|
// See Abramowitz and Stegun, 17.6.7. |
|
function ellipticF(phi, m) { |
|
if (!m) return phi; |
|
if (m === 1) return log(tan(phi / 2 + quarterPi)); |
|
var a = 1, |
|
b = sqrt(1 - m), |
|
c = sqrt(m); |
|
for (var i = 0; abs(c) > epsilon; i++) { |
|
if (phi % pi) { |
|
var dPhi = atan(b * tan(phi) / a); |
|
if (dPhi < 0) dPhi += pi; |
|
phi += dPhi + ~~(phi / pi) * pi; |
|
} else phi += phi; |
|
c = (a + b) / 2; |
|
b = sqrt(a * b); |
|
c = ((a = c) - b) / 2; |
|
} |
|
return phi / (pow(2, i) * a); |
|
} |
|
|
|
function guyouRaw(lambda, phi) { |
|
var k_ = (sqrt2 - 1) / (sqrt2 + 1), |
|
k = sqrt(1 - k_ * k_), |
|
K = ellipticF(halfPi, k * k), |
|
f = -1, |
|
psi = log(tan(pi / 4 + abs(phi) / 2)), |
|
r = exp(f * psi) / sqrt(k_), |
|
at = guyouComplexAtan(r * cos(f * lambda), r * sin(f * lambda)), |
|
t = ellipticFi(at[0], at[1], k * k); |
|
return [-t[1], (phi >= 0 ? 1 : -1) * (0.5 * K - t[0])]; |
|
} |
|
|
|
function guyouComplexAtan(x, y) { |
|
var x2 = x * x, |
|
y_1 = y + 1, |
|
t = 1 - x2 - y * y; |
|
return [ |
|
0.5 * ((x >= 0 ? halfPi : -halfPi) - atan2(t, 2 * x)), |
|
-0.25 * log(t * t + 4 * x2) +0.5 * log(y_1 * y_1 + x2) |
|
]; |
|
} |
|
|
|
function guyouComplexDivide(a, b) { |
|
var denominator = b[0] * b[0] + b[1] * b[1]; |
|
return [ |
|
(a[0] * b[0] + a[1] * b[1]) / denominator, |
|
(a[1] * b[0] - a[0] * b[1]) / denominator |
|
]; |
|
} |
|
|
|
guyouRaw.invert = function(x, y) { |
|
var k_ = (sqrt2 - 1) / (sqrt2 + 1), |
|
k = sqrt(1 - k_ * k_), |
|
K = ellipticF(halfPi, k * k), |
|
f = -1, |
|
j = ellipticJi(0.5 * K - y, -x, k * k), |
|
tn = guyouComplexDivide(j[0], j[1]), |
|
lambda = atan2(tn[1], tn[0]) / f; |
|
return [ |
|
lambda, |
|
2 * atan(exp(0.5 / f * log(k_ * tn[0] * tn[0] + k_ * tn[1] * tn[1]))) - halfPi |
|
]; |
|
}; |
|
|
|
var guyou = function() { |
|
return d3Geo.geoProjection(squareRaw(guyouRaw)) |
|
.scale(151.496); |
|
}; |
|
|
|
function hammerRaw(A, B) { |
|
if (arguments.length < 2) B = A; |
|
if (B === 1) return d3Geo.geoAzimuthalEqualAreaRaw; |
|
if (B === Infinity) return hammerQuarticAuthalicRaw; |
|
|
|
function forward(lambda, phi) { |
|
var coordinates = d3Geo.geoAzimuthalEqualAreaRaw(lambda / B, phi); |
|
coordinates[0] *= A; |
|
return coordinates; |
|
} |
|
|
|
forward.invert = function(x, y) { |
|
var coordinates = d3Geo.geoAzimuthalEqualAreaRaw.invert(x / A, y); |
|
coordinates[0] *= B; |
|
return coordinates; |
|
}; |
|
|
|
return forward; |
|
} |
|
|
|
function hammerQuarticAuthalicRaw(lambda, phi) { |
|
return [ |
|
lambda * cos(phi) / cos(phi /= 2), |
|
2 * sin(phi) |
|
]; |
|
} |
|
|
|
hammerQuarticAuthalicRaw.invert = function(x, y) { |
|
var phi = 2 * asin(y / 2); |
|
return [ |
|
x * cos(phi / 2) / cos(phi), |
|
phi |
|
]; |
|
}; |
|
|
|
var hammer = function() { |
|
var B = 2, |
|
m = d3Geo.geoProjectionMutator(hammerRaw), |
|
p = m(B); |
|
|
|
p.coefficient = function(_) { |
|
if (!arguments.length) return B; |
|
return m(B = +_); |
|
}; |
|
|
|
return p |
|
.scale(169.529); |
|
}; |
|
|
|
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) |
|
]; |
|
}; |
|
} |
|
|
|
var hammerRetroazimuthal = function() { |
|
var phi0 = 0, |
|
m = d3Geo.geoProjectionMutator(hammerRetroazimuthalRaw), |
|
p = m(phi0), |
|
rotate_ = p.rotate, |
|
stream_ = p.stream, |
|
circle = d3Geo.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$$1 = 1e-2, |
|
ring = circle.radius(90 - epsilon$$1)().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$$1)().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); |
|
}; |
|
|
|
var healpixParallel = 41 + 48 / 36 + 37 / 3600; |
|
var healpixLambert = cylindricalEqualAreaRaw(0); |
|
|
|
function healpixRaw(H) { |
|
var phi0 = healpixParallel * radians, |
|
dx = collignonRaw(pi, phi0)[0] - collignonRaw(-pi, phi0)[0], |
|
y0 = healpixLambert(0, phi0)[1], |
|
y1 = collignonRaw(0, phi0)[1], |
|
dy1 = sqrtPi - y1, |
|
k = tau / H, |
|
w = 4 / tau, |
|
h = y0 + (dy1 * dy1 * 4) / tau; |
|
|
|
function forward(lambda, phi) { |
|
var point, |
|
phi2 = abs(phi); |
|
if (phi2 > phi0) { |
|
var i = min(H - 1, max(0, floor((lambda + pi) / k))); |
|
lambda += pi * (H - 1) / H - i * k; |
|
point = collignonRaw(lambda, phi2); |
|
point[0] = point[0] * tau / dx - tau * (H - 1) / (2 * H) + i * tau / H; |
|
point[1] = y0 + (point[1] - y1) * 4 * dy1 / tau; |
|
if (phi < 0) point[1] = -point[1]; |
|
} else { |
|
point = healpixLambert(lambda, phi); |
|
} |
|
point[0] *= w, point[1] /= h; |
|
return point; |
|
} |
|
|
|
forward.invert = function(x, y) { |
|
x /= w, y *= h; |
|
var y2 = abs(y); |
|
if (y2 > y0) { |
|
var i = min(H - 1, max(0, floor((x + pi) / k))); |
|
x = (x + pi * (H - 1) / H - i * k) * dx / tau; |
|
var point = collignonRaw.invert(x, 0.25 * (y2 - y0) * tau / dy1 + y1); |
|
point[0] -= pi * (H - 1) / H - i * k; |
|
if (y < 0) point[1] = -point[1]; |
|
return point; |
|
} |
|
return healpixLambert.invert(x, y); |
|
}; |
|
|
|
return forward; |
|
} |
|
|
|
function sphere(step) { |
|
return { |
|
type: "Polygon", |
|
coordinates: [ |
|
d3Array.range(-180, 180 + step / 2, step).map(function(x, i) { return [x, i & 1 ? 90 - 1e-6 : healpixParallel]; }) |
|
.concat(d3Array.range(180, -180 - step / 2, -step).map(function(x, i) { return [x, i & 1 ? -90 + 1e-6 : -healpixParallel]; })) |
|
] |
|
}; |
|
} |
|
|
|
var healpix = function() { |
|
var H = 4, |
|
m = d3Geo.geoProjectionMutator(healpixRaw), |
|
p = m(H), |
|
stream_ = p.stream; |
|
|
|
p.lobes = function(_) { |
|
return arguments.length ? m(H = +_) : H; |
|
}; |
|
|
|
p.stream = function(stream) { |
|
var rotate = p.rotate(), |
|
rotateStream = stream_(stream), |
|
sphereStream = (p.rotate([0, 0]), stream_(stream)); |
|
p.rotate(rotate); |
|
rotateStream.sphere = function() { d3Geo.geoStream(sphere(180 / H), sphereStream); }; |
|
return rotateStream; |
|
}; |
|
|
|
return p |
|
.scale(239.75); |
|
}; |
|
|
|
function hillRaw(K) { |
|
var L = 1 + K, |
|
sinBt = sin(1 / L), |
|
Bt = asin(sinBt), |
|
A = 2 * sqrt(pi / (B = pi + 4 * Bt * L)), |
|
B, |
|
rho0 = 0.5 * A * (L + sqrt(K * (2 + K))), |
|
K2 = K * K, |
|
L2 = L * L; |
|
|
|
function forward(lambda, phi) { |
|
var t = 1 - sin(phi), |
|
rho, |
|
omega; |
|
if (t && t < 2) { |
|
var theta = halfPi - phi, i = 25, delta; |
|
do { |
|
var sinTheta = sin(theta), |
|
cosTheta = cos(theta), |
|
Bt_Bt1 = Bt + atan2(sinTheta, L - cosTheta), |
|
C = 1 + L2 - 2 * L * cosTheta; |
|
theta -= delta = (theta - K2 * Bt - L * sinTheta + C * Bt_Bt1 -0.5 * t * B) / (2 * L * sinTheta * Bt_Bt1); |
|
} while (abs(delta) > epsilon2 && --i > 0); |
|
rho = A * sqrt(C); |
|
omega = lambda * Bt_Bt1 / pi; |
|
} else { |
|
rho = A * (K + t); |
|
omega = lambda * Bt / pi; |
|
} |
|
return [ |
|
rho * sin(omega), |
|
rho0 - rho * cos(omega) |
|
]; |
|
} |
|
|
|
forward.invert = function(x, y) { |
|
var rho2 = x * x + (y -= rho0) * y, |
|
cosTheta = (1 + L2 - rho2 / (A * A)) / (2 * L), |
|
theta = acos(cosTheta), |
|
sinTheta = sin(theta), |
|
Bt_Bt1 = Bt + atan2(sinTheta, L - cosTheta); |
|
return [ |
|
asin(x / sqrt(rho2)) * pi / Bt_Bt1, |
|
asin(1 - 2 * (theta - K2 * Bt - L * sinTheta + (1 + L2 - 2 * L * cosTheta) * Bt_Bt1) / B) |
|
]; |
|
}; |
|
|
|
return forward; |
|
} |
|
|
|
var hill = function() { |
|
var K = 1, |
|
m = d3Geo.geoProjectionMutator(hillRaw), |
|
p = m(K); |
|
|
|
|