/*
 * Licensed to the Apache Software Foundation (ASF) under one or more
 * contributor license agreements.  See the NOTICE file distributed with
 * this work for additional information regarding copyright ownership.
 * The ASF licenses this file to You under the Apache License, Version 2.0
 * (the "License"); you may not use this file except in compliance with
 * the License.  You may obtain a copy of the License at
 *
 *     http://www.apache.org/licenses/LICENSE-2.0
 *
 * Unless required by applicable law or agreed to in writing, software
 * distributed under the License is distributed on an "AS IS" BASIS,
 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
 * See the License for the specific language governing permissions and
 * limitations under the License.
 */
package org.apache.sis.referencing;

import org.opengis.geometry.coordinate.Position;
import org.opengis.referencing.datum.Ellipsoid;
import org.opengis.referencing.crs.CoordinateReferenceSystem;
import org.apache.sis.internal.referencing.Resources;
import org.apache.sis.internal.referencing.Formulas;
import org.apache.sis.math.MathFunctions;

import static java.lang.Math.*;


/**
 * Performs geodetic calculations on an ellipsoid. This class overrides the spherical
 * formulas implemented in the parent class, replacing them by ellipsoidal formulas.
 * The methods for direct and inverse geodesic problem use the formulas described in
 * the following publication:
 *
 * <blockquote>
 * Charles F. F. Karney, 2013.
 * <a href="https://doi.org/10.1007/s00190-012-0578-z">Algorithms for geodesics</a>, SRI International.
 * </blockquote>
 *
 * The following symbols are used with the same meaning than in Karney's article,
 * except λ₁₂ which is represented by ∆λ:
 *
 * <ul>
 *   <li><var>a</var>:  equatorial radius of the ellipsoid of revolution.</li>
 *   <li><var>b</var>:  polar semi-axis of the ellipsoid of revolution.</li>
 *   <li><var>ℯ</var>:  first  eccentricity: ℯ  = √[(a²-b²)/a²].</li>
 *   <li><var>ℯ′</var>: second eccentricity: ℯ′ = √[(a²-b²)/b²].</li>
 *   <li><var>n</var>:  third flattening:    n  = (a-b)/(a+b).</li>
 *   <li><var>E</var>:  the point at which the geodesic crosses the equator in the northward direction.</li>
 *   <li><var>P</var>:  the start point (P₁) or end point (P₂).</li>
 *   <li><var>∆λ</var>: longitude difference between start point and end point (λ₁₂ in Karney).</li>
 *   <li><var>β</var>:  reduced latitude, related to φ geodetic latitude on the ellipsoid.</li>
 *   <li><var>ω</var>:  spherical longitude, related to λ geodetic longitude on the ellipsoid.</li>
 *   <li><var>σ</var>:  spherical arc length, related to distance <var>s</var> on the ellipsoid.</li>
 * </ul>
 *
 * Suffix 1 in variable names denotes values computed using starting point (P₁) and starting azimuth (α₁)
 * while suffix 2 denotes values computed using ending point (P₂) and ending azimuth (α₂).
 * All angular values stored in this class are in radians.
 *
 * <p><b>Limitations:</b>
 * Current implementation is still unable to compute the geodesics in some cases.
 * In particular, calculation may fail for antipodal points.
 * See <a href="https://issues.apache.org/jira/browse/SIS-467">SIS-467</a>.</p>
 *
 * <p>If the following cases where more than one geodesics exist, current implementation returns an arbitrary one:</p>
 * <ul>
 *   <li>Coincident points (distance is zero but azimuths can be anything).</li>
 *   <li>Starting point and ending points are at opposite poles (there are infinitely many geodesics).</li>
 *   <li>φ₁ = -φ₂ and ∆λ is close to 180° (two geodesics may exist).</li>
 *   <li>∆λ = ±180° (two geodesics may exist).</li>
 * </ul>
 *
 * @author  Matthieu Bastianelli (Geomatys)
 * @author  Martin Desruisseaux (Geomatys)
 * @version 1.0
 * @since   1.0
 * @module
 */
class GeodesicsOnEllipsoid extends GeodeticCalculator {
    /**
     * Whether to include code used for JUnit tests only. This field should be
     * set to {@code true} during development and to {@code false} in releases.
     *
     * @see #snapshot()
     */
    static final boolean STORE_LOCAL_VARIABLES = false;

    /**
     * Accuracy threshold iterative computations, in radians.
     * We take a finer accuracy than default SIS configuration in order to met the accuracy of numbers
     * published in Karney (2013). If this value is modified, the effect can be verified by executing
     * the {@code GeodesicsOnEllipsoidTest} methods that compare computed values against Karney's tables.
     * Remember to update {@link GeodeticCalculator} class javadoc if this value is changed.
     *
     * <p><b>Note:</b> when the iteration loop detects that it reached this requested accuracy, the loop
     * completes the iteration step which was in progress. Consequently the final accuracy is one iteration
     * better than the accuracy computed from this value.</p>
     */
    static final double ITERATION_TOLERANCE = Formulas.ANGULAR_TOLERANCE * (PI/180) / 20;

    /**
     * Difference between ending point and antipode of starting point for considering them as nearly antipodal.
     * This is used only for finding an approximate position before to start iteration using Newton's method,
     * so it is okay if the antipodal approximation has some inaccuracy.
     */
    private static final double NEARLY_ANTIPODAL_Δλ = 0.25 * (PI/180);

    /**
     * The square of eccentricity: ℯ² = (a²-b²)/a² where
     * <var>ℯ</var> is the eccentricity,
     * <var>a</var> is the <cite>semi-major</cite> axis length and
     * <var>b</var> is the <cite>semi-minor</cite> axis length.
     */
    final double eccentricitySquared;

    /**
     * The square of the second eccentricity: ℯ′² = (a²-b²)/b².
     */
    final double secondEccentricitySquared;

    /**
     * Third flattening of the ellipsoid: n = (a-b)/(a+b).
     * Not to be confused with the third eccentricity squared ℯ″² = (a²-b²)/(a²+b²) which is not used in this class.
     */
    final double thirdFlattening;

    /**
     * Ration between the semi-minor and semi-major axis b/a.
     * This is related to flattening <var>f</var> = 1 - b/a.
     */
    final double axisRatio;

    /**
     * Length of semi-minor axis.
     *
     * @see org.opengis.referencing.datum.Ellipsoid#getSemiMinorAxis()
     */
    private double semiMinorAxis() {
        return semiMajorAxis * axisRatio;
    }

    /**
     * The α value computed from the starting point and starting azimuth.
     * We use the sine and cosine instead than the angles because those components are more frequently used than angles.
     * Those values can be kept constant when computing many end points and end azimuths at different geodesic distances.
     * The {@link #COEFFICIENTS_FOR_START_POINT} flag specifies whether those fields need to be recomputed.
     */
    private double sinα0, cosα0;

    /**
     * Longitude angle from the equatorial point <var>E</var> to starting point P₁ on the ellipsoid.
     * This longitude is computed from the ω₁ longitude on auxiliary sphere, which is itself computed
     * from α₀, α₁, β₁ and ω₁ values computed from the starting point and starting azimuth.
     * The {@link #COEFFICIENTS_FOR_START_POINT} flag specifies whether this field needs to be recomputed.
     *
     * @see #sphericalToGeodeticLongitude(double, double)
     */
    private double λ1E;

    /**
     * Ellipsoidal arc length <var>s₁</var>/<var>b</var> computed from the spherical arc length <var>σ₁</var>.
     * The <var>σ₁</var> value is an arc length on the auxiliary sphere between equatorial point <var>E</var>
     * (the point on equator with forward direction toward azimuth α₀) and starting point P₁. This is computed
     * by I₁(σ₁) in Karney equation 15 and is saved for reuse when the starting point and azimuth do not change.
     * The {@link #COEFFICIENTS_FOR_START_POINT} flag specifies whether this field needs to be recomputed.
     *
     * @see #sphericalToEllipsoidalAngle(double, boolean)
     */
    private double I1_σ1;

    /**
     * The term to be raised to powers (ε⁰, ε¹, ε², ε³, …) in series expansions. Defined in Karney equation 16 as
     * ε = (√[k²+1] - 1) / (√[k²+1] + 1) where k² = {@linkplain #secondEccentricitySquared ℯ′²}⋅cos²α₀ (Karney 9).
     * The {@link #COEFFICIENTS_FOR_START_POINT} flag specifies whether this field needs to be recomputed.
     */
    private double ε;

    /**
     * Coefficients in series expansion of the form A × (σ + ∑Cₓ⋅sin(…⋅σ)).
     * There is 3 series expansions used in this class:
     *
     * <ul>
     *   <li><var>A₁</var> and <var>C₁ₓ</var> (Karney 15) are used for arc length conversions from auxiliary sphere to ellipsoid.</li>
     *   <li><var>A₂</var> and <var>C₂ₓ</var> (Karney 41) are used in calculation of reduced length.</li>
     *   <li><var>A₃</var> and <var>C₃ₓ</var> (Karney 23) are used for calculation of geodetic longitude from spherical angles.</li>
     * </ul>
     *
     * The <var>Cₓₓ</var> coefficients are hard-coded, except the <var>C₃ₓ</var> coefficients
     * (used together with <var>A₃</var>) which depend on {@link #sinα0} and {@link #cosα0}.
     * Note that the <var>C</var> coefficients differ from the ones published by Karney because
     * they have been combined using Clenshaw summation.
     *
     * <p>All those coefficients must be recomputed when the starting point or starting azimuth change.
     * The {@link #COEFFICIENTS_FOR_START_POINT} flag specifies whether those fields need to be recomputed.</p>
     *
     * @see #computeSeriesExpansionCoefficients()
     */
    private double A1, A2, A3, C31, C32, C33, C34, C35;

    /**
     * Coefficients for Rhumb line calculation from Bennett (1996) equation 2, modified with Clenshaw summation.
     */
    private final double R0, R2, R4, R6;

    /**
     * Constructs a new geodetic calculator expecting coordinates in the supplied CRS.
     *
     * @param  crs         the referencing system for the {@link Position} arguments and return values.
     * @param  ellipsoid   ellipsoid associated to the geodetic component of given CRS.
     */
    GeodesicsOnEllipsoid(final CoordinateReferenceSystem crs, final Ellipsoid ellipsoid) {
        super(crs, ellipsoid);
        final double a = semiMajorAxis;
        final double b = ellipsoid.getSemiMinorAxis();
        final double Δ2 = a*a - b*b;
        eccentricitySquared = Δ2 / (a*a);
        secondEccentricitySquared = Δ2 / (b*b);
        thirdFlattening = (a - b) / (a + b);
        axisRatio = b / a;
        /*
         * For rhumb-line calculation from Bennett (1996) equation 2, modified with Clenshaw summation.
         */
        double fe;
        R0 = 1  -  (eccentricitySquared)*(  1./4   + eccentricitySquared*( 3./64 + eccentricitySquared*( 5./256)));
        R2 = (fe  = eccentricitySquared)*( -3./8   - eccentricitySquared*( 3./32 + eccentricitySquared*(25./768)));
        R4 = (fe *= eccentricitySquared)*( 15./128 + eccentricitySquared*(45./512));
        R6 = (fe *  eccentricitySquared)*(-35./768);
    }

    /**
     * Computes series expansions coefficients.
     *
     * <p><b>Preconditions:</b> The {@link #sinα0} and {@link #cosα0} fields shall be set before to invoke this method.
     * It is caller's responsibility to ensure that sin(α₀)² + cos(α₀)² ≈ 1 (this is verified in assertion).</p>
     *
     * <p><b>Post-conditions:</b> this method sets the {@link #ε}, {@link #A1}, {@link #A2}, {@link #A3},
     * {@link #C31}, {@link #C32}, {@link #C33}, {@link #C34}, {@link #C35} fields.</p>
     *
     * It is caller's responsibility to invoke {@code setValid(COEFFICIENTS_FOR_START_POINT)} after it has updated
     * {@link #λ1E} and {@link #I1_σ1}.
     *
     * @return k² = ℯ′²⋅cos²α₀ term (Karney 9).
     */
    private double computeSeriesExpansionCoefficients() {
        assert abs(sinα0*sinα0 + cosα0*cosα0 - 1) <= 1E-10;                     // Arbitrary threshold.
        final double k2 = secondEccentricitySquared * (cosα0 * cosα0);          // (Karney 9)
        final double h = sqrt(1 + k2);
        ε = (h - 1) / (h + 1);                                                  // (Karney 16)
        final double ε2 = ε*ε;
        A1 = (((( 1./256)*ε2 + 1./64)*ε2 + 1./4)*ε2 + 1) / (1 - ε);             // (Karney 17)
        A2 = ((((25./256)*ε2 + 9./64)*ε2 + 1./4)*ε2 + 1) * (1 - ε);             // (Karney 42)
        /*
         * Compute A₃ from Karney equation 24 with signs redistributed in a different but equivalent way.
         * We use + or - operations in a way where multiplication factors of n are positive. By increasing
         * the number of times where the same value is used (for example a single 1/8 instead of two values
         * +1/8 and -1/8), our hope is that the compiler can reuse more often values that are already loaded
         * in registers. Maybe some compilers could even detect when the same multiplication is done two times.
         *
         * Note: for each line below, we could compute at construction time the parts at the right of ε.
         * However since they are only a few additions and multiplications, it may not be worth to vastly
         * increase the number of fields in `GeodesicOnEllipsoid` for that purpose, especially if compiler
         * can optimize itself those duplicated multiplications.
         */
        final double n = thirdFlattening;
        A3 = 1 - ε*(1./2    -   n*(1./2)
               + ε*(1./4    +   n*(1./8    -   n*(3./8 ))
               + ε*(1./16   +   n*(3./16   +   n*(1./16))
               + ε*(3./64   +   n*(1./32)
               + ε*(3./128)))));
        /*
         * Karney equation 25 with fε = εⁿ where n = 1, 2, 3, …. The coefficients below differ from
         * the ones published by Karney because they have been combined using Clenshaw summation
         * (see sphericalToEllipsoidalAngle(…) below for a simpler example for Clenshaw summation).
         */
        double fε;
        C31 = (fε  = ε) * ( 1./4     -   n*(1./4)
                   + ε  * ( 1./8     -   n*(            n*(1./8 ))
                   + ε  * ( 1./48    +   n*(3./32   -   n*(1./24))
                   + ε  * ( 1./64    +   n*(1./24)
                   + ε  * (23./1280)))));
        C32 = (fε *= ε) * ( 1./8     -   n*(3./16   -   n*(1./16))
                   + ε  * ( 3./32    -   n*(1./16   +   n*(3./32))
                   + ε  * (-1./128   +   n*(1./8)
                   + ε  * (-1./64))));
        C33 = (fε *= ε) * ( 5./48    -   n*(3./16   -   n*(5./48))
                   + ε  * ( 3./32    -   n*(5./48)
                   + ε  * (-7./160)));
        C34 = (fε *= ε) * ( 7./64    -   n*(7./32)
                   + ε  * ( 7./64));
        C35 = (fε *  ε) * (21./160);
        return k2;
    }

    /**
     * Computes a term in calculation of ellipsoidal arc length <var>s</var> from the given spherical arc length <var>σ</var>.
     * The <var>σ</var> value is an arc length on the auxiliary sphere between equatorial point <var>E</var>
     * (the point on equator with forward direction toward azimuth α₀) and another point.
     *
     * <p>If the {@code minusI2} argument is {@code false}, then this method computes the ellipsoidal arc length <var>s</var>
     * from the given spherical arc length <var>σ</var>. This value is computed by Karney's equation 7:
     *
     * <blockquote>s/b = I₁(σ)</blockquote>
     *
     * If the {@code minusI2} argument is {@code true}, then this method is modified for computing
     * Karney's equation 40 instead (a term in calculation of reduced length <var>m</var>):
     *
     * <blockquote>J(σ) = I₁(σ) - I₂(σ)</blockquote>
     *
     * <p><b>Precondition:</b> {@link #computeSeriesExpansionCoefficients()} must have been invoked before this method.</p>
     *
     * @param  σ        arc length on the auxiliary sphere since equatorial point <var>E</var>.
     * @param  minusI2  whether to subtract I₂(σ) after calculation of I₁(σ).
     * @return I₁(σ) if {@code minusI2} is false, or J(σ) if {@code minusI2} is true.
     */
    private double sphericalToEllipsoidalAngle(final double σ, final boolean minusI2) {
        /*
         * Derived from Karney (2013) equation 18 with θ = 2σ:
         *
         *   ∑Cₓ⋅sin(x⋅θ)  ≈  (-1/2 ⋅ε   +  3/16 ⋅ε³  +  -1/32  ⋅ε⁵)⋅sin(1θ) +                  // C₁₁
         *                    (-1/16⋅ε²  +  1/32 ⋅ε⁴  +  -9/2048⋅ε⁶)⋅sin(2θ) +                  // C₁₂
         *                    (            -1/48 ⋅ε³  +   3/256 ⋅ε⁵)⋅sin(3θ) +                  // C₁₃
         *                    (            -5/512⋅ε⁴  +   3/512 ⋅ε⁶)⋅sin(4θ) +                  // C₁₄
         *                    (                          -7/1280⋅ε⁵)⋅sin(5θ) +                  // C₁₅
         *                    (                          -7/2048⋅ε⁶)⋅sin(6θ)                    // C₁₆
         *
         * For performance and simplicity, we rewrite the equations using trigonometric identities
         * as described in Snyder 3-34 and 3-35, expanded with two more terms. Given:
         *
         *   s  =  A⋅sin(θ) + B⋅sin(2θ) + C⋅sin(3θ) + D⋅sin(4θ) + E⋅sin(5θ) + F⋅sin(6θ)
         *
         * We rewrite as:
         *
         *   s  =  sinθ⋅(A′ + cosθ⋅(B′ + cosθ⋅(C′ + cosθ⋅(D′ + cosθ⋅(E′ + cosθ⋅F′)))))
         *   A′ =  A - C + E
         *   B′ =  2B - 4D + 6F
         *   C′ =  4C - 12E
         *   D′ =  8D - 32F
         *   E′ = 16E
         *   F′ = 32F
         *
         * See https://svn.apache.org/repos/asf/sis/analysis/Map%20projection%20formulas.ods for calculations.
         * Additions are done with smallest values first for better precision.
         *
         * See Karney 2011, Geodesics on an ellipsoid of revolution, https://arxiv.org/pdf/1102.1215.pdf
         * if more coefficients (up to ε⁸) are desired. Those coefficients may be needed in a future SIS
         * version if we want to support celestial bodies with higher flattening factor.
         * Issue tracker: https://issues.apache.org/jira/browse/SIS-465
         */
        final double ε2 = ε*ε;
        final double ε3 = ε*ε2;
        final double ε4 = ε2*ε2;
        final double ε5 = ε2*ε3;
        final double ε6 = ε3*ε3;
        final double θ  = σ * 2;
        final double sinθ = sin(θ);
        final double cosθ = cos(θ);
        double m = A1 * (σ + sinθ * (-31./640 * ε5   +    5./24  * ε3   +   -1./2 * ε
                           + cosθ * (-27./512 * ε6   +   13./128 * ε4   +   -1./8 * ε2
                           + cosθ * (  9./80  * ε5   +   -1./12  * ε3
                           + cosθ * (  5./32  * ε6   +   -5./64  * ε4
                           + cosθ * ( -7./80  * ε5
                           + cosθ * ( -7./64  * ε6)))))));
        if (minusI2) {
            m  -=  A2 * (σ + sinθ * ( 39./640 * ε5   +   -1./24  * ε3   +    1./2 * ε
                           + cosθ * (105./512 * ε6   +  -27./128 * ε4   +    3./8 * ε2
                           + cosθ * ( -41./80 * ε5   +    5./12  * ε3
                           + cosθ * ( -35./32 * ε6   +   35./64  * ε4
                           + cosθ * (  63./80 * ε5
                           + cosθ * (  77./64 * ε6)))))));
        }
        return m;
    }

    /**
     * Computes the spherical arc length <var>σ</var> from the an ellipsoidal arc length <var>s</var>.
     * This method is the converse of {@link #sphericalToEllipsoidalAngle(double, boolean)}.
     * Input is τ = s/(b⋅A₁).
     *
     * <p><b>Precondition:</b> {@link #computeSeriesExpansionCoefficients()} must have been invoked before this method.</p>
     *
     * @param  τ  arc length on the ellipsoid since equatorial point <var>E</var> divided by {@link #A1}.
     * @return arc length <var>σ</var> on the auxiliary sphere.
     *
     * @see #sphericalToEllipsoidalAngle(double, boolean)
     */
    private double ellipsoidalToSphericalAngle(final double τ) {
        assert !isInvalid(COEFFICIENTS_FOR_START_POINT);
        /*
         * Derived from Karney (2013) equation 21 with θ = 2τ:
         *
         *   ∑C′ₓ⋅sin(x⋅θ)  ≈  (1/2 ⋅ε   +   -9/32  ⋅ε³  +    205/1536 ⋅ε⁵)⋅sin(1θ) +           // C′₁₁
         *                     (5/16⋅ε²  +  -37/96  ⋅ε⁴  +   1335/4096 ⋅ε⁶)⋅sin(2θ) +           // C′₁₂
         *                     (             29/96  ⋅ε³  +    -75/128  ⋅ε⁵)⋅sin(3θ) +           // C′₁₃
         *                     (            539/1536⋅ε⁴  +  -2391/2560 ⋅ε⁶)⋅sin(4θ) +           // C′₁₄
         *                     (                             3467/7680 ⋅ε⁵)⋅sin(5θ) +           // C′₁₅
         *                     (                            38081/61440⋅ε⁶)⋅sin(6θ)             // C′₁₆
         *
         * For performance and simplicity, we rewrite the equations using trigonometric identities
         * using the same technic than the one documented above in sphericalToEllipsoidalAngle(σ).
         */
        final double ε2 = ε*ε;
        final double ε3 = ε*ε2;
        final double ε4 = ε2*ε2;
        final double ε5 = ε2*ε3;
        final double ε6 = ε3*ε3;
        final double θ  = τ * 2;
        final double cosθ = cos(θ);
        return τ + sin(θ) * (   281./240  * ε5  +    -7./12  * ε3  +  1./2 * ε
                 + cosθ   * ( 20753./2560 * ε6  +  -835./384 * ε4  +  5./8 * ε2
                 + cosθ   * ( -4967./640  * ε5  +    29./24  * ε3
                 + cosθ   * (-52427./1920 * ε6  +   539./192 * ε4
                 + cosθ   * (  3467./480  * ε5
                 + cosθ   * ( 38081./1920 * ε6))))));
    }

    /**
     * Computes the longitude angle λ from the equatorial point <var>E</var> to a point on the ellipsoid
     * specified by the ω longitude on auxiliary sphere.
     *
     * <p><b>Precondition:</b> {@link #computeSeriesExpansionCoefficients()} must have been invoked before this method.</p>
     *
     * @param  ω  longitude on the auxiliary sphere.
     * @param  σ  spherical arc length from equatorial point E to point on auxiliary sphere
     *            along a geodesic with azimuth α₀ at equator.
     * @return geodetic longitude λ from the equatorial point E to the point.
     */
    private double sphericalToGeodeticLongitude(final double ω, final double σ) {
        /*
         * Derived from Karney (2013) equation 23:
         *
         *   I₃(σ) = A₃⋅(σ + ∑C₃ₓ⋅sin(x⋅θ))
         *
         * but with the sum of sine terms replaced by Clenshaw summation.
         */
        final double θ    = σ * 2;
        final double cosθ = cos(θ);
        final double I3   = A3*(sin(θ)*(C31 + cosθ*(C32 + cosθ*(C33 + cosθ*(C34 + cosθ*C35)))) + σ);
        if (STORE_LOCAL_VARIABLES) store("I₃(σ)", I3);
        return ω - sinα0 * I3 * (1 - axisRatio);
    }

    /**
     * Sets the {@link #sinα0} and {@link #cosα0} terms. Note that the {@link #msinα2} field is always set
     * to the {@code sinα0} value in this class. But those two fields may nevertheless have different values
     * if {@code msinα2} has been set independently, for example by {@link #createGeodesicPath2D(double)}.
     */
    private void α0(final double sinα1, final double cosα1, final double sinβ1, final double cosβ1) {
        sinα0 = sinα1 * cosβ1;                  // Azimuth at equator (Karney 5) as geographic angle.
        cosα0 = hypot(cosα1, sinα1*sinβ1);      // = sqrt(1 - sinα0*sinα0) with less rounding errors.
    }

    /**
     * Computes the end point from the start point, the azimuth and the geodesic distance.
     * This method should be invoked if the end point or ending azimuth is requested while
     * {@link #END_POINT} validity flag is not set.
     *
     * <p>This implementation computes {@link #φ2}, {@link #λ2} and azimuths using the method
     * described in Karney (2013) <cite>Algorithms for geodesics</cite>. The coefficients of
     * Fourier and Taylor series expansions are given by equations 15 and 17.</p>
     *
     * @throws IllegalStateException if the start point, azimuth or distance has not been set.
     */
    @Override
    final void computeEndPoint() {
        canComputeEndPoint();                                       // May throw IllegalStateException.
        if (isInvalid(COEFFICIENTS_FOR_START_POINT)) {
            final double m     = hypot(msinα1, mcosα1);
            final double sinα1 = msinα1 / m;                        // α is a geographic (not arithmetic) angle.
            final double cosα1 = mcosα1 / m;
            final double tanβ1 = axisRatio * tan(φ1);               // β₁ is reduced latitude (Karney 6).
            final double cosβ1 = 1 / sqrt(1 + tanβ1*tanβ1);
            final double sinβ1 = tanβ1 * cosβ1;
            α0(sinα1, cosα1, sinβ1, cosβ1);
            /*
             * Note:  Karney said that for equatorial geodesics (φ₁=0 and α₁=π/2), calculation of σ₁ is indeterminate
             * and σ₁=0 should be taken. The indetermination appears as atan2(0,0). However this expression evaluates
             * to 0 in Java, as required by Math.atan2(y,x) specification, which is what we want. So there is no need
             * for a special case.
             */
            final double cosα1_cosβ1 = cosα1*cosβ1;
            final double σ1 = atan2(sinβ1,       cosα1_cosβ1);      // Arc length on great circle (Karney 11).
            final double ω1 = atan2(sinβ1*sinα0, cosα1_cosβ1);
            final double k2 = computeSeriesExpansionCoefficients();
            λ1E   = sphericalToGeodeticLongitude(ω1, σ1);
            I1_σ1 = sphericalToEllipsoidalAngle(σ1, false);
            setValid(COEFFICIENTS_FOR_START_POINT);
            if (STORE_LOCAL_VARIABLES) {                            // For comparing values with Karney table 2.
                store("β₁", atan2(sinβ1, cosβ1));
                store("σ₁", σ1);
                store("ω₁", ω1);
                store("k²", k2);
                snapshot();
            }
        }
        /*
         * Distance from equatorial point E to ending point P₂ along the geodesic: s₂ = s₁ + ∆s.
         */
        final double s2b = I1_σ1 + geodesicDistance / semiMinorAxis();      // (Karney 18) + ∆s/b
        final double σ2  = ellipsoidalToSphericalAngle(s2b / A1);           // (Karney 21)
        final double sinσ2 = sin(σ2);
        final double cosσ2 = cos(σ2);
        /*
         * Azimuth at ending point α₂ = atan2(sinα₀, cosα₀⋅cosσ₂)    (Karney 14 adapted for ending point).
         * We do not need atan2(y,x) because we keep x and y components separated. It is not necessary to
         * normalize to a vector of length 1.
         */
        msinα2 = sinα0;
        mcosα2 = cosα0 * cosσ2;
        /*
         * Ending point coordinates on auxiliary sphere: Latitude β is given by Karney equation 13:
         *
         *   β₂ = atan2(cos(α₀)⋅sin(σ₂), hypot(cos(α₀)⋅cos(σ₂), sin(α₀))
         *
         * We replace cos(α₀)⋅cos(σ₂) by mcosα2 since we computed it above. Then we avoid the call to
         * atan2(y,x) by storing directly the y and x values. Note that `sinβ2` and `cosβ2` are not
         * really sine and cosine since we do not normalize them to sin² + cos² = 1. We do not need
         * to normalize because we use either a ratio of those 2 quantities or give them to atan2(…).
         */
        final double sinβ2 = cosα0 * sinσ2;
        final double cosβ2 = hypot(msinα2, mcosα2);                             // m⋅sin(α₂) = sin(α₀) in this class.
        final double ω2    = atan2(sinα0*sinσ2, cosσ2);                         // (Karney 12).
        /*
         * Convert reduced longitude ω and latitude β on auxiliary sphere
         * to geodetic longitude λ and latitude φ.
         */
        final double λ2E = sphericalToGeodeticLongitude(ω2, σ2);
        λ2 = IEEEremainder(λ2E - λ1E + λ1, 2*PI);
        φ2 = atan(sinβ2 / (cosβ2 * axisRatio));                                 // (Karney 6).
        setValid(END_POINT | ENDING_AZIMUTH);
        if (STORE_LOCAL_VARIABLES) {                // For comparing values with Karney table 2.
            store("s₂", s2b * semiMinorAxis());
            store("τ₂", s2b / A1);
            store("σ₂", σ2);
            store("α₂", atan2(msinα2, mcosα2));
            store("β₂", atan2(sinβ2, cosβ2));
            store("ω₂", ω2);
            store("λ₂", λ2E);
            store("Δλ", λ2E - λ1E);
        }
    }

    /**
     * Computes the geodesic distance and azimuths from the start point and end point.
     * This method should be invoked if the distance or an azimuth is requested while
     * {@link #STARTING_AZIMUTH}, {@link #ENDING_AZIMUTH} or {@link #GEODESIC_DISTANCE}
     * validity flag is not set.
     *
     * <p>Reminder: given <var>P₁</var> the starting point and <var>E</var> the intersection of the geodesic with equator:</p>
     * <ul>
     *   <li><var>α₁</var> is the azimuth (0° oriented North) of the geodesic from <var>E</var> to <var>P₁</var>.</li>
     *   <li><var>σ₁</var> is the spherical arc length (in radians) between <var>E</var> and <var>P₁</var>.</li>
     *   <li><var>ω₁</var> is the spherical longitude (in radians) on the auxiliary sphere at <var>P₁</var>.
     *       Spherical longitude is the angle formed by the meridian of <var>E</var> and the meridian of <var>P₁</var>.</li>
     * </ul>
     *
     * @throws IllegalStateException if the distance or azimuth has not been set.
     * @throws GeodeticException if an azimuth or the distance can not be computed.
     */
    @Override
    final void computeDistance() {
        canComputeDistance();
        /*
         * The algorithm in this method requires the following canonical configuration:
         *
         *   Negative latitude of starting point:         φ₁ ≦ 0
         *   Ending point latitude smaller in magnitude:  φ₁ ≦ φ₂ ≦ -φ₁
         *   Positive longitude difference:               0 ≦ ∆λ ≦ π
         *   (Consequence of above):                      0 ≦ α₀ ≦ π/2
         *
         * If the given points do not met above conditions, then we need to swap start and end points or to
         * swap coordinate signs. We apply those changes on local variables only, not on the class fields.
         * We will need to apply the converse of those changes in the final results.
         */
        double φ1 = this.φ1;
        double φ2 = this.φ2;
        double Δλ = IEEEremainder(λ2 - λ1, 2*PI);           // In [-π … +π] range.
        final boolean swapPoints = abs(φ2) > abs(φ1);
        if (swapPoints) {
            φ1 = this.φ2;
            φ2 = this.φ1;
            Δλ = -Δλ;
        }
        final boolean inverseLatitudeSigns = φ1 > 0;
        if (inverseLatitudeSigns) {
            φ1 = -φ1;
            φ2 = -φ2;
        }
        final boolean inverseLongitudeSigns = Δλ < 0;
        if (inverseLongitudeSigns) {
            Δλ = -Δλ;
        }
        /*
         * Compute an approximation of the azimuth α₁ at starting point. This estimation will be refined by iteration
         * in the loop later, but that iteration will not converge if the first α₁ estimation is not good enough. The
         * general formula does not give good α₁ initial value for antipodal points, so we need to check for special
         * cases first:
         *
         *   1) Nearly antipodal points with φ₁ = -φ₂.
         *   2) Nearly antipodal points with φ₁ slightly different than -φ₂.
         *   3) Equatorial case: φ₁ = φ₂ = 0 but restricted to ∆λ ≤ (1-f)⋅π.
         *   4) Meridional case: ∆λ = 0 or ∆λ = π (handled by general case in this method).
         */
        if (φ1 > -LATITUDE_THRESHOLD) {                         // Sufficient test because φ₁ ≦ 0 and |φ₂| ≦ φ₁
            /*
             * Points on equator but not nearly anti-podal. The geodesic is an arc on equator and the azimuths
             * are α₁ = α₂ = ±90°. We need this special case because when φ = 0, the general case get sinβ = 0
             * then σ = atan2(sinβ, cosα⋅cosβ) = 0, which result in a distance of 0. I have not yet understood
             * how to use the general formulas in such case. This code is a workaround in the meantime.
             *
             * See https://issues.apache.org/jira/browse/SIS-467
             */
            if (Δλ > axisRatio * PI) {
                // Karney's special case documented before equation 45.
                throw new GeodeticException("Can not compute geodesics for antipodal points on equator.");
            }
            super.computeDistance();
            return;
        }
        /*
         * Reduced latitudes β (Karney 6). Actually we don't need the β angles
         * (except for a special case), but rather their sine and cosine values.
         */
        final double tanβ1 = axisRatio * tan(φ1);
        final double tanβ2 = axisRatio * tan(φ2);
        final double cosβ1 = 1 / sqrt(1 + tanβ1*tanβ1);
        final double cosβ2 = 1 / sqrt(1 + tanβ2*tanβ2);
        final double sinβ1 = tanβ1 * cosβ1;
        final double sinβ2 = tanβ2 * cosβ2;
        double α1;
        if (Δλ >= PI - NEARLY_ANTIPODAL_Δλ && abs(φ1 + φ2) <= NEARLY_ANTIPODAL_Δλ) {
            /*
             * Nearly antipodal points. Karney's equations 53 are reproduced on the left side
             * (with f = 1 - b/a), which we replace by the equations on the right side below:
             *
             *   Δ  = f⋅a⋅π⋅cos²β₁              ┃      Δ₁ = f⋅π⋅cosβ₁
             *   ∆λ = π + Δ/(a⋅cosβ₁)⋅x         ┃      ∆λ = π + Δ₁⋅x
             *   β₂ = -β₁ + (Δ/a)⋅y             ┃      β₂ = -β₁ + (Δ₁⋅cosβ₁)⋅y
             */
            final double β1 = atan2(sinβ1, cosβ1);
            final double β2 = atan2(sinβ2, cosβ2);
            final double Δ1  = (1 - axisRatio) * PI * cosβ1;                // Differ from Karney by a⋅cosβ₁ factor.
            final double y  = (β2 + β1) / (Δ1*cosβ1);
            final double x  = (PI - Δλ) / Δ1;                               // Opposite sign of Karney. We have x ≧ 0.
            final double x2 = x*x;
            final double y2 = y*y;
            if (STORE_LOCAL_VARIABLES) {                        // For comparing values with Karney table 4.
                store("x", -x);                                 // Sign used by Karney.
                store("y",  y);
            }
            if (y2 < 1E-12) {                                   // Empirical threshold. See μ(…) for more information.
                α1 = (x2 > 1) ? PI/2 : atan(x / sqrt(1 - x2));  // (Karney 57) with opposite sign of x. Result in α₁ ≧ 0.
            } else {
                final double μ = μ(x2, y2);
                α1 = atan2(x*μ, (y*(1+μ)));                     // (Karney 56) with opposite sign of x.
                if (STORE_LOCAL_VARIABLES) {
                    store("μ", μ);
                }
            }
        } else {
            /*
             * Usual case (non-antipodal). Estimation is based on variation of geodetic longitude λ
             * and spherical longitude ω on the auxiliary sphere. Karney makes a special case for
             * Δλ = 0 and Δλ = π by defining α₁ = Δλ. However formulas below produce the same result.
             */
            double w = (cosβ1 + cosβ2) / 2;
            w = sqrt(1 - eccentricitySquared * (w*w));
            final double Δω = Δλ / w;
            final double cω = cos(Δω);
            final double sω = sin(Δω);
            final double αx = cosβ1*sinβ2 - sinβ1*cosβ2*cω;
            final double αy = cosβ2*sω;
            α1 = atan2(αy, αx);                                                                 // (Karney 49)
            if (STORE_LOCAL_VARIABLES) {        // For comparing values with Karney table 3.
                store("ωb", w);
                store("Δω", Δω);
                store("α₂", atan2(cosβ1*sω,      cosβ1*cω*sinβ2 - sinβ1*cosβ2));                // (Karney 50)
                store("Δσ", atan2(hypot(αy, αx), cosβ1*cω*cosβ2 + sinβ1*sinβ2));                // (Karney 51)
                // For small distances we could stop here. But it is easier to let iteration does its work.
            }
        }
        /*
         * Stores locale variables for comparison against Karney tables 4, 5 and 6. Values β₁ and β₂ are keep
         * constant during all this method. Value α₁ is a first estimation and will be updated during iteration.
         * Not that because Karney separate calculation of α₁ and remaining calculation in two separated tables,
         * we need to truncate α₁ to the same number of digits than Karney in order to get the same numbers in
         * the rest of this method.
         */
        if (STORE_LOCAL_VARIABLES) {
            store("β₁", atan(tanβ1));
            store("β₂", atan(tanβ2));
            store("α₁", α1);
            α1 = computedToGiven(α1);
        }
        /*
         * Refine α₁ using Newton's method. Karney proposes an hybrid approach: approximate α₁,
         * then compute σ₁, σ₂, ω₁ and ω₂ (actually the sine and cosine of those angles in our
         * implementation).
         */
        double σ1, σ2;
        int moreRefinements = Formulas.MAXIMUM_ITERATIONS;
        do {
            α0(msinα1 = sin(α1), mcosα1 = cos(α1), sinβ1, cosβ1);
            final double k2 = computeSeriesExpansionCoefficients();
            /*
             * Compute α₂ from Karney equation 5: sin(α₀) = sin(α₁)⋅cos(β₁) = sin(α₂)⋅cos(β₂)
             * The cos(α₂) term could be computed from sin(α₂), but Karnay recommends instead
             * equation 45. An older publication (Karney 2010) went one step further with the
             * replacement of (cos²β₂ - cos²β₁) by:
             *
             *     (cosβ₂ - cosβ₁)⋅(cosβ₂ + cosβ₁)  if β₁ < -π/4        (Reminder: β₁ is always negative)
             *     (sinβ₁ - sinβ₂)⋅(sinβ₁ + sinβ₂)  otherwise.
             *
             * Actually we don't need α values directly, but instead cos(α)⋅cos(β).
             * Note that cos(α₁) can be negative because α₁ ∈ [0…2π].
             */
            final double cosα1_cosβ1 = mcosα1 * cosβ1;
            final double cosα2_cosβ2 = sqrt(cosα1_cosβ1*cosα1_cosβ1 +
                    (cosβ1 <= 1/MathFunctions.SQRT_2 ? (cosβ2 - cosβ1)*(cosβ2 + cosβ1)
                                                     : (sinβ1 - sinβ2)*(sinβ1 + sinβ2)));
            msinα2 = sinα0;
            mcosα2 = cosα2_cosβ2;
            /*
             * Karney gives the following formulas:
             *
             *   σ  =  atan2(sinβ, cosα⋅cosβ)        —     spherical arc length.
             *   ω  =  atan2(sin(α₀)⋅sinσ, cosσ)     —     spherical longitude.
             *
             * We perform the following substitutions where c is an unknown constant:
             * That unknown coefficient will disaspear in atan2(c⋅y, c⋅x) expressions:
             *
             *   sin(σ) = c⋅sinβ
             *   cos(σ) = c⋅cosα⋅cosβ
             */
            final double ω1, ω2;
            σ1 = atan2(sinβ1,       cosα1_cosβ1);               // (Karney 11)
            σ2 = atan2(sinβ2,       cosα2_cosβ2);
            ω1 = atan2(sinβ1*sinα0, cosα1_cosβ1);               // (Karney 12) with above-cited substitutions.
            ω2 = atan2(sinβ2*sinα0, cosα2_cosβ2);
            /*
             * Compute the difference in longitude using the current start point and end point.
             * If this difference is close enough to the requested accuracy, we are almost done.
             * Otherwise refine the results. Note that if we detect that accuracy is good enough,
             * we still complete the computation in order to not waste what we have computed so
             * far in current iteration.
             */
            final double λ2E;
            λ1E = sphericalToGeodeticLongitude(ω1, σ1);
            λ2E = sphericalToGeodeticLongitude(ω2, σ2);
            final double Δλ_error = IEEEremainder(λ2E - λ1E - Δλ, 2*PI);
            if (abs(Δλ_error) <= ITERATION_TOLERANCE) {
                moreRefinements = 0;
            } else if (--moreRefinements == 0) {
                throw new GeodeticException(Resources.format(Resources.Keys.NoConvergence));
            }
            /*
             * Special case for α₁ = π/2 and β₂ = ±β₁ (Karney's equation 47). We replace the β₂ = ±β₁
             * condition by |β₂| - |β₁| ≈ 0. Assuming tan(θ) ≈ θ for small angles we take the tangent
             * of above difference and use tan(β₂ - β₁) = (tanβ₂ - tanβ₁)/(1 + tanβ₂⋅tanβ₁) identity.
             * Note that tanβ₁ ≦ 0 and |tanβ₂| ≦ |tanβ₁| in this method.
             */
            final double dΔλ_dα1;
            if (abs(mcosα1) < LATITUDE_THRESHOLD && (-tanβ1 - abs(tanβ2)) < (1 + abs(tanβ1*tanβ2)) * LATITUDE_THRESHOLD) {
                dΔλ_dα1 = -2 * sqrt(1 - eccentricitySquared * (cosβ1*cosβ1)) / sinβ1;
            } else {
                /*
                 * Karney's equation 38 combined with equation 46. The substitutions
                 * for sin(σ) and cos(σ) described above are applied again here.
                 */
                final double h1 = hypot(sinβ1, cosα1_cosβ1);
                final double h2 = hypot(sinβ2, cosα2_cosβ2);
                final double sinσ1 = sinβ1 / h1;
                final double sinσ2 = sinβ2 / h2;
                final double cosσ1 = cosα1_cosβ1 / h1;
                final double cosσ2 = cosα2_cosβ2 / h2;
                final double J2 = sphericalToEllipsoidalAngle(σ2, true);
                final double J1 = sphericalToEllipsoidalAngle(σ1, true);
                final double Δm = (sqrt(1 + k2*(sinσ2*sinσ2))*cosσ1*sinσ2
                                - sqrt(1 + k2*(sinσ1*sinσ1))*sinσ1*cosσ2
                                - cosσ1*cosσ2*(J2 - J1));
                dΔλ_dα1 = Δm * axisRatio / cosα2_cosβ2;
                if (STORE_LOCAL_VARIABLES) {
                    store("J(σ₁)", J1);
                    store("J(σ₂)", J2);
                    store("Δm",    Δm * semiMinorAxis());
                }
            }
            final double dα1 = Δλ_error / dΔλ_dα1;                  // Opposite sign of Karney δα₁ term.
            α1 -= dα1;
            if (STORE_LOCAL_VARIABLES) {                            // For comparing values against Karney table 5 and 6.
                final double I1_σ2;
                I1_σ1 = sphericalToEllipsoidalAngle(σ1, false);     // Required for computation of s₁ in `snapshot()`.
                I1_σ2 = sphericalToEllipsoidalAngle(σ2, false);
                snapshot();
                store("k²",    k2);
                store("σ₁",    σ1);
                store("σ₂",    σ2);
                store("ω₁",    ω1);
                store("ω₂",    ω2);
                store("α₂",    atan2(msinα2, mcosα2));
                store("λ₂",    λ2E);
                store("Δλ",    λ2E - λ1E);
                store("δλ",    Δλ_error);
                store("dλ/dα", dΔλ_dα1);
                store("δσ₁",  -dα1);
                store("α₁",    α1);
                store("s₂",    I1_σ2 * semiMinorAxis());
                store("Δs",    (I1_σ2 - I1_σ1) * semiMinorAxis());
            }
        } while (moreRefinements != 0);
        final double I1_σ2;
        I1_σ1 = sphericalToEllipsoidalAngle(σ1, false);
        I1_σ2 = sphericalToEllipsoidalAngle(σ2, false);
        geodesicDistance = (I1_σ2 - I1_σ1) * semiMinorAxis();
        /*
         * Restore the coordinate sign and order to the original configuration.
         */
        if (swapPoints) {
            double t;
            t = msinα1; msinα1 = msinα2; msinα2 = t;
            t = mcosα1; mcosα1 = mcosα2; mcosα2 = t;
        }
        if (inverseLongitudeSigns ^ swapPoints) {
            msinα1 = -msinα1;
            msinα2 = -msinα2;
        }
        if (inverseLatitudeSigns ^ swapPoints) {
            mcosα1 = -mcosα1;
            mcosα2 = -mcosα2;
        }
        setValid(STARTING_AZIMUTH | ENDING_AZIMUTH | GEODESIC_DISTANCE);
        if (!(swapPoints | inverseLongitudeSigns | inverseLatitudeSigns)) {
            setValid(COEFFICIENTS_FOR_START_POINT);
        }
    }

    /**
     * Computes the positive root of quartic equation for estimation of α₁ in nearly antipodal case.
     * Formula is given in appendix B of <a href="https://arxiv.org/pdf/1102.1215.pdf">C.F.F Karney (2011)</a>
     * given <var>x</var> and <var>y</var> the coordinates on a plane coordinate system centered on the antipodal point:
     *
     * {@preformat math
     *     μ⁴ + 2μ³ + (1−x²-y²)μ² − 2y²μ - y² = 0
     * }
     *
     * The results should have only one positive root {@literal (μ > 0)}.
     *
     * <div class="section">Condition on <var>y</var> value</div>
     * This method is indeterminate when <var>y</var> → 0 (it returns {@link Double#NaN}). For values too close to zero,
     * the result may be non-significative because of rounding errors. For choosing a threshold value for <var>y</var>,
     * {@code GeodesicsOnEllipsoidTest.Calculator} compares the value computed by this method against the value computed
     * by {@link org.apache.sis.math.MathFunctions#polynomialRoots(double...)}. If the values differ too much, we presume
     * that they are mostly noise caused by a <var>y</var> value too low.
     *
     * @param  x2  the square of <var>x</var>.
     * @param  y2  the square of <var>y</var>.
     */
    private static double μ(final double x2, final double y2) {
        final double r  = (x2 + y2 - 1)/6;
        final double r3 = r*r*r;
        final double S  = x2*y2/4;
        final double d  = S*(S + 2*r3);
        final double u;
        if (d < 0) {
            u = r*(1 + 2*cos(atan2(sqrt(-d), -(S + r3))/3));
        } else {
            final double T = cbrt(S + r3 + copySign(sqrt(d), S + r3));
            u = (T == 0) ? 0 : r + T + r*r/T;
        }
        final double v = sqrt(u*u + y2);
        final double w = (v + u - y2) / (2*v);
        return (v + u) / (sqrt(v + u + w*w) + w);
    }

    /**
     * Computes (∂y/∂φ)⁻¹ at the given latitude on an ellipsoid with semi-major axis length of 1.
     * This derivative is close to cos(φ) for a slightly flattened sphere.
     *
     * @see org.apache.sis.referencing.operation.projection.ConformalProjection#dy_dφ
     */
    @Override
    double dφ_dy(final double φ) {
        final double sinφ = sin(φ);
        final double cosφ = cos(φ);
        return cosφ/(1  -  eccentricitySquared * (cosφ*cosφ) / (1 - eccentricitySquared * (sinφ*sinφ)));
    }

    /**
     * Takes a snapshot of the current fields in this class. This is used for JUnit tests only. During development phases,
     * {@link #STORE_LOCAL_VARIABLES} should be {@code true} for allowing {@code GeodesicsOnEllipsoidTest} to verify the
     * values of a large range of local variables. But when the storage of locale variables is disabled, this method allows
     * {@code GeodesicsOnEllipsoidTest} to still verify a few variables.
     */
    final void snapshot() {
        store("ε",      ε);
        store("A₁",     A1);
        store("A₂",     A2);
        store("A₃",     A3);
        store("α₀",     atan2(sinα0, cosα0));
        store("I₁(σ₁)", I1_σ1);
        store("s₁",     I1_σ1 * semiMinorAxis());
        store("λ₁",     λ1E);
    }

    /**
     * Stores the value of a local variable or a field. This is used in JUnit tests only.
     *
     * @param  name   name of the local variable.
     * @param  value  value of the local variable.
     */
    void store(String name, double value) {
    }

    /**
     * Replaces a computed value by the value given in Karney table. This is used when the result published in Karney
     * table 4 is used as input for Karney table 5. We need to truncate the value to the same numbers of digits than
     * Karney, otherwise computation results will differ.
     */
    double computedToGiven(double α1) {
        return α1;
    }

    /**
     * Computes rhumb line using series expansion.
     *
     * <p><b>Source:</b> G.G. Bennett, 1996. <a href="https://doi.org/10.1017/S0373463300013151">
     * Practical Rhumb Line Calculations on the Spheroid</a>. J. Navigation 49(1), 112-119.</p>
     */
    @Override
    final void computeRhumbLine() {
        canComputeDistance();
        /*
         * Bennett (1996) equation 1 computes isometric latitudes Ψ for given geodetic latitudes φ:
         *
         *     Ψ(φ) = log(tan(PI/4 + φ/2) * pow((1 - ℯsinφ) / (1 + ℯsinφ), ℯ/2));
         *
         * (ℯ is the eccentricity, not squared). However we need only the isometric latitudes difference:
         *
         *     ΔΨ = Ψ(φ₂) - Ψ(φ₁)
         *
         * We rewrite the equation using log(Ψ₁) - log(Ψ₂) = log(Ψ₁/Ψ₂) and other identities:
         *
         *     ΔΨ = log(tan(PI/4 + φ₂/2)) + log(pow((1 - ℯsinφ₂) / (1 + ℯsinφ₂), ℯ/2))
         *        - log(tan(PI/4 + φ₁/2)) - log(pow((1 - ℯsinφ₁) / (1 + ℯsinφ₁), ℯ/2))
         *
         *        = log(tan(PI/4 + φ₂/2) /
         *              tan(PI/4 + φ₁/2) *
         *              pow(((1 - ℯsinφ₂)*(1 + ℯsinφ₁)) /
         *                  ((1 + ℯsinφ₂)*(1 - ℯsinφ₁)), ℯ/2))
         *
         * The code below combines ℯsinφ terms otherwise. Note that we could also use product-to-sum
         * identities for rewriting the  tan(¼π + ½φ₂) / tan(¼π + ½φ₁)  expression as  (a + b) / (a - b)
         * where  a = cos((φ₂+φ₁)/2)  and  b = sin((φ₂-φ₁)/2), but the amount of trigonometric method calls
         * would be about the same and result may be less accurate.
         */
        final double eccentricity = sqrt(eccentricitySquared);      // TODO: avoid computing on each invocation.
        final double sinφ1 = sin(φ1);
        final double sinφ2 = sin(φ2);
        final double sd = eccentricity * (sinφ1 - sinφ2);
        final double sm = 1 - eccentricitySquared * (sinφ1 * sinφ2);
        final double ΔΨ = log(tan(PI/4 + φ2/2) / tan(PI/4 + φ1/2) * pow((sm+sd)/(sm-sd), eccentricity/2));
        final double Δλ = IEEEremainder(λ2 - λ1, 2*PI);
        final double h  = hypot(Δλ, ΔΨ);
        final double S;
        if (abs(φ1 - φ2) < LATITUDE_THRESHOLD) {
            final double φm = (φ1 + φ2)/2;
            final double sinφ = sin(φm);
            S = cos(φm) / sqrt(1 - eccentricitySquared*(sinφ*sinφ));        // Bennett equation 4 with sin(α) = Δλ/h.
        } else {
            final double m1 = m(φ1, sinφ1);
            final double m2 = m(φ2, sinφ2);
            S = (m2 - m1) / ΔΨ;                                             // Bennett (1996) with cos(α) = ΔΨ/h.
            if (STORE_LOCAL_VARIABLES) {
                store("m₁", m1);
                store("m₂", m2);
                store("Δm", m2 - m1);
            }
        }
        rhumblineLength = S * h * semiMajorAxis;
        rhumblineAzimuth = atan2(Δλ, ΔΨ);
        if (STORE_LOCAL_VARIABLES) {
            store("Δλ", Δλ);
            store("ΔΨ", ΔΨ);
        }
    }

    /**
     * Computes Bennett (1996) equation 2 modified with Clenshaw summation.
     */
    private double m(final double φ, final double sinφ) {
        final double cosφ = cos(φ);
        final double sinθ = 2*sinφ*cosφ;                        // sin(2φ)
        final double cosθ = (cosφ + sinφ) * (cosφ - sinφ);      // cos(2φ)  =  cos²φ - sin²φ
        return R0*φ + sinθ*(R2 + cosθ*(R4 + cosθ*R6));
    }
}