Source: x-latlon-geodesic-karney-b4-integral-functions.js

'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.initialBearingTo = function(point) {
    if (!(point instanceof LatLon)) throw new TypeError('point is not LatLon object');
    return this.inverse(point).initialBearing;
};

LatLon.prototype.finalBearingTo = function(point) {
    if (!(point instanceof LatLon)) throw new TypeError('point is not LatLon object');
    return this.inverse(point).finalBearing;
};

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;
    let φ1 = p1.lat.toRadians();
    let φ2 = p2.lat.toRadians();
    let λ12 = p2.lon.toRadians() - p1.lon.toRadians();

    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 π = Math.PI;

    // canonicalise configuration so that φ₁ ≤ 0, φ₁ ≤ φ₂ ≤ −φ₁, 0 ≤ λ₁₂ ≤ π
    let txLon = λ12>=0 ? 1 : -1;
    λ12 *= txLon;
    let txPts = Math.abs(φ1) < Math.abs(φ2) ? -1 : 1;
    if (txPts < 0) {
        txLon = -txLon;
        [φ1,φ2] = [φ2,φ1]; // swap φ1 & φ2
    }
    let txLat = φ1<0 ? 1 : -1;
    φ1 *= txLat, φ2 *= txLat;

    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 sinλ12 = Math.sin(λ12), cosλ12 = Math.cos(λ12);


    // solve asteroid problem; define plane coordinate system (x,y) centered on antipodal point
    // where Δ = f·a·π·cos²β1 is the unit of length
    const Δ = f * a * π * cosβ1*cosβ1;
    const x = (λ12-π) * a*cosβ1 / Δ;
    const β1 = Math.atan2(sinβ1, cosβ1), β2 = Math.atan2(sinβ2, cosβ2); // TODO: fudge for y
    const y = (β1+β2) * a / Δ;
    const μ = astroid(x, y);

    const α1 = Math.atan2(-x/(1+μ), y/μ); // initial guess
    let sinα1 = Math.sin(α1), cosα1 = Math.cos(α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);
    const σ1 = Math.atan2(sinσ1, cosσ1);               // TODO: where to declare?

    // α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 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 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 sinα2 = sinα0 / cosβ2;                                // Clairaut’s relation
    const cosα2 = Math.sqrt(cosα1*cosα1 * cosβ1*cosβ1 + (cosβ2*cosβ2 - cosβ1*cosβ1)) / cosβ2;
    const α2 = Math.atan2(sinα2, cosα2)

    const σ2 = Math.atan2(sinβ2, cosα2*cosβ2);
    const sinσ2 = Math.sin(σ2), cosσ2 = Math.cos(σ2);

    //const ω2 = Math.atan2(sinα0*sinσ2, cosσ2);

    // ε 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 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 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 I1σ1 = A1 * (σ1 + ΣC1σ1);                             // I₁(σ): distance integral for s/b

    const A2 = (1 + 1/4*ε2 + 9/64*ε4 + 25/256*ε6) * (1-ε);
    const C2 = [ null,
         1/2*ε  + 1/16*ε3 +    1/32*ε5,
        3/16*ε2 + 1/32*ε4 + 35/2048*ε6,
                  5/48*ε3 +   5/256*ε5,
                35/512*ε4 +   7/512*ε6,
                            63/1280*ε5,
                            77/2048*ε6,
    ];
    const ΣC2σ1 = clenshawSinSeries(sinσ1, cosσ1, C2);          // Σ C₁ₗsin2lσ
    const I2σ1 = A1 * (σ2 + ΣC1σ1);                             // I₁(σ): distance integral for s/b

    const Jσ1 = I1σ1 - I2σ1;














    //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;


    // TODO cosα2 = ???



    console.log('txPts', txPts, 'txLat', txLat, 'txLon', txLon)
    console.log('x', x);
    console.log('y', y);
    console.log('μ', μ);
    console.log('α1', α1.toDegrees());

    console.log('β1', Math.atan(tanβ1).toDegrees());
    console.log('α0', Math.atan2(sinα0, cosα0).toDegrees());
    console.log('σ1', Math.atan2(sinσ1, cosσ1).toDegrees());
    console.log('ω1', ω1.toDegrees());

    console.log('β2', Math.atan(tanβ2).toDegrees());
    console.log('α2', α2.toDegrees());
    console.log('σ2', σ2.toDegrees());
    console.log('ω2', ω2.toDegrees());

    console.log('k2', k2);
    console.log('ε', ε);

    console.log('λ1', λ1.toDegrees());
    console.log('λ2', λ2.toDegrees());
    console.log('λ12', (λ2-λ1).toDegrees());

    console.log('δλ12', (λ2-λ1-λ12).toDegrees());
    console.log('Jσ1', Jσ1);





    //console.log('ω12', ω12.toDegrees());
    //console.log('σ12', σ12.toDegrees());
    //console.log('s12', s12);



    if (txPts < 0) {
        [sinα1,sinα2] = [sinα2,sinα1]; // swap sinα1 & sinα2
        [cosα1,cosα2] = [cosα2,cosα1]; // swap cosα1 & cosα2
    }
    //sinα1 *= txPts * txLon, cosα1 *= txPts * txLon;
    //sinα2 *= txPts * txLon, cosα2 *= txPts * txLon;


    //let s = s12;
    let s = 0; // TODO tmp
    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