'use strict';
if (typeof module!='undefined' && module.exports) var LatLon = require('./latlon-ellipsoidal.js'); // CommonJS (Node)
// http://geographiclib.sourceforge.net/html/js/tutorial-2-interface.html
/**
* Direct (simplistic) implementation of Algorithms for geodesics, Charles F F Karney, Journal of Geodesy 2012
* (see also Geodesics on an ellipsoid of revolution, Charles F F Karney, 2011 for further explanations).
*
* Notable historical background to methods used.
*
* @param point
* @returns {number|*}
*/
LatLon.prototype.distanceTo = function(point) {
if (!(point instanceof LatLon)) throw new TypeError('point is not LatLon object');
return this.inverse(point).distance;
};
LatLon.prototype.destinationPoint = function(distance, initialBearing) {
return this.direct(Number(distance), Number(initialBearing)).point;
};
LatLon.prototype.finalBearingOn = function(distance, initialBearing) {
return this.direct(Number(distance), Number(initialBearing)).finalBearing;
};
LatLon.prototype.direct = function(distance, initialBearing) {
const φ1 = this.lat.toRadians();
const α1 = initialBearing.toRadians();
const s = distance;
const a = this.datum.ellipsoid.a, b = this.datum.ellipsoid.b, f = this.datum.ellipsoid.f;
const n = f / (2-f); // third flattening
const e2 = f * (2-f); // eccentricity squared
const eʹ2 = e2 / (1-e2); // second eccentricity squared
const sinα1 = Math.sin(α1), cosα1 = Math.cos(α1);
const sinφ1 = Math.sin(φ1), cosφ1 = Math.cos(φ1);
// β is the reduced latitude (on the auxiliary sphere); tanβ = (1-f)·tanφ
const sinβ1 = sinφ1*(1-f) / Math.hypot(sinφ1*(1-f), cosφ1);
const cosβ1 = cosφ1 / Math.hypot(sinφ1*(1-f), cosφ1);
// α0 is the azimuth of the geodesic at the equator
const sinα0 = sinα1 * cosβ1; // Clairaut’s relation
const cosα0 = Math.hypot(cosα1, sinα1*sinβ1); // [1]
// σ is the arc length from the intersection of the geodesic with the equator
const sinσ1 = sinβ1 / Math.hypot(sinβ1, cosβ1*cosα1);
const cosσ1 = cosβ1*cosα1 / Math.hypot(sinβ1, cosβ1*cosα1);
// ω is the longitude on the auxiliary sphere
const sinω1 = sinα0 * sinβ1;
const cosω1 = cosβ1 * cosα1;
// ε is series expansion parameter
const k2 = eʹ2 * cosα0 *cosα0; // k²; k = eʹ·cosα0
const ε = k2 / (2 * (1 + Math.sqrt(1+k2)) + k2); // ε = k² / (√(1+k²) + 1)²
const ε2 = ε*ε, ε3 = ε*ε2, ε4 = ε*ε3, ε5 = ε*ε4, ε6 = ε*ε5; // powers of ε
const A1 = (1 + ε2/4 + ε4/64 + ε6/256) / (1 - ε); // A₁; eq.17
const C1 = [ null, // C₁; eq.18
-1/2*ε + 3/16*ε3 - 1/32*ε5,
-1/16*ε2 + 1/32*ε4 - 9/2048*ε6,
-1/48*ε3 + 3/256*ε5,
-5/512*ε4 + 3/512*ε6,
-7/1280*ε5,
-7/2048*ε6,
];
const ΣC1σ1 = clenshawSinSeries(sinσ1, cosσ1, C1); // Σ C₁ₗsin2lσ
const σ1 = Math.atan2(sinσ1, cosσ1);
const I1σ1 = A1 * (σ1 + ΣC1σ1); // I₁(σ): distance integral for s/b
// s₁, s₂ are distances from intersection of the geodesic with the equator
const s1 = I1σ1 * b; // from distance integral for s/b
const s2 = s1 + s;
const τ2 = s2 / (b*A1);
const sinτ2 = Math.sin(τ2), cosτ2 = Math.cos(τ2);
console.log('sinτ2', sinτ2, 'cosτ2', cosτ2)
const Cʹ1 = [ null, // Cʹ₁; eq.21
1/2*ε - 9/32*ε3 + 205/1536*ε5,
5/16*ε2 - 37/96*ε4 + 1335/4096*ε6,
29/96*ε3 - 75/128*ε5,
539/1536*ε4 - 2391/2560*ε6,
3467/7680*ε5,
38081/61440*ε6,
];
const ΣCʹ1τ2 = clenshawSinSeries(sinτ2, cosτ2, Cʹ1); // Σ Cʹ₁ₗsin2lτ2
const σ2 = τ2 + ΣCʹ1τ2;
const sinσ2 = Math.sin(σ2), cosσ2 = Math.cos(σ2);
const tanβ2 = cosα0*sinσ2 / Math.hypot(cosα0*cosσ2, sinα0); // β = Arg(|cosα0cosσ +i·sinα0| + i·cosα0sinσ)
const A3 = 1 - (1/2 - 1/2*n)*ε - (1/4 + 1/8*n - 3/8*n*n)*ε2 - (1/16 + 3/16*n + 1/16*n*n)*ε3 - (3/64 + 1/32*n)*ε4 - 3/128*ε5;
const C3 = [ null, // eq.24, 25
(1/4 - 1/4*n)*ε + (1/8 + 1/8*n*n)*ε2 + (3/64 + 3/64*n - 1/64*n*n)*ε3 + (5/128 + 1/64*n)*ε4 + 3/128*ε5,
(1/16 - 3/32*n + 1/32*n*n)*ε2 + (3/64 - 1/32*n - 3/64*n*n)*ε3 + (3/128 + 1/128*n)*ε4 + 5/256*ε5,
(5/192 - 3/64*n + 5/192*n*n)*ε3 + (3/128 - 5/192*n)*ε4 + 7/512*ε5,
(7/512 - 7/256*n)*ε4 + 7/512*ε5,
21/2560*ε5
];
const ΣC3σ1 = clenshawSinSeries(sinσ1, cosσ1, C3); // Σ C₃ₗsin2lσ
const I3σ1 = A3 * (σ1 + ΣC3σ1); // I₃(σ) = A₃·(σ + Σ C₃ₗsin2lσ)
const ΣC3σ2 = clenshawSinSeries(sinσ2, cosσ2, C3); // Σ C₃ₗsin2lσ
const I3σ2 = A3 * (σ2 + ΣC3σ2); // I₃(σ) = A₃·(σ + Σ C₃ₗsin2lσ)
// ω is the longitude on the auxiliary sphere
const ω1 = Math.atan2(sinα0*sinσ1, cosσ1); // ω₁ = Arg(cosσ + i·sinα₀sinσ)
const ω2 = Math.atan2(sinα0*sinσ2, cosσ2); // ω₂ = Arg(cosσ + i·sinα₀sinσ)
const λ1 = ω1 - f*sinα0*I3σ1; // λ = ω − f·sinα₀·0I₃(σ)
const λ2 = ω2 - f*sinα0*I3σ2; // λ = ω − f·sinα₀·0I₃(σ)
const λ12 = λ2 - λ1;
const φ2 = Math.atan(tanβ2 / (1 - f)); // tanβ = (1-f)·tanφ
const α2 = Math.atan2(sinα0, cosα0*cosσ2); // α₂ = final bearing
return {
point: new LatLon(φ2.toDegrees(), this.lon+λ12.toDegrees(), this.datum),
finalBearing: wrap360(α2.toDegrees()),
};
// [1]: qv Geodesics on an ellipsoid of revolution (Karney 2011), p8
};
/**
*
*/
LatLon.prototype.inverse = function(point) {
const p1 = this, p2 = point;
const φ1 = p1.lat.toRadians(), λ1 = p1.lon.toRadians();
const φ2 = p2.lat.toRadians(), λ2 = p2.lon.toRadians();
const a = this.datum.ellipsoid.a, b = this.datum.ellipsoid.b, f = this.datum.ellipsoid.f;
const e2 = f * (2-f); // eccentricity squared
const π = Math.PI;
// canonicalise configuration so that φ₁ ≤ 0, φ₁ ≤ φ₂ ≤ −φ₁, 0 ≤ λ₁₂ ≤ π
// TODO
const tanβ1 = (1-f) * Math.tan(φ1); // β = reduced latitude
const cosβ1 = 1 / Math.sqrt((1 + tanβ1*tanβ1)); // trig.identity
const sinβ1 = tanβ1 * cosβ1; // trig.identity
// from code:
//var sinβ1 = Math.sin(φ1) * (1-f);
//var cosβ1 = Math.cos(φ1);
//sinβ1 /= Math.hypot(sinβ1, cosβ1);
//cosβ1 /= Math.hypot(sinβ1, cosβ1);
const tanβ2 = (1-f) * Math.tan(φ2); // β = reduced latitude
const cosβ2 = 1 / Math.sqrt((1 + tanβ2*tanβ2)); // trig.identity
const sinβ2 = tanβ2 * cosβ2; // trig.identity
// from code:
//var sinβ2 = Math.sin(φ2) * (1-f);
//var cosβ2 = Math.cos(φ2);
//sinβ2 /= Math.hypot(sinβ1, cosβ1);
//cosβ2 /= Math.hypot(sinβ1, cosβ1);
const λ12 = λ2 - λ1;
const sinλ12 = Math.sin(λ12), cosλ12 = Math.cos(λ12);
// solve asteroid problem
const Δ = f * a * π * cosβ1*cosβ1;
const x = (λ12-π) * a*cosβ1 / Δ;
const β1 = Math.asin(sinβ1), β2 = Math.asin(sinβ2); // TODO: fudge
const y = (β1+β2) * a / Δ;
const μ = astroid(x, y);
//const α1 = Math.atan2(x/(1+μ), y/μ); // initial guess
const meridional = sinλ12 == 0 || p1.lat == -90; // code
const equatorial = sinβ1 == 0 && sinβ2 == 0; // code
// const meridional = λ12==0 || λ12==π; // paper
// const equatorial = φ1==0 && φ2==0 && λ12<=(1-f)*π; // paper
if (meridional) {
// TODO: α1 = λ12
}
if (equatorial) {
// TODO: α1 = π/2;
}
if (!meridional && !equatorial) { // general case
// TODO
}
//const sinα1 = Math.sin(α1);
//const cosα1 = Math.cos(α1);
//const sinα0 = sinα1*cosβ1;
const ώ = Math.sqrt(1 - e2 * Math.pow((cosβ1+cosβ2)/2, 2)); // √( 1 − e² ((cosβ +cosβ₂)/2)² ) TODO use ** (v8 5.1)
const ω12 = λ12 / ώ; // TODO table 3 assumption
const sinω12 = Math.sin(ω12), cosω12 = Math.cos(ω12);
const z1 = { x: cosβ1*sinβ2 - sinβ1*cosβ2*cosω12, y: cosβ2*sinω12 };
const z2 = { x: -sinβ1*cosβ2 + cosβ1*sinβ2*cosω12, y: cosβ2*sinω12 };
const α1 = Math.atan2(z1.y, z1.x); // TODO: refine by Newton's method...
const α2 = Math.atan2(z2.y, z2.x);
const σ12 = Math.atan(Math.hypot(z1.x, z1.y), sinβ1*sinβ2 + cosβ1*cosβ2*cosω12);
// TODO: provisionally, ...
const s12 = a*ώ*σ12;
const sinα1 = Math.sin(α1), cosα1 = Math.cos(α1); // TODO: !!
// α0 is the azimuth of the geodesic at the equator
const sinα0 = sinα1 * cosβ1; // Clairaut’s relation
const cosα0 = Math.hypot(cosα1, sinα1*sinβ1); // [1]
const sinα2 = sinα0 / cosβ2; // Clairaut’s relation
// TODO cosα2 = ???
console.log('β1', Math.atan(tanβ1).toDegrees());
console.log('β2', Math.atan(tanβ2).toDegrees());
console.log('ώ', ώ);
console.log('ω12', ω12.toDegrees());
console.log('σ12', σ12.toDegrees());
console.log('α1', α1.toDegrees());
console.log('α2', α2.toDegrees());
console.log('s12', s12);
console.log('x', x);
console.log('y', y);
console.log('μ', μ);
let s = s12;
s = Number(s.toFixed(6)); // round to μm precision
return { distance: s, initialBearing: α1.toDegrees(), finalBearing: α2.toDegrees() };
};
// www.researchgate.net/publication/242330657_SOME_APPLICATIONS_OF_CLENSHAW'S_RECURRENCE_FORMULA_IN_MAP_PROJECTIONS
// geographiclib.sourceforge.net/html/Geodesic_8cpp_source.html
// github.com/devbharat/gtsam/blob/master/gtsam/3rdparty/GeographicLib/matlab/private/SinCosSeries.m
// return σ12, sinα1, cosα1, sinα2, cosα2, dnm
function inverseStartPoint(sinβ1, cosβ1, dn1,
sinβ2, cosβ2, dn2, λ12,
C1a, C2a) {
// Return a starting point for Newton's method in sinα1 and cosα1
// (function value is -1). If Newton's method doesn't need to be
// used, return also sinα2 and cosα2 and function value is σ12.
// sinα2, cosα2 only updated if return val >= 0.
var vals = {};
// β12 = β2 - β1 in [0, pi); β12a = β2 + β1 in (-pi, 0]
const sinβ12 = sinβ2 * cosβ1 - cosβ2 * sinβ1;
const cosβ12 = cosβ2 * cosβ1 + sinβ2 * sinβ1;
const sinβ12ʹ = sinβ2 * cosβ1 + cosβ2 * sinβ1; // TODO ??
var shortline, ω12, sinβm2, sinω12, cosω12, t, sinσ12, cosσ12,
x, y, lamscale, betscale, k2, eps, cosβ12a, β12a, m12b, m0, nvals,
k, ω12a;
vals.σ12 = -1; // Return value
// Volatile declaration needed to fix inverse cases
// 88.202499451857 0 -88.202499451857 179.981022032992859592
// 89.262080389218 0 -89.262080389218 179.992207982775375662
// 89.333123580033 0 -89.333123580032997687 179.99295812360148422
// which otherwise fail with g++ 4.4.4 x86 -O3
shortline = cosβ12 >= 0 && sinβ12 < 0.5 && cosβ2 * λ12 < 0.5;
ω12 = λ12;
if (shortline) {
// sin((β1+β2)/2)² = (sinβ1 + sinβ2)² / ((sinβ1 + sinβ2)² + (cosβ1 + cosβ2)²)
sinβm2 = Math.pow(sinβ1 + sinβ2, 2) / ( Math.pow(sinβ1 + sinβ2, 2) + Math.pow(cosβ1 + cosβ2, 2) );
vals.dnm = Math.sqrt(1 + this._ep2 * sinβm2);
ω12 = λ12 / (1-f) * vals.dnm;
}
sinω12 = Math.sin(ω12); cosω12 = Math.cos(ω12);
vals.sinα1 = cosβ2 * sinω12;
vals.cosα1 = cosω12 >= 0 ?
sinβ12 + cosβ2 * sinβ1 * sinω12*sinω12 / (1 + cosω12) :
sinβ12ʹ - cosβ2 * sinβ1 * sinω12*sinω12 / (1 - cosω12);
sinσ12 = m.hypot(vals.sinα1, vals.cosα1);
cosσ12 = sinβ1 * sinβ2 + cosβ1 * cosβ2 * cosω12;
if (shortline && sinσ12 < this._etol2) {
// really short lines
vals.sinα2 = cosβ1 * sinω12;
vals.cosα2 = sinβ12 - cosβ1 * sinβ2 *
(cosω12 >= 0 ? m.sq(sinω12) / (1 + cosω12) : 1 - cosω12);
// norm(vals.sinα2, vals.cosα2);
t = m.hypot(vals.sinα2, vals.cosα2); vals.sinα2 /= t; vals.cosα2 /= t;
// Set return value
vals.σ12 = Math.atan2(sinσ12, cosσ12);
} else if (Math.abs(this._n) > 0.1 || // Skip astroid calc if too eccentric
cosσ12 >= 0 ||
sinσ12 >= 6 * Math.abs(this._n) * Math.PI * m.sq(cosβ1)) {
// Nothing to do, zeroth order spherical approximation is OK
} else {
// Scale λ12 and β2 to x, y coordinate system where antipodal
// point is at origin and singular point is at y = 0, x = -1.
if (this.f >= 0) { // In fact f == 0 does not get here
// x = dlong, y = dlat
k2 = m.sq(sinβ1) * this._ep2;
eps = k2 / (2 * (1 + Math.sqrt(1 + k2)) + k2);
lamscale = this.f * cosβ1 * this.A3f(eps) * Math.PI;
betscale = lamscale * cosβ1;
x = (λ12 - Math.PI) / lamscale;
y = sinβ12ʹ / betscale;
} else { // f < 0
// x = dlat, y = dlong
cosβ12a = cosβ2 * cosβ1 - sinβ2 * sinβ1;
β12a = Math.atan2(sinβ12ʹ, cosβ12a);
// In the case of lon12 = 180, this repeats a calculation made
// in Inverse.
nvals = this.Lengths(this._n, Math.PI + β12a,
sinβ1, -cosβ1, dn1, sinβ2, cosβ2, dn2,
cosβ1, cosβ2, g.REDUCEDLENGTH, C1a, C2a);
m12b = nvals.m12b; m0 = nvals.m0;
x = -1 + m12b / (cosβ1 * cosβ2 * m0 * Math.PI);
betscale = x < -0.01 ? sinβ12ʹ / x :
-this.f * m.sq(cosβ1) * Math.PI;
lamscale = betscale / cosβ1;
y = (λ12 - Math.PI) / lamscale;
}
if (y > -tol1_ && x > -1 - xthresh_) {
// strip near cut
if (this.f >= 0) {
vals.sinα1 = Math.min(1, -x);
vals.cosα1 = - Math.sqrt(1 - m.sq(vals.sinα1));
} else {
vals.cosα1 = Math.max(x > -tol1_ ? 0 : -1, x);
vals.sinα1 = Math.sqrt(1 - m.sq(vals.cosα1));
}
} else {
// Estimate α1, by solving the astroid problem.
//
// Could estimate alpha1 = theta + pi/2, directly, i.e.,
// cosα1 = y/k; sinα1 = -x/(1+k); for f >= 0
// cosα1 = x/(1+k); sinα1 = -y/k; for f < 0 (need to check)
//
// However, it's better to estimate ω12 from astroid and use
// spherical formula to compute α1. This reduces the mean number of
// Newton iterations for astroid cases from 2.24 (min 0, max 6) to 2.12
// (min 0 max 5). The changes in the number of iterations are as
// follows:
//
// change percent
// 1 5
// 0 78
// -1 16
// -2 0.6
// -3 0.04
// -4 0.002
//
// The histogram of iterations is (m = number of iterations estimating
// α1 directly, n = number of iterations estimating via ω12, total
// number of trials = 148605):
//
// iter m n
// 0 148 186
// 1 13046 13845
// 2 93315 102225
// 3 36189 32341
// 4 5396 7
// 5 455 1
// 6 56 0
//
// Because ω12 is near pi, estimate work with ω12a = pi - ω12
k = Astroid(x, y);
ω12a = lamscale * ( this.f >= 0 ? -x * k/(1 + k) : -y * (1 + k)/k );
sinω12 = Math.sin(ω12a); cosω12 = -Math.cos(ω12a);
// Update spherical estimate of α1 using ω12 instead of
// λ12
vals.sinα1 = cosβ2 * sinω12;
vals.cosα1 = sinβ12ʹ -
cosβ2 * sinβ1 * m.sq(sinω12) / (1 - cosω12);
}
}
// Sanity check on starting guess. Backwards check allows NaN through.
if (!(vals.sinα1 <= 0)) {
// norm(vals.sinα1, vals.cosα1);
t = Math.hypot(vals.sinα1, vals.cosα1); vals.sinα1 /= t; vals.cosα1 /= t;
} else {
vals.sinα1 = 1; vals.cosα1 = 0;
}
return vals;
}
/**
* Clenshaw summation to sum trigonometric sin series ∑ Cₙsin2nθ
*
* @param {number} sinθ - Value of sin θ.
* @param {number} cosθ - Value of cos θ.
* @param {number[]} C - Series expansion coefficients (1-based array).
* @returns {number} Series sum.
*/
function clenshawSinSeries(sinθ, cosθ, C) {
const a = 2 * (cosθ-sinθ) * (cosθ+sinθ); // 2·cos(2θ)
let N = C.length-1; // O(N)
let y1 = N%2==1 ? C[N--] : 0;
let y2 = 0;
for (let k=N; k>0; k-=2) {
y2 = a*y1 - y2 + C[k];
y1 = a*y2 - y1 + C[k-1];
}
return 2 * sinθ * cosθ * y1;
}
/**
* Clenshaw summation to sum trigonometric cos series ∑ Cₙcos(2n+1)θ
*
* @param {number} sinθ - Value of sin θ.
* @param {number} cosθ - Value of cos θ.
* @param {number[]} C - Series expansion coefficients (0-based array).
* @returns {number} Series sum.
*/
function clenshawCosSeries(sinθ, cosθ, C) {
const a = 2 * (cosθ-sinθ) * (cosθ+sinθ); // 2·cos(2θ)
let N = C.length-1; // O(N)
let y1 = N%2==0 ? C[N--] : 0;
let y2 = 0;
for (let k=N; k>0; k-=2) {
y2 = a*y1 - y2 + C[k];
y1 = a*y2 - y1 + C[k-1];
}
return cosθ * (y1-y2);
}
/*
* Solve μ⁴ + 2μ³ + (1-x²-y²)μ² − 2y²μ − y² = 0 for +ve root μ.
*
* From GeographicLib Geocentric::Reverse; qv dlmf.nist.gov/1.11#ii.
*/
function astroid(x, y) {
const p = x*x;
const q = y*y;
const r = (p + q - 1) / 6;
let μ = null;
if ( !(q == 0 && r <= 0) ) {
// avoid possible division by zero when r = 0 by multiplying equations for s and t by r³ and r, resp.
const S = p * q / 4; // S = r³·s
const r2 = r * r; // r²
const r3 = r * r2; // r³
// discriminant of the quadratic equation for T3; this is zero on the evolute curve p^⅓+q^⅓ = 1
const discriminant = S * (S + 2*r3);
let u = r;
if (discriminant >= 0) {
// pick the sign on the sqrt to maximize |T³|, to minimize loss of precision due to
// cancellation; result is unchanged because of the way the T is used in definition of u.
const T3 = S+r3 < 0 ? S+r3 - Math.sqrt(discriminant) : S+r3 + Math.sqrt(discriminant); // T³ = (r·t)³
const T = Math.cbrt(T3); // T = r·t; T can be zero; but then r²/T -> 0.
u += T + (T != 0 ? r2/T : 0);
} else {
// T is complex, but the way u is defined the result is real
const θ = Math.atan2(Math.sqrt(-discriminant), -(S+r3));
// there are three possible cube roots; choose the root which avoids cancellation
// note that discriminant < 0 implies that r < 0.
u += 2 * r * Math.cos(θ/3);
}
const v = Math.sqrt(u*u + q); // (guaranteed positive)
// avoid loss of accuracy when u < 0.
const uv = u < 0 ? q / (v - u) : u + v; // u+v, guaranteed positive
const w = (uv - q) / (2 * v); // positive?
// rearrange expression for μ to avoid loss of accuracy due to subtraction;
// division by 0 not possible because uv > 0, w >= 0
μ = uv / (Math.sqrt(uv + w*w) + w); // guaranteed positive
} else { // q == 0 && r <= 0
// y = 0 with |x| <= 1; handle this case directly; for y small, +ve root is μ = |y|/√(1-x²)
μ = 0;
}
return μ;
}
function wrap360(degrees) { return (degrees%360+360) % 360; }
function wrap180(degrees) { return (degrees%360+180) % 360 - 180; }
/* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
/** Extend Number object with method to convert numeric degrees to radians */
if (Number.prototype.toRadians === undefined) {
Number.prototype.toRadians = function() { return this * Math.PI / 180; };
}
/** Extend Number object with method to convert radians to numeric (signed) degrees */
if (Number.prototype.toDegrees === undefined) {
Number.prototype.toDegrees = function() { return this * 180 / Math.PI; };
}
/* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
if (typeof module != 'undefined' && module.exports) module.exports = LatLon; // CommonJS (Node)
if (typeof define == 'function' && define.amd) define([], function() { return LatLon; }); // AMD