sis-commits mailing list archives

Site index · List index
Message view « Date » · « Thread »
Top « Date » · « Thread »
From desruisse...@apache.org
Subject [sis] 01/01: Creation of GeodesicsOnEllipsoid class for computing geodescis on ellipsoid of revolution. Based on the algorithm published in Karney C.F.F. (2013), Algorithms for geodesics, SRI International, https://doi.org/10.1007/s00190-012-0578-z Deprecate DefaultEllipsoid.orthodromicDistance(…), which is replaced by GeodeticCalculator. https://issues.apache.org/jira/browse/SIS-386
Date Sun, 28 Jul 2019 07:34:17 GMT
This is an automated email from the ASF dual-hosted git repository.

desruisseaux pushed a commit to branch geoapi-4.0
in repository https://gitbox.apache.org/repos/asf/sis.git

commit 55e1313451379b1f030e4d1a2087dfa87bccec8a
Author: Matthieu_Bastianelli <matthieu.bastianelli@geomatys.com>
AuthorDate: Tue Jun 11 19:04:21 2019 +0200

    Creation of GeodesicsOnEllipsoid class for computing geodescis on ellipsoid of revolution.
    Based on the algorithm published in Karney C.F.F. (2013), Algorithms for geodesics, SRI International, https://doi.org/10.1007/s00190-012-0578-z
    Deprecate DefaultEllipsoid.orthodromicDistance(…), which is replaced by GeodeticCalculator.
    https://issues.apache.org/jira/browse/SIS-386
---
 .../apache/sis/internal/referencing/Formulas.java  |   5 +-
 .../sis/referencing/GeodesicsOnEllipsoid.java      | 913 +++++++++++++++++++++
 .../apache/sis/referencing/GeodeticCalculator.java | 174 ++--
 .../sis/referencing/datum/DefaultEllipsoid.java    |   7 +-
 .../org/apache/sis/referencing/datum/Sphere.java   |   1 +
 .../referencing/operation/GeodesicException.java   |  56 ++
 .../operation/projection/Mollweide.java            |   4 +-
 .../sis/referencing/GeodesicsOnEllipsoidTest.java  | 514 ++++++++++++
 .../sis/referencing/GeodeticCalculatorTest.java    | 312 ++++---
 .../sis/test/suite/ReferencingTestSuite.java       |   1 +
 .../java/org/apache/sis/math/MathFunctions.java    | 249 +++++-
 .../org/apache/sis/math/MathFunctionsTest.java     |  43 +
 .../test/java/org/apache/sis/test/TestCase.java    |   3 +
 .../java/org/apache/sis/index/tree/QuadTree.java   |  12 +-
 14 files changed, 2106 insertions(+), 188 deletions(-)

diff --git a/core/sis-referencing/src/main/java/org/apache/sis/internal/referencing/Formulas.java b/core/sis-referencing/src/main/java/org/apache/sis/internal/referencing/Formulas.java
index fe86928..b3b6afe 100644
--- a/core/sis-referencing/src/main/java/org/apache/sis/internal/referencing/Formulas.java
+++ b/core/sis-referencing/src/main/java/org/apache/sis/internal/referencing/Formulas.java
@@ -88,8 +88,11 @@ public final class Formulas extends Static {
      * Maximum number of iterations for iterative computations. Defined in this {@code Formulas} class as a default value,
      * but some classes may use a derived value (for example twice this amount). This constant is mostly useful for identifying
      * places where iterations occur.
+     *
+     * <p>In <cite>Algorithms for geodesics</cite>, Karney (2013) reports that in a tiny fraction of cases up to 16 iterations
+     * is required for the Newton method published in his article. We set {@code MAXIMUM_ITERATIONS} to that maximal value.</p>
      */
-    public static final int MAXIMUM_ITERATIONS = 15;
+    public static final int MAXIMUM_ITERATIONS = 16;
 
     /**
      * Do not allow instantiation of this class.
diff --git a/core/sis-referencing/src/main/java/org/apache/sis/referencing/GeodesicsOnEllipsoid.java b/core/sis-referencing/src/main/java/org/apache/sis/referencing/GeodesicsOnEllipsoid.java
new file mode 100644
index 0000000..7ccf8ef
--- /dev/null
+++ b/core/sis-referencing/src/main/java/org/apache/sis/referencing/GeodesicsOnEllipsoid.java
@@ -0,0 +1,913 @@
+/*
+ * 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.referencing.operation.GeodesicException;
+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.*;
+import static org.apache.sis.internal.metadata.ReferencingServices.NAUTICAL_MILE;
+
+
+/**
+ * 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.
+ * 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 value should be {@code false} in releases.
+     *
+     * @see #snapshot()
+     */
+    static final boolean STORE_LOCAL_VARIABLES = true;
+
+    /**
+     * 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);
+
+    /**
+     * Maximal difference (in radians) between two latitudes for enabling the use of simplified formulas.
+     * This is used in two contexts:
+     *
+     * <ul>
+     *   <li>Maximal difference between latitude φ₁ and equator for using the equatorial approximation.</li>
+     *   <li>Maximal difference between |β₁| and |β₂| for enabling the use of Karney's equation 47.</li>
+     * </ul>
+     *
+     * Those special cases are needed when general formulas produce indeterminations like 0/0.
+     * Current angular value corresponds to a distance of 1 millimetre on a planet of the size
+     * of Earth, which is about 1.57E-10 radians. This value is chosen empirically by trying to
+     * minimize the amount of "No convergence errors" reported by {@code GeodesicsOnEllipsoidTest}
+     * in the {@code compareAgainstDataset()} method.
+     *
+     * <p><b>Note:</b> this is an angular tolerance threshold, but is also used with sine and cosine values
+     * because sin(θ) ≈ θ for small angles.</p>
+     */
+    private static final double LATITUDE_THRESHOLD = 0.001 / (NAUTICAL_MILE*60) * (PI/180);
+
+    /**
+     * Desired accuracy 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.
+     *
+     * <p><b>Note:</b> when the iteration loop detects that it reached this accuracy, the loop completes the
+     * iteration step which was in progress. Consequently the final accuracy is one iteration better than this
+     * accuracy.</p>
+     */
+    private static final double ACCURACY = Formulas.ANGULAR_TOLERANCE * (PI/180) / 20;
+
+    /**
+     * 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 final double semiMinor;
+
+    /**
+     * 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.
+     *
+     * @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.
+     *
+     * @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).
+     */
+    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}.
+     *
+     * <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;
+
+    /**
+     * 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 = ellipsoid.getSemiMajorAxis();
+        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;
+        semiMinor = b;
+    }
+
+    /**
+     * 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, …
+         *
+         * TODO: combine using the technic described in sphericalToEllipsoidalAngle,
+         * then make constant positive.
+         */
+        double fε;
+        C31 = (fε  = ε) * ( 1./4     +   n*(-1./4)
+                   + ε  * ( 1./8     +   n*(             n*(-1./8 ))
+                   + ε  * ( 3./64    +   n*(3./64    +   n*(-1./64))
+                   + ε  * ( 5./128   +   n*(1./64)
+                   + ε  * ( 3./128)))));
+        C32 = (fε *= ε) * ( 1./16    +   n*(-3./32   +   n*( 1./32))
+                   + ε  * ( 3./64    +   n*(-1./32   +   n*(-3./64))
+                   + ε  * ( 3./128   +   n*( 1./128)
+                   + ε  * ( 5./256))));
+        C33 = (fε *= ε) * ( 5./192   +   n*(-3./64   +   n*( 5./192))
+                   + ε  * ( 3./128   +   n*(-5./192)
+                   + ε  * ( 7./512)));
+        C34 = (fε *= ε) * ( 7./512   +   n*(-7./256)
+                   + ε  * ( 7./512));
+        C35 = (fε *  ε) * (21./2560);
+        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.
+     *
+     * @todo combine using the technic described in sphericalToEllipsoidalAngle,
+     */
+    private double sphericalToGeodeticLongitude(final double ω, final double σ) {
+        /*
+         * Derived from Karney (2013) equation 23.
+         *
+         *   I₃(σ) = A₃⋅(σ + ∑C₃ₓ⋅sin(x⋅θ))
+         */
+        final double I3 = A3*(C35 * sin(2 * 5 * σ)
+                            + C34 * sin(2 * 4 * σ)
+                            + C33 * sin(2 * 3 * σ)
+                            + C32 * sin(2 * 2 * σ)
+                            + C31 * sin(2 * 1 * σ) + σ);
+        if (STORE_LOCAL_VARIABLES) store("I₃(σ)", I3);
+        return ω - sinα0 * I3 * (1 - axisRatio);
+    }
+
+    /**
+     * 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 vm    = hypot(dφ1, dλ1);
+            final double sinα1 = dλ1 / vm;                          // α is a geographic (not arithmetic) angle.
+            final double cosα1 = dφ1 / vm;
+            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;
+            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.
+            /*
+             * 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 / semiMinor;            // (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.
+         */
+        dλ2 = sinα0;                                // cos(π/2 − θ)  =  sin(θ)
+        dφ2 = cosα0 * cosσ2;                        // sin(π/2 − θ)  =  cos(θ)
+        /*
+         * 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 dφ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(dφ2, sinα0);
+        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 * semiMinor);
+            store("τ₂", s2b / A1);
+            store("σ₂", σ2);
+            store("α₂", atan2(sinα0, cosα0*cosσ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 GeodesicException if an azimuth or the distance can not be computed.
+     */
+    @Override
+    final void computeDistance() throws GeodesicException {
+        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 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;
+        /*
+         * 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 GeodesicException("Can not compute geodesics for antipodal points on equator.");
+            }
+            final double Δφ = φ2 - φ1;
+            geodesicDistance = hypot(Δλ, Δφ) * ellipsoid.getSemiMajorAxis();
+            dλ1 = dλ2 = (inverseLongitudeSigns ^ swapPoints) ? -Δλ : Δλ;
+            dφ1 = dφ2 = (inverseLatitudeSigns  ^ swapPoints) ? -Δφ : Δφ;
+            setValid(STARTING_AZIMUTH | ENDING_AZIMUTH | GEODESIC_DISTANCE);
+            return;
+        }
+        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;
+        for (int nbIter = Formulas.MAXIMUM_ITERATIONS;;) {
+            final double sinα1 = sin(α1);
+            final double cosα1 = cos(α1);
+            sinα0 = sinα1 * cosβ1;
+            cosα0 = hypot(cosα1, sinα1*sinβ1);                  // = sqrt(1 - sinα0*sinα0) with less rounding errors.
+            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 = cosα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)));
+            /*
+             * 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);
+            final boolean done = (abs(Δλ_error) <= ACCURACY);
+            if (--nbIter < 0 && !done) {
+                throw new GeodesicException(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(cosα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 * semiMinor);
+                }
+            }
+            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(sinα0, cosα2_cosβ2));
+                store("λ₂",    λ2E);
+                store("Δλ",    λ2E - λ1E);
+                store("δλ",    Δλ_error);
+                store("dλ/dα", dΔλ_dα1);
+                store("δσ₁",  -dα1);
+                store("α₁",    α1);
+                store("s₂",    I1_σ2 * semiMinor);
+                store("Δs",    (I1_σ2 - I1_σ1) * semiMinor);
+            }
+            if (done) {
+                dλ1 = sinα1;
+                dφ1 = cosα1;
+                dλ2 = sinα0;
+                dφ2 = cosα2_cosβ2;
+                break;
+            }
+        }
+        final double I1_σ2;
+        I1_σ1 = sphericalToEllipsoidalAngle(σ1, false);
+        I1_σ2 = sphericalToEllipsoidalAngle(σ2, false);
+        geodesicDistance = (I1_σ2 - I1_σ1) * semiMinor;
+        /*
+         * Restore the coordinate sign and order to the original configuration.
+         */
+        if (swapPoints) {
+            double t;
+            t = dφ1; dφ1 = dφ2; dφ2 = t;
+            t = dλ1; dλ1 = dλ2; dλ2 = t;
+        }
+        if (inverseLongitudeSigns ^ swapPoints) {
+            dλ1 = -dλ1;
+            dλ2 = -dλ2;
+        }
+        if (inverseLatitudeSigns ^ swapPoints) {
+            dφ1 = -dφ1;
+            dφ2 = -dφ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);
+    }
+
+    /**
+     * 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 * semiMinor);
+        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;
+    }
+}
diff --git a/core/sis-referencing/src/main/java/org/apache/sis/referencing/GeodeticCalculator.java b/core/sis-referencing/src/main/java/org/apache/sis/referencing/GeodeticCalculator.java
index 7687d69..f99dd1c 100644
--- a/core/sis-referencing/src/main/java/org/apache/sis/referencing/GeodeticCalculator.java
+++ b/core/sis-referencing/src/main/java/org/apache/sis/referencing/GeodeticCalculator.java
@@ -37,6 +37,7 @@ import org.apache.sis.measure.AngleFormat;
 import org.apache.sis.measure.Latitude;
 import org.apache.sis.measure.Units;
 import org.apache.sis.geometry.CoordinateFormat;
+import org.apache.sis.referencing.operation.GeodesicException;
 import org.apache.sis.internal.referencing.PositionTransformer;
 import org.apache.sis.internal.referencing.ReferencingUtilities;
 import org.apache.sis.internal.referencing.j2d.ShapeUtilities;
@@ -74,6 +75,7 @@ import static java.lang.Math.*;
  * This class is not thread-safe. If geodetic calculations are needed in a multi-threads environment,
  * then a distinct instance of {@code GeodeticCalculator} needs to be created for each thread.
  *
+ * @author  Martin Desruisseaux (Geomatys)
  * @version 1.0
  *
  * @see <a href="https://en.wikipedia.org/wiki/Great-circle_navigation">Great-circle navigation on Wikipedia</a>
@@ -108,7 +110,7 @@ public class GeodeticCalculator {
      *
      * @see org.apache.sis.referencing.datum.DefaultEllipsoid#getAuthalicRadius()
      */
-    private final double radius;
+    final double authalicRadius;
 
     /**
      * The (<var>latitude</var>, <var>longitude</var>) coordinates of the start point <strong>in radians</strong>.
@@ -116,7 +118,7 @@ public class GeodeticCalculator {
      *
      * @see #START_POINT
      */
-    private double φ1, λ1;
+    double φ1, λ1;
 
     /**
      * The (<var>latitude</var>, <var>longitude</var>) coordinates of the end point <strong>in radians</strong>.
@@ -124,7 +126,7 @@ public class GeodeticCalculator {
      *
      * @see #END_POINT
      */
-    private double φ2, λ2;
+    double φ2, λ2;
 
     /**
      * The azimuth at start point and end point as vector components.
@@ -138,11 +140,17 @@ public class GeodeticCalculator {
      *
      * Those vectors may not be normalized to unitary vectors. This is often not needed because most formulas are
      * written in a way that cancel the magnitude. If nevertheless needed, normalization will be applied in formulas.
+     * Those quantities are related to α as below, with α a <em>geographic</em> (not arithmetic) angle:
+     *
+     * <ul>
+     *   <li>{@code sinα = dλ / hypot(dφ, dλ)}</li>
+     *   <li>{@code cosα = dφ / hypot(dφ, dλ)}</li>
+     * </ul>
      *
      * @see #STARTING_AZIMUTH
      * @see #ENDING_AZIMUTH
      */
-    private double dφ1, dλ1, dφ2, dλ2;
+    double dφ1, dλ1, dφ2, dλ2;
 
     /**
      * The shortest distance from the starting point ({@link #φ1},{@link #λ1}) to the end point ({@link #φ2},{@link #λ2}).
@@ -151,7 +159,7 @@ public class GeodeticCalculator {
      * @see #GEODESIC_DISTANCE
      * @see #getDistanceUnit()
      */
-    private double geodesicDistance;
+    double geodesicDistance;
 
     /**
      * Length of the rhumb line from the starting point ({@link #φ1},{@link #λ1}) to the end point ({@link #φ2},{@link #λ2}).
@@ -160,21 +168,24 @@ public class GeodeticCalculator {
      * @see #RHUMBLINE_LENGTH
      * @see #getDistanceUnit()
      */
-    private double rhumblineLength;
+    double rhumblineLength;
 
     /**
      * A bitmask specifying which information are valid. For example if the {@link #END_POINT} bit is not set,
      * then {@link #φ2} and {@link #λ2} need to be computed, which implies the computation of ∂φ/∂λ as well.
      * If the {@link #GEODESIC_DISTANCE} bit is not set, then {@link #geodesicDistance} needs to be computed,
      * which implies recomputation of ∂φ/∂λ as well.
+     *
+     * @see #isInvalid(int)
+     * @see #setValid(int)
      */
     private int validity;
 
     /**
      * Bitmask specifying which information are valid.
      */
-    private static final int START_POINT = 1, END_POINT = 2, STARTING_AZIMUTH = 4, ENDING_AZIMUTH = 8,
-            GEODESIC_DISTANCE = 16, RHUMBLINE_LENGTH = 32;
+    static final int START_POINT = 1, END_POINT = 2, STARTING_AZIMUTH = 4, ENDING_AZIMUTH = 8,
+            GEODESIC_DISTANCE = 16, RHUMBLINE_LENGTH = 32, COEFFICIENTS_FOR_START_POINT = 64;
 
     /**
      * Constructs a new geodetic calculator expecting coordinates in the supplied CRS.
@@ -186,16 +197,16 @@ public class GeodeticCalculator {
      * <p>This class is currently not designed for sub-classing outside this package. If in a future version we want to
      * relax this restriction, we should revisit the package-private API in order to commit to a safer protected API.</p>
      *
-     * @param  crs  the reference system for the {@link Position} arguments and return values.
+     * @param  crs         the reference system for the {@link Position} arguments and return values.
+     * @param  ellipsoid   ellipsoid associated to the geodetic component of given CRS.
      */
-    GeodeticCalculator(final CoordinateReferenceSystem crs) {
-        ArgumentChecks.ensureNonNull("crs", crs);
+    GeodeticCalculator(final CoordinateReferenceSystem crs, final Ellipsoid ellipsoid) {
         final GeographicCRS geographic = ReferencingUtilities.toNormalizedGeographicCRS(crs, true, true);
         if (geographic == null) {
             throw new IllegalArgumentException(Errors.format(Errors.Keys.IllegalCRSType_1, ReferencingUtilities.getInterface(crs)));
         }
-        ellipsoid      = ReferencingUtilities.getEllipsoid(crs);
-        radius         = Formulas.getAuthalicRadius(ellipsoid);
+        this.ellipsoid = ellipsoid;
+        authalicRadius = Formulas.getAuthalicRadius(ellipsoid);
         userToGeodetic = new PositionTransformer(crs, geographic, null);
     }
 
@@ -213,17 +224,33 @@ public class GeodeticCalculator {
      * @return a new geodetic calculator using the specified CRS.
      */
     public static GeodeticCalculator create(final CoordinateReferenceSystem crs) {
-        return new GeodeticCalculator(crs);
+        ArgumentChecks.ensureNonNull("crs", crs);
+        final Ellipsoid ellipsoid = ReferencingUtilities.getEllipsoid(crs);
+        if (ellipsoid == null) {
+            throw new IllegalArgumentException(Errors.format(Errors.Keys.IllegalCRSType_1, ReferencingUtilities.getInterface(crs)));
+        }
+        if (ellipsoid.isSphere()) {
+            return new GeodeticCalculator(crs, ellipsoid);
+        } else {
+            return new GeodesicsOnEllipsoid(crs, ellipsoid);
+        }
     }
 
     /**
      * Returns {@code true} if at least one of the properties identified by the given mask is invalid.
      */
-    private boolean isInvalid(final int mask) {
+    final boolean isInvalid(final int mask) {
         return (validity & mask) != mask;
     }
 
     /**
+     * Sets the properties specified by the given bitmask as valid.
+     */
+    final void setValid(final int mask) {
+        validity |= mask;
+    }
+
+    /**
      * Returns the Coordinate Reference System (CRS) in which {@code Position}s are represented, unless otherwise specified.
      * This is the CRS of all {@link Position} instances returned by methods in this class. This is also the default CRS
      * assumed by methods receiving a {@link Position} argument when the given position does not specify its CRS.
@@ -328,7 +355,7 @@ public class GeodeticCalculator {
      * and the current azimuth and distance.
      *
      * @return the destination (end point) represented in the CRS specified at construction time.
-     * @throws TransformException if the coordinates can not be transformed to {@linkplain #getPositionCRS() position CRS}.
+     * @throws TransformException if the coordinates can not be computed.
      * @throws IllegalStateException if the destination point, azimuth or distance have not been set.
      *
      * @see #getStartPoint()
@@ -372,8 +399,9 @@ public class GeodeticCalculator {
         ArgumentChecks.ensureFinite("longitude", longitude);
         φ2 = toRadians(max(Latitude.MIN_VALUE, min(Latitude.MAX_VALUE, latitude)));
         λ2 = toRadians(longitude);
-        validity |= END_POINT;
-        validity &= ~(STARTING_AZIMUTH | ENDING_AZIMUTH | GEODESIC_DISTANCE | RHUMBLINE_LENGTH);
+        setValid(END_POINT);
+        validity &= ~(STARTING_AZIMUTH | ENDING_AZIMUTH | GEODESIC_DISTANCE | RHUMBLINE_LENGTH | COEFFICIENTS_FOR_START_POINT);
+        // Coefficients for starting point are invalidated because the starting azimuth changed.
     }
 
     /**
@@ -385,8 +413,9 @@ public class GeodeticCalculator {
      *
      * @return the azimuth in degrees from -180° to +180°. 0° is toward North and values are increasing clockwise.
      * @throws IllegalStateException if the end point, azimuth or distance have not been set.
+     * @throws GeodesicException if the azimuth can not be computed.
      */
-    public double getStartingAzimuth() {
+    public double getStartingAzimuth() throws GeodesicException {
         if (isInvalid(STARTING_AZIMUTH)) {
             computeDistance();
         }
@@ -409,8 +438,8 @@ public class GeodeticCalculator {
         azimuth = toRadians(azimuth);
         dφ1 = cos(azimuth);                                 // sin(π/2 − θ)  =  cos(θ)
         dλ1 = sin(azimuth);                                 // cos(π/2 − θ)  =  sin(θ)
-        validity |= STARTING_AZIMUTH;
-        validity &= ~(END_POINT | ENDING_AZIMUTH | RHUMBLINE_LENGTH);
+        setValid(STARTING_AZIMUTH);
+        validity &= ~(END_POINT | ENDING_AZIMUTH | RHUMBLINE_LENGTH | COEFFICIENTS_FOR_START_POINT);
     }
 
     /**
@@ -421,13 +450,14 @@ public class GeodeticCalculator {
      *
      * @return the azimuth in degrees from -180° to +180°. 0° is toward North and values are increasing clockwise.
      * @throws IllegalStateException if the destination point, azimuth or distance have not been set.
+     * @throws GeodesicException if the azimuth can not be computed.
      */
-    public double getEndingAzimuth() {
+    public double getEndingAzimuth() throws GeodesicException {
         if (isInvalid(ENDING_AZIMUTH)) {
             if (isInvalid(END_POINT)) {
                 computeEndPoint();                      // Compute also ending azimuth from start point and distance.
             } else {
-                computeDistance();                         // Compute also ending azimuth from start point and end point.
+                computeDistance();                      // Compute also ending azimuth from start point and end point.
             }
         }
         return toDegrees(atan2(dλ2, dφ2));
@@ -441,10 +471,11 @@ public class GeodeticCalculator {
      *
      * @return the shortest distance in the unit of measurement given by {@link #getDistanceUnit()}.
      * @throws IllegalStateException if a point has not been set.
+     * @throws GeodesicException if the distance can not be computed.
      *
      * @see #getDistanceUnit()
      */
-    public double getGeodesicDistance() {
+    public double getGeodesicDistance() throws GeodesicException {
         if (isInvalid(GEODESIC_DISTANCE)) {
             computeDistance();
         }
@@ -463,7 +494,7 @@ public class GeodeticCalculator {
     public void setGeodesicDistance(final double distance) {
         ArgumentChecks.ensurePositive("distance", distance);
         geodesicDistance = distance;
-        validity |= GEODESIC_DISTANCE;
+        setValid(GEODESIC_DISTANCE);
         validity &= ~(END_POINT | ENDING_AZIMUTH | RHUMBLINE_LENGTH);
     }
 
@@ -509,17 +540,27 @@ public class GeodeticCalculator {
     }
 
     /**
-     * Computes the length of rhumb line from start point to end point.
-     *
-     * @see <a href="https://en.wikipedia.org/wiki/Rhumb_line">Rhumb line on Wikipedia</a>
+     * Ensures that the start point and end point are set.
+     * This method should be invoked at the beginning of {@link #computeDistance()}.
      *
-     * @todo if {@literal Δλ > 180}, must split in two segments.
+     * @throws IllegalStateException if the start point or end point has not been set.
      */
-    private void computeRhumbLine() {
+    final void canComputeDistance() {
         if (isInvalid(START_POINT | END_POINT)) {
             throw new IllegalStateException(Resources.format(
                     Resources.Keys.StartOrEndPointNotSet_1, Integer.signum(validity & START_POINT)));
         }
+    }
+
+    /**
+     * Computes the length of rhumb line from start point to end point.
+     *
+     * @see <a href="https://en.wikipedia.org/wiki/Rhumb_line">Rhumb line on Wikipedia</a>
+     *
+     * @todo if {@literal Δλ > 180}, must split in two segments.
+     */
+    void computeRhumbLine() {
+        canComputeDistance();
         final double Δλ = λ2 - λ1;
         final double Δφ = φ2 - φ1;
         final double factor;
@@ -539,12 +580,12 @@ public class GeodeticCalculator {
             final double α = atan(Δλ / ΔG);
             factor = Δφ / cos(α);
         }
-        rhumblineLength = radius * abs(factor);
-        validity |= RHUMBLINE_LENGTH;
+        rhumblineLength = authalicRadius * abs(factor);
+        setValid(RHUMBLINE_LENGTH);
     }
 
     /**
-     * Computes the geodetic distance and azimuths from the start point and end point.
+     * 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.
@@ -556,12 +597,10 @@ public class GeodeticCalculator {
      * </ul>
      *
      * @throws IllegalStateException if the distance or azimuth has not been set.
+     * @throws GeodesicException if an azimuth or the distance can not be computed.
      */
-    private void computeDistance() {
-        if (isInvalid(START_POINT | END_POINT)) {
-            throw new IllegalStateException(Resources.format(
-                    Resources.Keys.StartOrEndPointNotSet_1, Integer.signum(validity & START_POINT)));
-        }
+    void computeDistance() throws GeodesicException {
+        canComputeDistance();
         final double Δλ    = λ2 - λ1;           // No need to reduce to −π … +π range.
         final double sinΔλ = sin(Δλ);
         final double cosΔλ = cos(Δλ);
@@ -581,9 +620,23 @@ public class GeodeticCalculator {
          * Δσ = atan2(…) computes the same value but with better accuracy.
          */
         double Δσ = sinφ1*sinφ2 + cosφ1*cosφ2*cosΔλ;        // Actually Δσ = acos(…).
-        Δσ = atan2(hypot(dλ1, dφ1), Δσ);
-        geodesicDistance = radius * Δσ;
-        validity |= (STARTING_AZIMUTH | ENDING_AZIMUTH | GEODESIC_DISTANCE);
+        Δσ = atan2(hypot(dλ1, dφ1), Δσ);                    // Δσ ∈ [0…π].
+        geodesicDistance = authalicRadius * Δσ;
+        setValid(STARTING_AZIMUTH | ENDING_AZIMUTH | GEODESIC_DISTANCE);
+    }
+
+    /**
+     * Ensures that the start point, starting azimuth and geodesic distance are set.
+     * This method should be invoked at the beginning of {@link #computeEndPoint()}.
+     *
+     * @throws IllegalStateException if the start point, azimuth or distance has not been set.
+     */
+    final void canComputeEndPoint() {
+        if (isInvalid(START_POINT | STARTING_AZIMUTH | GEODESIC_DISTANCE)) {
+            throw new IllegalStateException(isInvalid(START_POINT)
+                    ? Resources.format(Resources.Keys.StartOrEndPointNotSet_1, 0)
+                    : Resources.format(Resources.Keys.AzimuthAndDistanceNotSet));
+        }
     }
 
     /**
@@ -594,21 +647,18 @@ public class GeodeticCalculator {
      * <p>The default implementation computes {@link #φ2}, {@link #λ2} and ∂φ/∂λ derivatives using
      * spherical formulas. Subclasses should override if they can provide ellipsoidal formulas.</p>
      *
-     * @throws IllegalStateException if the azimuth and the distance have not been set.
+     * @throws IllegalStateException if the start point, azimuth or distance has not been set.
+     * @throws GeodesicException if the end point or azimuth can not be computed.
      */
-    void computeEndPoint() {
-        if (isInvalid(START_POINT | STARTING_AZIMUTH | GEODESIC_DISTANCE)) {
-            throw new IllegalStateException(isInvalid(START_POINT)
-                    ? Resources.format(Resources.Keys.StartOrEndPointNotSet_1, 0)
-                    : Resources.format(Resources.Keys.AzimuthAndDistanceNotSet));
-        }
+    void computeEndPoint() throws GeodesicException {
+        canComputeEndPoint();                   // May throw IllegalStateException.
         final double vm    = hypot(dφ1, dλ1);
-        final double Δσ    = geodesicDistance / radius;
+        final double Δσ    = geodesicDistance / authalicRadius;
         final double sinΔσ = sin(Δσ);
         final double cosΔσ = cos(Δσ);
         final double sinφ1 = sin(φ1);
         final double cosφ1 = cos(φ1);
-        final double sinα1 = dλ1 / vm;          // α₁ is the azimuth at starting point.
+        final double sinα1 = dλ1 / vm;          // α₁ is the azimuth at starting point as a geographic (not arithmetic) angle.
         final double cosα1 = dφ1 / vm;
         final double sinΔσ_cosα1 = sinΔσ * cosα1;
         final double Δλy = sinΔσ * sinα1;
@@ -620,7 +670,7 @@ public class GeodeticCalculator {
         λ2  = IEEEremainder(λ1 + Δλ, 2*PI);
         dφ2 = cosΔσ*cosα1 - sinφ1/cosφ1 * sinΔσ;
         dλ2 = sinα1;
-        validity |= END_POINT;
+        setValid(END_POINT | ENDING_AZIMUTH);
     }
 
     /**
@@ -630,14 +680,19 @@ public class GeodeticCalculator {
      * some of them will need to be specified again.
      *
      * @see #setStartGeographicPoint(double, double)
+     * @throws GeodesicException if the end point or azimuth can not be computed.
      */
-    public void moveToEndPoint() {
+    public void moveToEndPoint() throws GeodesicException {
         if (isInvalid(END_POINT)) {
             computeEndPoint();
         }
-        φ1  = φ2;  dφ1 = dφ2;
-        λ1  = λ2;  dλ1 = dλ2;
-        validity = START_POINT | STARTING_AZIMUTH;
+        φ1 = φ2;  dφ1 = dφ2;
+        λ1 = λ2;  dλ1 = dλ2;
+        /*
+         * Following assumes that ENDING_AZIMUTH >>> 1 == STARTING_AZIMUTH.
+         * This is verified by GeodeticCalculatorTest.testMoveToEndPoint().
+         */
+        validity = ((validity & ENDING_AZIMUTH) >>> 1) | START_POINT;
     }
 
     /**
@@ -756,7 +811,7 @@ public class GeodeticCalculator {
         PathBuilder(final double εx) {
             super(ReferencingUtilities.getDimension(userToGeodetic.defaultCRS));
             dφf = dφ2;  dλf = dλ2;  φf = φ2;  λf = λ2;
-            tolerance = toDegrees(εx / radius);
+            tolerance = toDegrees(εx / authalicRadius);
             distance  = geodesicDistance;
             length    = rhumblineLength;
             flags     = validity;
@@ -780,7 +835,7 @@ public class GeodeticCalculator {
                 λ2 = λf;  dλ2 = dλf;
             } else {
                 geodesicDistance = distance * t;
-                validity |= GEODESIC_DISTANCE;
+                setValid(GEODESIC_DISTANCE);
                 computeEndPoint();
             }
             evaluateAtEndPoint();
@@ -863,7 +918,8 @@ public class GeodeticCalculator {
             final double α1 = t * (2*PI);
             dλ1 = sin(α1);
             dφ1 = cos(α1);
-            validity |= STARTING_AZIMUTH;
+            setValid(STARTING_AZIMUTH);
+            validity &= ~COEFFICIENTS_FOR_START_POINT;
             computeEndPoint();
             evaluateAtEndPoint();
             if (depth <= 1) {
@@ -964,7 +1020,7 @@ public class GeodeticCalculator {
                 resources.appendLabel(Vocabulary.Keys.GeodesicDistance, buffer);
                 buffer.append(' ').append(nf.format(distance))
                       .append(' ').append(unit).append(lineSeparator);
-            } catch (IllegalStateException e) {
+            } catch (IllegalStateException | GeodesicException e) {
                 // Ignore.
             }
         } catch (IOException e) {
diff --git a/core/sis-referencing/src/main/java/org/apache/sis/referencing/datum/DefaultEllipsoid.java b/core/sis-referencing/src/main/java/org/apache/sis/referencing/datum/DefaultEllipsoid.java
index 31f07f4..5be0bd7 100644
--- a/core/sis-referencing/src/main/java/org/apache/sis/referencing/datum/DefaultEllipsoid.java
+++ b/core/sis-referencing/src/main/java/org/apache/sis/referencing/datum/DefaultEllipsoid.java
@@ -111,7 +111,7 @@ import static org.apache.sis.util.ArgumentChecks.ensureNonNull;
  *
  * @author  Martin Desruisseaux (IRD, Geomatys)
  * @author  Cédric Briançon (Geomatys)
- * @version 0.8
+ * @version 1.0
  *
  * @see org.apache.sis.referencing.CommonCRS#ellipsoid()
  * @see org.apache.sis.referencing.factory.GeodeticAuthorityFactory#createEllipsoid(String)
@@ -536,7 +536,12 @@ public class DefaultEllipsoid extends AbstractIdentifiedObject implements Ellips
      * @param  λ2  longitude of second point (in decimal degrees).
      * @param  φ2  latitude  of second point (in decimal degrees).
      * @return the orthodromic distance (in the units of this ellipsoid's axis).
+     *
+     * @deprecated Replaced by {@link org.apache.sis.referencing.GeodeticCalculator}.
+     *
+     * @see <a href="https://issues.apache.org/jira/browse/SIS-386">SIS-386</a>
      */
+    @Deprecated
     public double orthodromicDistance(double λ1, double φ1, double λ2, double φ2) {
         λ1 = toRadians(λ1);
         φ1 = toRadians(φ1);
diff --git a/core/sis-referencing/src/main/java/org/apache/sis/referencing/datum/Sphere.java b/core/sis-referencing/src/main/java/org/apache/sis/referencing/datum/Sphere.java
index b9171fa..436aa5c 100644
--- a/core/sis-referencing/src/main/java/org/apache/sis/referencing/datum/Sphere.java
+++ b/core/sis-referencing/src/main/java/org/apache/sis/referencing/datum/Sphere.java
@@ -110,6 +110,7 @@ final class Sphere extends DefaultEllipsoid {
      * @return the orthodromic distance (in the units of this ellipsoid's axis).
      */
     @Override
+    @Deprecated
     public double orthodromicDistance(double λ1, double φ1, double λ2, double φ2) {
         φ1 = toRadians(φ1);
         φ2 = toRadians(φ2);
diff --git a/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/GeodesicException.java b/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/GeodesicException.java
new file mode 100644
index 0000000..ca6ce6a
--- /dev/null
+++ b/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/GeodesicException.java
@@ -0,0 +1,56 @@
+/*
+ * 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.operation;
+
+import org.opengis.referencing.operation.TransformException;
+
+
+/**
+ * Thrown when an error occurred while computing geodesic between two points.
+ *
+ * <div class="note"><b>API note:</b>
+ * defined as a sub-type of {@link TransformException} because some kind of coordinate operations are involved in
+ * geodesics calculation (e.g. transformation of latitudes and longitudes to coordinates on an auxiliary sphere).
+ * The starting and ending points can also be given in a CRS that require additional coordinate operations.</div>
+ *
+ * @author  Martin Desruisseaux (Geomatys)
+ * @version 1.0
+ * @since   1.0
+ * @module
+ */
+public class GeodesicException extends TransformException {
+    /**
+     * For cross-version compatibility.
+     */
+    private static final long serialVersionUID = -7844524593329236175L;
+
+    /**
+     * Constructs a new exception with no message.
+     */
+    public GeodesicException() {
+        super();
+    }
+
+    /**
+     * Constructs a new exception with the specified detail message.
+     *
+     * @param message  the detail message, or {@code null} if none.
+     */
+    public GeodesicException(final String message) {
+        super(message);
+    }
+}
diff --git a/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/projection/Mollweide.java b/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/projection/Mollweide.java
index 5f96038..e48903b 100644
--- a/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/projection/Mollweide.java
+++ b/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/projection/Mollweide.java
@@ -114,10 +114,10 @@ public class Mollweide extends NormalizedProjection {
          */
         if (abs(sinφ) != 1) {
             final double πsinφ = PI * sinφ;
-            int nbIte = MAXIMUM_ITERATIONS;
+            int nbIter = MAXIMUM_ITERATIONS;
             double Δθ;
             do {
-                if (--nbIte < 0) {
+                if (--nbIter < 0) {
                     throw new ProjectionException(Resources.format(Resources.Keys.NoConvergence));
                 }
                 Δθ = (θp + sin(θp) - πsinφ) / (1 + cos(θp));
diff --git a/core/sis-referencing/src/test/java/org/apache/sis/referencing/GeodesicsOnEllipsoidTest.java b/core/sis-referencing/src/test/java/org/apache/sis/referencing/GeodesicsOnEllipsoidTest.java
new file mode 100644
index 0000000..3d17a62
--- /dev/null
+++ b/core/sis-referencing/src/test/java/org/apache/sis/referencing/GeodesicsOnEllipsoidTest.java
@@ -0,0 +1,514 @@
+/*
+ * 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 java.util.Map;
+import java.util.HashMap;
+import java.util.Arrays;
+import org.opengis.geometry.DirectPosition;
+import org.opengis.referencing.crs.GeographicCRS;
+import org.opengis.referencing.operation.TransformException;
+import org.apache.sis.referencing.operation.GeodesicException;
+import org.apache.sis.internal.referencing.ReferencingUtilities;
+import org.apache.sis.internal.referencing.Formulas;
+import org.apache.sis.referencing.crs.HardCodedCRS;
+import org.apache.sis.math.MathFunctions;
+import org.apache.sis.measure.Units;
+import org.apache.sis.test.DependsOn;
+import org.apache.sis.test.DependsOnMethod;
+import org.junit.Test;
+
+import static java.lang.StrictMath.*;
+import static org.junit.Assert.*;
+
+
+/**
+ * Tests {@link GeodesicsOnEllipsoid}. Some of the tests in this class use data published by
+ * <a href="https://link.springer.com/content/pdf/10.1007%2Fs00190-012-0578-z.pdf">Algorithms
+ * for geodesics from Charles F. F. Karney (SRI International)</a>.
+ *
+ * @author  Matthieu Bastianelli (Geomatys)
+ * @author  Martin Desruisseaux (Geomatys)
+ * @version 1.0
+ * @since   1.0
+ * @module
+ */
+@DependsOn(GeodeticCalculatorTest.class)
+public final strictfp class GeodesicsOnEllipsoidTest extends GeodeticCalculatorTest {
+    /**
+     * The instance to be tested.
+     */
+    private GeodesicsOnEllipsoid testedEarth;
+
+    /**
+     * Values of local variables in {@link GeodesicsOnEllipsoid} methods. If values for the same key are added
+     * many times, the new values are appended after the current ones. This happen for example during iterations.
+     *
+     * @see GeodesicsOnEllipsoid#STORE_LOCAL_VARIABLES
+     */
+    private Map<String,double[]> localVariables;
+
+    /**
+     * {@code true} if {@link GeodesicsOnEllipsoid#store(String, double)} shall verify consistency instead than
+     * storing local variables.
+     */
+    private boolean verifyConsistency;
+
+    /**
+     * Creates a new test case.
+     */
+    public GeodesicsOnEllipsoidTest() {
+    }
+
+    /**
+     * Returns the calculator to use for testing purpose.
+     *
+     * @param  normalized  whether to force (longitude, latitude) axis order.
+     */
+    @Override
+    GeodeticCalculator create(final boolean normalized) {
+        testedEarth = new Calculator(normalized ? HardCodedCRS.WGS84 : HardCodedCRS.WGS84_φλ);
+        return testedEarth;
+    }
+
+    /**
+     * Creates a calculator which will store locale variables in the {@link #localVariables} map.
+     * This is used for the test cases that compare intermediate calculations with values provided
+     * in tables of Karney (2013) <cite>Algorithms for geodesics</cite> publication.
+     */
+    private void createTracked() {
+        localVariables = new HashMap<>();
+        testedEarth = new Calculator(HardCodedCRS.WGS84) {
+            @Override double computedToGiven(final double α1) {
+                return (abs(TRUNCATED_α1 - toDegrees(α1)) < 1E-3) ? toRadians(TRUNCATED_α1) : α1;
+            }
+
+            @Override void store(final String name, final double value) {
+                super.store(name, value);
+                if (verifyConsistency) {
+                    assertValueEquals(name, 0, value, STRICT, false);
+                } else {
+                    double[] values = localVariables.putIfAbsent(name, new double[] {value});
+                    if (values != null) {
+                        final int i = values.length;
+                        values = Arrays.copyOf(values, i+1);
+                        values[i] = value;
+                        localVariables.put(name, values);
+                    }
+                }
+            }
+        };
+    }
+
+    /**
+     * The {@link GeodesicsOnEllipsoid} implementation to use for tests. This implementation compares the
+     * quartic root computed by {@link GeodesicsOnEllipsoid#μ(double, double)} against the values computed
+     * by {@link MathFunctions#polynomialRoots(double...)}.
+     */
+    private static class Calculator extends GeodesicsOnEllipsoid {
+        /** Values needed for computation of μ. */
+        private double x, y;
+
+        /** Creates a new calculator for the given coordinate reference system. */
+        Calculator(final GeographicCRS crs) {
+            super(crs, ReferencingUtilities.getEllipsoid(crs));
+        }
+
+        /** Invoked when {@link GeodesicsOnEllipsoid} computed an intermediate value. */
+        @Override void store(final String name, final double value) {
+            if (name.length() == 1) {
+                switch (name.charAt(0)) {
+                    case 'x': x = value; break;
+                    case 'y': y = value; break;
+                    case 'μ': {
+                        double μ = 0;
+                        final double y2 = y*y;
+                        for (final double r : MathFunctions.polynomialRoots(-y2, -2*y2, (1 - x*x)-y2, 2, 1)) {
+                            if (r > μ) μ = r;
+                        }
+                        if (μ > 0) {
+                            /*
+                             * This assertion fails if y² is too close to zero. We that as a criterion for choosing
+                             * the threshold value that determine when to use Karney (2013) equation 57 instead.
+                             * That threshold is in `GeodesicOnEllipsoid.computeDistance()` method.
+                             */
+                            assertEquals("μ(x²,y²)", μ, value, abs(μ) * 1E-8);
+                        }
+                    }
+                }
+            }
+        }
+
+        /** Temporary implementation for allowing tests to pass. */
+        @Override @Deprecated final void computeRhumbLine() {
+            super.computeRhumbLine();
+            rhumblineLength *= ellipsoid.getSemiMajorAxis() / authalicRadius;
+            if (φ1 != 0 || φ2 != 0) rhumblineLength *= 1.1;     // DELETE ME!
+        }
+    }
+
+    /**
+     * Asserts that variable of the given name is equal to the given value. This is used for comparing an
+     * intermediate value computed by Apache SIS against the value published in a table of Karney (2013)
+     * <cite>Algorithms for geodesics</cite>.
+     *
+     * @param  name       name of the variable to verify.
+     * @param  index      if the same variable has been set many times, its index with numbering starting at 1. Otherwise 0.
+     * @param  expected   the expected value.
+     * @param  tolerance  tolerance threshold for comparison.
+     * @param  angular    whether the stored value needs to be converted from radians to degrees.
+     */
+    private void assertValueEquals(final String name, int index, final double expected, final double tolerance, final boolean angular) {
+        final double[] values = localVariables.get(name);
+        if (values != null) {
+            if (index > 0) {
+                index--;
+            } else if (values.length != 1) {
+                fail("Expected exactly one value for " + name + " but got " + Arrays.toString(values));
+            }
+            double value = values[index];
+            if (angular) value = toDegrees(value);
+            assertEquals(name, expected, value, tolerance);
+        } else if (GeodesicsOnEllipsoid.STORE_LOCAL_VARIABLES) {
+            fail("Missing value: " + name);
+        }
+    }
+
+    /**
+     * Takes a snapshot of some intermediate calculations in {@link GeodesicsOnEllipsoid}.
+     * If {@link GeodesicsOnEllipsoid#STORE_LOCAL_VARIABLES} is {@code true}, then we already
+     * have snapshots of those variables and this method will instead verify consistency.
+     */
+    private void snapshot() {
+        verifyConsistency = true;
+        testedEarth.snapshot();
+    }
+
+    /**
+     * Verifies that the constant parameters set by constructor are equal to the ones used by Karney
+     * for his example. Those values do not depend on the point being tested.
+     *
+     * <p><b>Source:</b>Karney (2013) table 1.</p>
+     */
+    private void verifyParametersForWGS84() {
+        assertEquals("1/f",  298.257223563,       testedEarth.ellipsoid.getInverseFlattening(), 1E-9);
+        assertEquals("a",    6378137,             testedEarth.ellipsoid.getSemiMajorAxis(),     STRICT);
+        assertEquals("b",    6356752.314245,      testedEarth.ellipsoid.getSemiMinorAxis(),     1E-6);
+        assertEquals("c",    6371007.180918,      testedEarth.authalicRadius,                   1E-6);
+        assertEquals("ℯ²",   0.00669437999014132, testedEarth.eccentricitySquared,              1E-17);
+        assertEquals("ℯ′²",  0.00673949674227643, testedEarth.secondEccentricitySquared,        1E-17);
+        assertEquals("n",    0.00167922038638370, testedEarth.thirdFlattening,                  1E-16);
+        assertEquals("unit", Units.METRE,         testedEarth.getDistanceUnit());
+    }
+
+    /**
+     * Tests solving the direct geodesic problem by a call to {@link GeodesicsOnEllipsoid#computeEndPoint()}.
+     * The data used by this test are provided by the example published in Karney (2013) table 2.
+     * Input values are:
+     *
+     * <ul>
+     *   <li>Starting longitude:  λ₁  =  any</li>
+     *   <li>Starting latitude:   φ₁  =  40°</li>
+     *   <li>Starting azimuth:    α₁  =  30°</li>
+     *   <li>Geodesic distance:   s₁₂ =  10 000 000 m</li>
+     * </ul>
+     *
+     * About 20 intermediate values are published in Karney table 2 and copied in this method body.
+     * Those values are verified if {@link GeodesicsOnEllipsoid#STORE_LOCAL_VARIABLES} is {@code true}.
+     *
+     * <p><b>Source:</b> Karney (2013), <u>Algorithms for geodesics</u> table 2.</p>
+     *
+     * @throws TransformException should never happen since this method does not apply map projections.
+     */
+    @Test
+    public void testComputeEndPoint() throws TransformException {
+        createTracked();
+        verifyParametersForWGS84();
+        testedEarth.setStartGeographicPoint(40, 10);
+        testedEarth.setStartingAzimuth(30);
+        testedEarth.setGeodesicDistance(10000000);
+        assertEquals("λ₁", toDegrees(testedEarth.λ1), 10.0, Formulas.ANGULAR_TOLERANCE);
+        assertEquals("φ₁", toDegrees(testedEarth.φ1), 40.0, Formulas.ANGULAR_TOLERANCE);
+        assertEquals("∆s", testedEarth.geodesicDistance, 10000000, STRICT);
+        /*
+         * The following method invocation causes calculation of all intermediate values.
+         * Some intermediate values are verified by `assertValueEquals` statements below.
+         * If some values are wrong, it is easier to check correctness with step-by-step
+         * debugging of `computeDistance()` method while keeping table 3 of Karney aside.
+         */
+        final DirectPosition endPoint = testedEarth.getEndPoint();
+        snapshot();
+        /*
+         * Start point on auxiliary sphere:
+         *
+         *   β₁  —  reduced latitude of starting point.
+         *   α₀  —  azimuth at equator in the direction given by α₁ azimuth at the β₁ reduced latitude.
+         *   σ₁  —  arc length from equatorial point E to starting point on auxiliary sphere.
+         *   ω₁  —  spherical longitude of the starting point (from equatorial point E) on the auxiliary sphere.
+         *   s₁  —  distance from equatorial point E to starting point along the geodesic.
+         */
+        assertValueEquals("β₁", 0, 39.90527714601, 1E-11, true);
+        assertValueEquals("α₀", 0, 22.55394020262, 1E-11, true);
+        assertValueEquals("σ₁", 0, 43.99915364500, 1E-11, true);
+        assertValueEquals("ω₁", 0, 20.32371827837, 1E-11, true);
+        /*
+         * End point on auxiliary sphere.
+         *
+         *   σ₂  —  arc lentgth from equator to the ending point σ₂ on auxiliary sphere.
+         *   β₂  —  reduced latitude of ending point.
+         *   ω₂  —  spherical longitude of the ending point (from equatorial point E) on the auxiliary sphere.
+         */
+        assertValueEquals("k²",     0, 0.00574802962857, 1E-14, false);
+        assertValueEquals("ε",      0, 0.00143289220416, 1E-14, false);
+        assertValueEquals("A₁",     0, 1.00143546236207, 1E-14, false);
+        assertValueEquals("I₁(σ₁)", 0, 0.76831538886412, 1E-14,  false);
+        assertValueEquals("s₁",     0,  4883990.626232,  1E-6,  false);
+        assertValueEquals("s₂",     0, 14883990.626232,  1E-6,  false);
+        assertValueEquals("τ₂",     0, 133.96266050208,  1E-11, true);
+        assertValueEquals("σ₂",     0, 133.92164083038,  1E-11, true);
+        assertValueEquals("α₂",     0, 149.09016931807,  1E-11, true);
+        assertValueEquals("β₂",     0,  41.69771809250,  1E-11, true);
+        assertValueEquals("ω₂",     0, 158.28412147112,  1E-11, true);
+        /*
+         * Conversion of end point coordinates to geodetic longitude.
+         *
+         *   λ₁  —  longitudes of the starting point, measured from equatorial point E.
+         *   λ₂  —  longitudes of the ending point, measured from Equatorial point E.
+         */
+        assertValueEquals("A₃",    0,   0.99928424306, 1E-11, false);
+        assertValueEquals("I₃(σ)", 1,   0.76773786069, 1E-11, false);
+        assertValueEquals("I₃(σ)", 2,   2.33534322170, 1E-11, false);
+        assertValueEquals("λ₁",    0,  20.26715038016, 1E-11, true);
+        assertValueEquals("λ₂",    0, 158.11205042393, 1E-11, true);
+        assertValueEquals("Δλ",    0, 137.84490004377, 1E-11, true);
+        /*
+         * Final result:
+         *
+         *   φ₂  —  end point latitude.
+         *   λ₂  —  end point longitude. Shifted by 10° compared to Karney because of λ₁ = 10°.
+         *   α₂  —  azimuth at end point.
+         */
+        assertEquals("φ₂",  41.79331020506, endPoint.getOrdinate(1), 1E-11);
+        assertEquals("λ₂", 147.84490004377, endPoint.getOrdinate(0), 1E-11);
+        assertEquals("α₂", 149.09016931807, testedEarth.getEndingAzimuth(), 1E-11);
+    }
+
+    /**
+     * Tests solving the inverse geodesic problem by a call to {@link GeodesicsOnEllipsoid#computeDistance()}
+     * for a short distance. The data used by this test are provided by the example published in Karney (2013)
+     * table 3. Input values are:
+     *
+     * <ul>
+     *   <li>Starting longitude:  λ₁ = any</li>
+     *   <li>Starting latitude:   φ₁ = -30.12345°</li>
+     *   <li>Ending   longitude:  λ₂ = λ₁ + 0.00005°</li>
+     *   <li>Ending   latitude:   φ₂ = -30.12344°</li>
+     * </ul>
+     *
+     * About 5 intermediate values are published in Karney table 3 and copied in this method body.
+     * Those values are verified if {@link GeodesicsOnEllipsoid#STORE_LOCAL_VARIABLES} is {@code true}.
+     *
+     * <p><b>Source:</b> Karney (2013), <u>Algorithms for geodesics</u> table 3.</p>
+     *
+     * @throws GeodesicException if a geodesic can not be computed.
+     */
+    @Test
+    @DependsOnMethod("testComputeEndPoint")
+    public void testComputeShortDistance() throws GeodesicException {
+        createTracked();
+        verifyParametersForWGS84();
+        testedEarth.setStartGeographicPoint(-30.12345, 2);
+        testedEarth.setEndGeographicPoint  (-30.12344, 2.00005);
+        /*
+         * The following method invocation causes calculation of all intermediate values.
+         * Some intermediate values are verified by `assertValueEquals` statements below.
+         * If some values are wrong, it is easier to check correctness with step-by-step
+         * debugging of `computeDistance()` method while keeping table 3 of Karney aside.
+         *
+         *   β₁  —  reduced latitude of starting point.
+         *   β₂  —  reduced latitude of ending point.
+         */
+        final double distance = testedEarth.getGeodesicDistance();
+        assertValueEquals("β₁", 0, -30.03999083821, 1E-11, true);
+        assertValueEquals("β₂", 0, -30.03998085491, 1E-11, true);
+        assertValueEquals("ωb", 0,   0.99748847744, 1E-11, false);
+        assertValueEquals("Δω", 0,   0.00005012589, 1E-11, true);
+        assertValueEquals("Δσ", 0,   0.00004452641, 1E-11, true);
+        assertValueEquals("α₁", 1,  77.04353354237, 1E-8,  true);
+        assertValueEquals("α₂", 1,  77.04350844913, 1E-8,  true);
+        /*
+         * Final result. Values are slightly different than the values published
+         * in Karney table 3 because GeodeticCalculator has done an iteration.
+         */
+        assertEquals("α₁", 77.04353354237, testedEarth.getStartingAzimuth(), 1E-8);     // Last 3 digits differ.
+        assertEquals("α₂", 77.04350844913, testedEarth.getEndingAzimuth(),   1E-8);
+        assertEquals("Δs", 4.944208,       distance,                         1E-6);
+    }
+
+    /**
+     * Result of Karney table 4, used as input in Karney table 5. We need to truncated that intermediate result
+     * to the same number of digits than Karney in order to get the numbers published in table 5 and 6.
+     * This value is used in {@link #testComputeNearlyAntipodal()} only.
+     */
+    private static final double TRUNCATED_α1 = 161.914;
+
+    /**
+     * Tests solving the inverse geodesic problem by a call to {@link GeodesicsOnEllipsoid#computeDistance()} for nearly
+     * antipodal points. The data used by this test are provided by the example published in Karney (2013) tables 4 to 6.
+     * Input values are:
+     *
+     * <ul>
+     *   <li>Starting longitude:  λ₁  =  any</li>
+     *   <li>Starting latitude:   φ₁  =  -30°</li>
+     *   <li>Ending   longitude:  λ₂  =  λ₁ + 179.8°</li>
+     *   <li>Ending   latitude:   φ₂  =  29.9°</li>
+     * </ul>
+     *
+     * Some intermediate values are published in Karney tables 4 to 6 and copied in this method body.
+     * Those values are verified if {@link GeodesicsOnEllipsoid#STORE_LOCAL_VARIABLES} is {@code true}.
+     *
+     * <p><b>Source:</b> Karney (2013), <u>Algorithms for geodesics</u> tables 4, 5 and 6.</p>
+     *
+     * @throws GeodesicException if a geodesic can not be computed.
+     */
+    @Test
+    @DependsOnMethod("testComputeShortDistance")
+    public void testComputeNearlyAntipodal() throws GeodesicException {
+        createTracked();
+        verifyParametersForWGS84();
+        testedEarth.setStartGeographicPoint(-30, 0);
+        testedEarth.setEndGeographicPoint(29.9, 179.8);
+        /*
+         * The following method invocation causes calculation of all intermediate values.
+         * Values β₁ and β₂ are keep constant during all iterations.
+         * Other values are given in Karney table 4.
+         */
+        final double distance = testedEarth.getGeodesicDistance();
+        assertValueEquals("β₁", 0, -29.91674771324, 1E-11, true);
+        assertValueEquals("β₂", 0,  29.81691642189, 1E-11, true);
+        assertValueEquals("x",  0,  -0.382344,      1E-6, false);
+        assertValueEquals("y",  0,  -0.220189,      1E-6, false);
+        assertValueEquals("μ",  0,   0.231633,      1E-6, false);
+        assertValueEquals("α₁", 1,   TRUNCATED_α1,  1E-3,  true);               // Initial value before iteration.
+        /*
+         * Following values are updated during iterations. Note that in order to get the same values than the ones
+         * published in Karney table 5, we need to truncate the α₁ initial value to the same number of digits than
+         * Karney. This is done automatically if GeodesicsOnEllipsoid.STORE_LOCAL_VARIABLES is true.
+         */
+        assertValueEquals("α₀",    1,    15.60939746414,    1E-11, true);
+        assertValueEquals("σ₁",    1,  -148.81253566596,    1E-11, true);
+        assertValueEquals("ω₁",    1,  -170.74896696128,    1E-11, true);
+        assertValueEquals("α₂",    1,    18.06728796231,    1E-11, true);
+        assertValueEquals("σ₂",    1,    31.08244976895,    1E-11, true);
+        assertValueEquals("ω₂",    1,     9.21345761110,    1E-11, true);
+        assertValueEquals("k²",    1,     0.00625153791662, 1E-14, false);
+        assertValueEquals("ε",     1,     0.00155801826780, 1E-14, false);
+        assertValueEquals("λ₁",    1,  -170.61483552458,    1E-11, true);
+        assertValueEquals("λ₂",    1,     9.18542009839,    1E-11, true);
+        assertValueEquals("Δλ",    1,   179.80025562297,    1E-11, true);
+        assertValueEquals("δλ",    1,     0.00025562297,    1E-11, true);
+        assertValueEquals("J(σ₁)", 1,    -0.00948040927640, 1E-14, false);
+        assertValueEquals("J(σ₂)", 1,     0.00031349128630, 1E-14, false);
+        assertValueEquals("Δm",    1, 57288.000110,         1E-6,  false);
+        assertValueEquals("dλ/dα", 1,     0.01088931716115, 1E-14, false);
+        assertValueEquals("δσ₁",   1,    -0.02347465519,    1E-11, true);
+        assertValueEquals("α₁",    2,   161.89052534481,    1E-11, true);
+        /*
+         * After second iteration.
+         */
+        assertValueEquals("δλ",    2,     0.00000000663,    1E-11, true);
+        assertValueEquals("α₁",    3,   161.89052473633,    1E-11, true);
+        /*
+         * After third iteration.
+         */
+        assertValueEquals("α₀",    3,    15.62947966537,    1E-11, true);
+        assertValueEquals("σ₁",    3,  -148.80913691776,    1E-11, true);
+        assertValueEquals("ω₁",    3,  -170.73634378066,    1E-11, true);
+        assertValueEquals("α₂",    3,    18.09073724574,    1E-11, true);
+        assertValueEquals("σ₂",    3,    31.08583447040,    1E-11, true);
+        assertValueEquals("ω₂",    3,     9.22602862110,    1E-11, true);
+        assertValueEquals("s₁",    3,  -16539979.064227,    1E-6, false);
+        assertValueEquals("s₂",    3,    3449853.763383,    1E-6, false);
+        assertValueEquals("Δs",    3,   19989832.827610,    1E-6, false);
+        assertValueEquals("λ₁",    3,  -170.60204712148,    1E-11, true);
+        assertValueEquals("λ₂",    3,     9.19795287852,    1E-11, true);
+        assertValueEquals("Δλ",    3,   179.80000000000,    1E-11, true);
+        /*
+         * Final result:
+         */
+        assertEquals("α₁", 161.89052473633, testedEarth.getStartingAzimuth(), 1E-11);
+        assertEquals("α₂",  18.09073724574, testedEarth.getEndingAzimuth(),   1E-11);
+        assertEquals("Δs", 19989832.827610, distance,                         1E-6);
+    }
+
+    /**
+     * Same test than the one defined in parent class, but with expected results modified for ellipsoidal formulas.
+     * Input points are from <a href="https://en.wikipedia.org/wiki/Great-circle_navigation#Example">Wikipedia</a>.
+     * Outputs were computed with GeographicLib.
+     *
+     * @throws TransformException if an error occurred while transforming coordinates.
+     */
+    @Test
+    @Override
+    @DependsOnMethod({"testComputeEndPoint", "testComputeNearlyAntipodal"})
+    public void testGeodesicDistanceAndAzimuths() throws TransformException {
+        final GeodeticCalculator c = create(false);
+        c.setStartGeographicPoint(-33.0, -71.6);            // Valparaíso
+        c.setEndGeographicPoint  ( 31.4, 121.8);            // Shanghai
+        /*
+         *                   α₁        α₂       distance
+         *   Spherical:     -94.41°   -78.42°   18743 km
+         *   Ellipsoidal:   -94.82°   -78.29°   18752 km
+         */
+        assertEquals(Units.METRE,                c.getDistanceUnit());
+        assertEquals("α₁",        -94.82071749,  c.getStartingAzimuth(),         1E-8);
+        assertEquals("α₂",        -78.28609385,  c.getEndingAzimuth(),           1E-8);
+        assertEquals("distance",   18752.493521, c.getGeodesicDistance() / 1000, 1E-6);
+        assertPositionEquals(31.4, 121.8,        c.getEndPoint(),                1E-12);    // Should be the specified value.
+        /*
+         * Keep start point unchanged, but set above azimuth and distance.
+         * Verify that we get the Shanghai coordinates.
+         */
+        c.setStartingAzimuth(-94.82);
+        c.setGeodesicDistance(18752494);
+        assertEquals("α₁",   -94.82,       c.getStartingAzimuth(), 1E-12);                  // Should be the specified value.
+        assertEquals("α₂",   -78.28678389, c.getEndingAzimuth(),   1E-8);
+        assertEquals("Δs", 18752.494,      c.getGeodesicDistance() / 1000, STRICT);         // Should be the specified value.
+        assertPositionEquals(31.4, 121.8,  c.getEndPoint(), 0.0002);
+    }
+
+    /**
+     * Tells whether failure to compute geodesic for the given data should cause the test case to fail.
+     * This is invoked by {@link #compareAgainstDataset()}.
+     *
+     * @param  expected  a row from the {@code $SIS_DATA/Tests/GeodTest.dat} file.
+     * @return whether the JUnit test should fail.
+     */
+    @Override
+    boolean isFailure(final double[] expected) {
+        final double φ1 = expected[COLUMN_φ1];
+        final double φ2 = expected[COLUMN_φ2];
+        final double Δλ = expected[COLUMN_λ2] - expected[COLUMN_λ1];
+        if (Δλ > 90 && max(abs(φ1), abs(φ2)) < 2E-4) {
+            return false;                           // Ignore equatorial case.
+        }
+        if (Δλ > 179 && abs(φ1 + φ2) < 0.002) {
+            return false;                           // Ignore antipodal case.
+        }
+        return super.isFailure(expected);
+    }
+}
diff --git a/core/sis-referencing/src/test/java/org/apache/sis/referencing/GeodeticCalculatorTest.java b/core/sis-referencing/src/test/java/org/apache/sis/referencing/GeodeticCalculatorTest.java
index b969d95..4d6e52f 100644
--- a/core/sis-referencing/src/test/java/org/apache/sis/referencing/GeodeticCalculatorTest.java
+++ b/core/sis-referencing/src/test/java/org/apache/sis/referencing/GeodeticCalculatorTest.java
@@ -28,6 +28,7 @@ import java.io.LineNumberReader;
 import org.opengis.geometry.DirectPosition;
 import org.opengis.referencing.cs.AxisDirection;
 import org.opengis.referencing.operation.TransformException;
+import org.apache.sis.referencing.operation.GeodesicException;
 import org.apache.sis.internal.referencing.j2d.ShapeUtilitiesExt;
 import org.apache.sis.internal.referencing.Formulas;
 import org.apache.sis.referencing.crs.HardCodedCRS;
@@ -58,11 +59,20 @@ import static org.opengis.test.Assert.*;
  *   <li>Charles Karney's <a href="https://geographiclib.sourceforge.io/">GeographicLib</a> implementation.</li>
  * </ul>
  *
+ * This base class tests calculator using spherical formulas.
+ * Subclass executes the same test but using ellipsoidal formulas.
+ *
  * @version 1.0
- * @since 1.0
+ * @since   1.0
  * @module
  */
-public final strictfp class GeodeticCalculatorTest extends TestCase {
+public strictfp class GeodeticCalculatorTest extends TestCase {
+    /**
+     * Creates a new test case.
+     */
+    public GeodeticCalculatorTest() {
+    }
+
     /**
      * Verifies that the given point is equals to the given latitude and longitude.
      *
@@ -71,7 +81,7 @@ public final strictfp class GeodeticCalculatorTest extends TestCase {
      * @param p  the actual position to verify.
      * @param ε  the tolerance threshold.
      */
-    private static void assertPositionEquals(final double φ, final double λ, final DirectPosition p, final double ε) {
+    static void assertPositionEquals(final double φ, final double λ, final DirectPosition p, final double ε) {
         assertEquals("φ", φ, p.getOrdinate(0), ε);
         assertEquals("λ", λ, p.getOrdinate(1), ε);
     }
@@ -90,21 +100,33 @@ public final strictfp class GeodeticCalculatorTest extends TestCase {
     }
 
     /**
-     * Returns the calculator to use for testing purpose. This classes uses a calculator for a sphere.
+     * Returns the calculator to use for testing purpose. Default implementation uses a calculator
+     * for a sphere. Subclasses should override for creating instances of the class to be tested.
      *
      * @param  normalized  whether to force (longitude, latitude) axis order.
      */
-    private static GeodeticCalculator create(final boolean normalized) {
-        return new GeodeticCalculator(normalized ? HardCodedCRS.SPHERE : HardCodedCRS.SPHERE_φλ);
+    GeodeticCalculator create(final boolean normalized) {
+        final GeodeticCalculator c = GeodeticCalculator.create(normalized ? HardCodedCRS.SPHERE : HardCodedCRS.SPHERE_φλ);
+        assertEquals(GeodeticCalculator.class, c.getClass());       // Expect the implementation with spherical formulas.
+        return c;
+    }
+
+    /**
+     * Returns a reference implementation for the given geodetic calculator.
+     */
+    private static Geodesic createReferenceImplementation(final GeodeticCalculator c) {
+        return new Geodesic(c.ellipsoid.getSemiMajorAxis(), 1/c.ellipsoid.getInverseFlattening());
     }
 
     /**
      * Tests some simple azimuth directions. The expected directions are approximately North, East,
      * South and West, but not exactly because of Earth curvature. The test verify merely that the
      * azimuths are approximately correct.
+     *
+     * @throws GeodesicException if a geodesic can not be computed.
      */
     @Test
-    public void testCardinalAzimuths() {
+    public void testCardinalAzimuths() throws GeodesicException {
         final GeodeticCalculator c = create(false);
         final double tolerance = 0.2;
         c.setStartGeographicPoint(20, 12);
@@ -116,9 +138,12 @@ public final strictfp class GeodeticCalculatorTest extends TestCase {
 
     /**
      * Tests azimuths at poles.
+     *
+     * @throws GeodesicException if a geodesic can not be computed.
      */
     @Test
-    public void testAzimuthAtPoles() {
+    @DependsOnMethod("testCardinalAzimuths")
+    public void testAzimuthAtPoles() throws GeodesicException {
         final GeodeticCalculator c = create(false);
         final double tolerance = 0.2;
         c.setStartGeographicPoint( 90,  30);
@@ -137,15 +162,20 @@ public final strictfp class GeodeticCalculatorTest extends TestCase {
 
     /**
      * Tests geodesic distances and rhumb line length on the equator.
+     *
+     * @throws GeodesicException if a geodesic can not be computed.
      */
     @Test
-    public void testDistanceAtEquator() {
-        final Random random = TestUtilities.createRandomNumberGenerator();
+    public void testDistanceAtEquator() throws GeodesicException {
+        final Random random = TestUtilities.createRandomNumberGenerator(86909845084528L);
         final GeodeticCalculator c = create(false);
-        final double r = c.ellipsoid.getSemiMajorAxis() * (PI / 180);
+        final double a = c.ellipsoid.getSemiMajorAxis();
+        final double b = c.ellipsoid.getSemiMinorAxis();
+        final double r = a * (PI / 180);
+        final double λmax = nextDown((b/a) * 180);
         c.setStartGeographicPoint(0, 0);
         for (int i=0; i<100; i++) {
-            final double x = 360 * random.nextDouble() - 180;
+            final double x = (2*λmax) * random.nextDouble() - λmax;
             c.setEndGeographicPoint(0, x);
             final double expected = abs(x) * r;
             assertEquals("Geodesic",   expected, c.getGeodesicDistance(), Formulas.LINEAR_TOLERANCE);
@@ -162,6 +192,7 @@ public final strictfp class GeodeticCalculatorTest extends TestCase {
      * @see <a href="https://en.wikipedia.org/wiki/Great-circle_navigation#Example">Great-circle navigation on Wikipedia</a>
      */
     @Test
+    @DependsOnMethod({"testCardinalAzimuths", "testDistanceAtEquator"})
     public void testGeodesicDistanceAndAzimuths() throws TransformException {
         final GeodeticCalculator c = create(false);
         c.setStartGeographicPoint(-33.0, -71.6);            // Valparaíso
@@ -208,10 +239,30 @@ public final strictfp class GeodeticCalculatorTest extends TestCase {
         final double λ = -71.6;
         c.setStartPoint(new DirectPosition2D(λ, φ));
         assertPositionEquals(λ, φ, c.getStartPoint(), Formulas.ANGULAR_TOLERANCE);
-
+        /*
+         * Test the "Valparaíso to Shanghai" geodesic given by Wikipedia, but with a larger tolerance threshold for allowing
+         * the test to pass with both spherical and ellipsoidal formula. It is not the purpose of this test to check accuracy;
+         * that check is done by `testGeodesicDistanceAndAzimuths()`.
+         */
         c.setStartingAzimuth(-94.41);
         c.setGeodesicDistance(18743000);
-        assertPositionEquals(121.8, 31.4, c.getEndPoint(), 0.01);
+        assertPositionEquals(121.8, 31.4, c.getEndPoint(), 0.2);
+    }
+
+    /**
+     * Tests {@link GeodeticCalculator#moveToEndPoint()}.
+     *
+     * @throws TransformException if an error occurred while transforming coordinates.
+     */
+    @Test
+    public void testMoveToEndPoint() throws TransformException {
+        // Following relationship is required by GeodeticCalculator.moveToEndPoint() implementation.
+        assertEquals(GeodeticCalculator.ENDING_AZIMUTH >>> 1, GeodeticCalculator.STARTING_AZIMUTH);
+        final GeodeticCalculator c = create(false);
+        c.setStartGeographicPoint(-33.0, -71.6);            // Valparaíso
+        c.setEndGeographicPoint  ( 31.4, 121.8);            // Shanghai
+        c.moveToEndPoint();
+        assertPositionEquals(31.4, 121.8, c.getStartPoint(), 1E-12);
     }
 
     /**
@@ -248,17 +299,18 @@ public final strictfp class GeodeticCalculatorTest extends TestCase {
     @DependsOnMethod("testUsingTransform")
     public void testGeodesicPath2D() throws TransformException {
         final GeodeticCalculator c = create(true);
-        final double tolerance = 0.05;
         c.setStartGeographicPoint(-33.0, -71.6);                                        // Valparaíso
         c.setEndGeographicPoint  ( 31.4, 121.8);                                        // Shanghai
         final Shape singleCurve = c.createGeodesicPath2D(Double.POSITIVE_INFINITY);
         final Shape multiCurves = c.createGeodesicPath2D(10000);                        // 10 km tolerance.
         /*
          * The approximation done by a single curve is not very good, but is easier to test.
+         * The coordinate at t=0.5 has a larger tolerance threshold because it varies slightly
+         * depending on whether we are testing with spherical or ellipsoidal formulas.
          */
-        assertPointEquals( -71.6, -33.0, ShapeUtilitiesExt.pointOnBezier(singleCurve, 0),   tolerance);
-        assertPointEquals(-238.2,  31.4, ShapeUtilitiesExt.pointOnBezier(singleCurve, 1),   tolerance);       // λ₂ = 121.8° - 360°
-        assertPointEquals(-159.2,  -6.8, ShapeUtilitiesExt.pointOnBezier(singleCurve, 0.5), tolerance);
+        assertPointEquals( -71.6, -33.0, ShapeUtilitiesExt.pointOnBezier(singleCurve, 0),   0.05);
+        assertPointEquals(-238.2,  31.4, ShapeUtilitiesExt.pointOnBezier(singleCurve, 1),   0.05);      // λ₂ = 121.8° - 360°
+        assertPointEquals(-159.2,  -6.9, ShapeUtilitiesExt.pointOnBezier(singleCurve, 0.5), 0.25);
         /*
          * The more accurate curve can not be simplified to a Java2D primitive.
          */
@@ -300,8 +352,8 @@ public final strictfp class GeodeticCalculatorTest extends TestCase {
     public void testBetweenRandomPoints() throws TransformException {
         final Random random = TestUtilities.createRandomNumberGenerator();
         final GeodeticCalculator c = create(false);
-        final Geodesic reference = new Geodesic(c.ellipsoid.getSemiMajorAxis(), 1/c.ellipsoid.getInverseFlattening());
-        final StatisticsFormat sf = VERBOSE ? StatisticsFormat.getInstance() : null;
+        final Geodesic reference = createReferenceImplementation(c);
+        Statistics[] errors = null;
         for (int i=0; i<100; i++) {
             final double φ1 = random.nextDouble() * 180 -  90;
             final double λ1 = random.nextDouble() * 360 - 180;
@@ -317,32 +369,45 @@ public final strictfp class GeodeticCalculatorTest extends TestCase {
             assertEquals("Starting azimuth",  expected.azi1, c.getStartingAzimuth(), Formulas.ANGULAR_TOLERANCE);
             assertEquals("Ending azimuth",    expected.azi2, c.getEndingAzimuth(),   Formulas.ANGULAR_TOLERANCE);
             assertTrue  ("Rhumb ≧ geodesic",  rhumbLine >= geodesic);
-            if (sf != null) {
+            if (VERBOSE) {
                 // Checks the geodesic path on only 10% of test data, because this computation is expensive.
                 if ((i % 10) == 0) {
-                    out.println(c);
-                    out.println(sf.format(geodesicPathFitness(c, 1000)));
+                    final Statistics[] stats = geodesicPathFitness(c, 1000);
+                    if (errors == null) {
+                        errors = stats;
+                    } else for (int j=0; j<errors.length; j++) {
+                        errors[j].combine(stats[j]);
+                    }
                 }
             }
         }
+        if (errors != null) {
+            out.println("Distance between points on Bézier curve and points on geodesic.");
+            out.println("∆x/r and ∆y/r should be smaller than 1:");
+            out.println(StatisticsFormat.getInstance().format(errors));
+        }
     }
 
     /**
      * Estimates the differences between the points on the Bézier curves and the points computed by geodetic calculator.
-     * This method estimates the length of the path returned by {@link GeodeticCalculator#createGeodesicPath2D(double)}
-     * and compares with the expected distance and azimuth at each point, by iterating over line segments and computing
-     * the geodesic distance of each segment. The state of the given calculator is modified by this method.
+     * This method approximates the Bézier curve by line segments. Then for each point of the approximated Bézier curve,
+     * this method computes the location of a close point on the geodesic (more specifically a point at the same geodesic
+     * distance from the start point). The distance in metres between the two points is measured and accumulated as a
+     * fraction of the expected resolution <var>r</var>. Consequently the values in ∆x/r and ∆y/r columns should be less
+     * than 1.
+     *
+     * <div class="note"><b>Note:</b> the state of the given calculator is modified by this method.</div>
      *
      * @param  resolution  tolerance threshold for the curve approximation, in metres.
      * @return statistics about errors relative to the resolution.
      */
     private static Statistics[] geodesicPathFitness(final GeodeticCalculator c, final double resolution) throws TransformException {
         final PathIterator iterator = c.createGeodesicPath2D(resolution).getPathIterator(null, Formulas.ANGULAR_TOLERANCE);
-        final Statistics   xError   = new Statistics("Δx/r");
-        final Statistics   yError   = new Statistics("Δy/r");
-        final Statistics   aErrors  = new Statistics("Δα (°)");
+        final Statistics   xError   = new Statistics("∆x/r");
+        final Statistics   yError   = new Statistics("∆y/r");
+        final Statistics   aErrors  = new Statistics("∆α (°)");
         final double       azimuth  = c.getStartingAzimuth();
-        final double       toMetres = (PI/180) * Formulas.getAuthalicRadius(c.ellipsoid);
+        final double       toMetres = (PI/180) * c.authalicRadius;
         final double[]     buffer   = new double[2];
         while (!iterator.isDone()) {
             switch (iterator.currentSegment(buffer)) {
@@ -367,6 +432,14 @@ public final strictfp class GeodeticCalculatorTest extends TestCase {
     }
 
     /**
+     * Column index for stating/ending latitude/longitude and geodesic distance.
+     * Used by {@link #compareAgainstDataset()} only.
+     */
+    static final int COLUMN_φ1 = 0, COLUMN_λ1 = 1, COLUMN_α1 = 2,
+                     COLUMN_φ2 = 3, COLUMN_λ2 = 4, COLUMN_α2 = 5,
+                     COLUMN_Δs = 6;
+
+    /**
      * Compares computations against values provided in <cite>Karney (2010) Test set for geodesics</cite>.
      * This is an optional test executed only if the {@code $SIS_DATA/Tests/GeodTest.dat} file is found.
      *
@@ -375,100 +448,117 @@ public final strictfp class GeodeticCalculatorTest extends TestCase {
      */
     @Test
     public void compareAgainstDataset() throws IOException, TransformException {
+        int noConvergenceCount = 0;
         try (LineNumberReader reader = OptionalTestData.GEODESIC.reader()) {
-            final GeodeticCalculator c = new GeodeticCalculator(HardCodedCRS.WGS84_φλ);
-            final Geodesic reference = new Geodesic(Formulas.getAuthalicRadius(c.ellipsoid), 0);
             final Random random = TestUtilities.createRandomNumberGenerator();
-            final double[] data = new double[7];
-            try {
-                String line;
-                while ((line = reader.readLine()) != null) {
-                    Arrays.fill(data, Double.NaN);
-                    final CharSequence[] split = CharSequences.split(line, ' ');
-                    for (int i=min(split.length, data.length); --i >= 0;) {
-                        data[i] = Double.parseDouble(split[i].toString());
+            final GeodeticCalculator c = create(false);
+            final boolean isSphere = c.ellipsoid.isSphere();
+            final double[] expected = new double[7];
+            String line;
+            while ((line = reader.readLine()) != null) {
+                Arrays.fill(expected, Double.NaN);
+                final CharSequence[] split = CharSequences.split(line, ' ');
+                for (int i=min(split.length, expected.length); --i >= 0;) {
+                    expected[i] = Double.parseDouble(split[i].toString());
+                }
+                /*
+                 * Choose randomly whether we will test the direct or inverse geodesic problem.
+                 * We execute only one test for each row instead than executing both tests,
+                 * for making sure that `GeodeticCalculator` never see the expected values.
+                 */
+                final boolean isTestingInverse = random.nextBoolean();
+                final double cosφ1 = abs(cos(toRadians(expected[COLUMN_φ1])));          // For adjusting longitude tolerance.
+                final double cosφ2 = abs(cos(toRadians(expected[COLUMN_φ2])));
+                double linearTolerance, latitudeTolerance, longitudeTolerance, azimuthTolerance;
+                if (isSphere) {
+                    /*
+                     * When spherical formulas are used instead than ellipsoidal formulas, an error up to 1% is expected
+                     * in distance calculations (source: Wikipedia). A yet larger error is observed for azimuth values,
+                     * especially near poles and between antipodal points. Following are empirical thresholds.
+                     */
+                    linearTolerance    = expected[COLUMN_Δs] * 0.01;
+                    latitudeTolerance  = toDegrees(linearTolerance / c.authalicRadius);
+                    longitudeTolerance = expected[COLUMN_φ2] > 89.5 ? 180 : latitudeTolerance / cosφ2;
+                    azimuthTolerance   = 0.5;                                   // About 8.8 metres at distance of 1 km.
+                    if (isTestingInverse) {
+                        final double Δλ = abs(expected[COLUMN_λ2] - expected[COLUMN_λ1]);
+                             if (Δλ > 179) azimuthTolerance = 100;
+                        else if (Δλ > 178) azimuthTolerance = 20;
+                        else if (Δλ > 175) azimuthTolerance = 10;
+                        else if (Δλ > 168) azimuthTolerance = 3;
+                        else if (Δλ > 158) azimuthTolerance = 1;
                     }
+                } else {
                     /*
-                     * We aim for an 1 cm accuracy. However when spherical formulas are used instead
-                     * than ellipsoidal formulas, an error up to 1% is expected (Wikipedia).
+                     * When ellipsoidal formulas are used, we aim for an 1 cm accuracy in coordinate values.
+                     * We also aim for azimuthd such as the error is less than 1 cm after the first 10 km.
+                     * If points are nearly antipodal, we relax the azimuth tolerance threshold to 1 meter.
                      */
-                    final double tolerance = data[6] * 0.01;                // 1% of distance.
-                    final double cosφ = abs(cos(toRadians(data[3])));       // For adjusting longitude tolerance.
-                    c.setStartGeographicPoint(data[0], data[1]);            // (φ₁, λ₁)
-                    if (random.nextBoolean()) {
-                        /*
-                         * Computes the end point from a distance and azimuth. The angular tolerance
-                         * is derived from the linear tolerance, except at pole where we disable the
-                         * check of longitude and azimuth values.
-                         */
-                        c.setStartingAzimuth (data[2]);
-                        c.setGeodesicDistance(data[6]);
-                        final double latitudeTolerance = tolerance * (1d/6371007 * (180/PI));
-                        final double longitudeTolerance;
-                        if (data[3] > 89.5) {
-                            longitudeTolerance = 180;               // TODO: remove after we use spherical formulas.
-                        } else {
-                            longitudeTolerance = latitudeTolerance / cosφ;
+                    linearTolerance    = Formulas.LINEAR_TOLERANCE;
+                    latitudeTolerance  = Formulas.ANGULAR_TOLERANCE;
+                    longitudeTolerance = Formulas.ANGULAR_TOLERANCE / cosφ2;
+                    azimuthTolerance   = Formulas.LINEAR_TOLERANCE * (180/PI) / 10000;
+                    if (isTestingInverse) {
+                        final double Δ = max(abs(180 - abs(expected[COLUMN_λ2] - expected[COLUMN_λ1])),
+                                                       abs(expected[COLUMN_φ1] + expected[COLUMN_φ2]));
+                        if (Δ < 1) {
+                            azimuthTolerance = 1 * (180/PI) / 10000;                // 1 meter for 10 km.
                         }
-                        final double azimuthTolerance = 0.5 / cosφ;
-                        compareGeodeticData(data, c, latitudeTolerance, longitudeTolerance,
-                                            azimuthTolerance, Formulas.LINEAR_TOLERANCE);
-                        /*
-                         * Replace the distance and azimuth values by values computed using spherical formulas,
-                         * then compare again with values computed by GeodeticCalculator but with low tolerance.
-                         */
-                        final GeodesicData gd = reference.Direct(data[0], data[1], data[2], data[6]);
-                        data[3] = gd.lat2;
-                        data[4] = gd.lon2;
-                        data[5] = gd.azi2;
-                    } else {
-                        /*
-                         * Compute the distance and azimuth values between two points. We perform
-                         * this test or the above test randomly instead of always executing both
-                         * of them for making sure that GeodeticCalculator never see the expected
-                         * values.
-                         */
-                        c.setEndGeographicPoint(data[3], data[4]);  // (φ₂, λ₂)
-                        compareGeodeticData(data, c,
-                                Formulas.ANGULAR_TOLERANCE,         // Latitude tolerance
-                                Formulas.ANGULAR_TOLERANCE,         // Longitude tolerance
-                                100 / cosφ, tolerance);             // Azimuth is inaccurate for reason not yet identified.
-                        /*
-                         * Replace the distance and azimuth values by values computed using spherical formulas,
-                         * then compare again with values computed by GeodeticCalculator but with low tolerance.
-                         */
-                        final GeodesicData gd = reference.Inverse(data[0], data[1], data[3], data[4]);
-                        data[2] = gd.azi1;
-                        data[5] = gd.azi2;
-                        data[6] = gd.s12;
                     }
-                    compareGeodeticData(data, c, Formulas.ANGULAR_TOLERANCE,
-                            Formulas.ANGULAR_TOLERANCE, 1E-4, Formulas.LINEAR_TOLERANCE);
                 }
-            } catch (AssertionError e) {
-                out.printf("Test failure at line %d%nGeodetic calculator is:%n", reader.getLineNumber());
-                out.println(c);
-                throw e;
+                /*
+                 * Set input values, compute then verify results. The azimuth tolerance is divided by cos(φ).
+                 * This is consistent with the ∂y/∂φ = 1/cos(φ) term in the Jacobian of Mercator projection.
+                 * An error of ε in the calculation of φ causes an error of ε/cos(φ) in calculation of tan(α).
+                 * Given that tan(α)+ε ≈ tan(α+ε) for small ε values, we assume that the error on α is also
+                 * proportional to 1/cos(φ).
+                 */
+                c.setStartGeographicPoint(expected[COLUMN_φ1], expected[COLUMN_λ1]);
+                if (isTestingInverse) {
+                    c.setEndGeographicPoint(expected[COLUMN_φ2], expected[COLUMN_λ2]);
+                } else {
+                    c.setStartingAzimuth (expected[COLUMN_α1]);
+                    c.setGeodesicDistance(expected[COLUMN_Δs]);
+                }
+                final DirectPosition start = c.getStartPoint();
+                final DirectPosition end   = c.getEndPoint();
+                try {
+                    assertEquals("φ₁", expected[COLUMN_φ1], start.getOrdinate(0),    Formulas.ANGULAR_TOLERANCE);
+                    assertEquals("λ₁", expected[COLUMN_λ1], start.getOrdinate(1),    Formulas.ANGULAR_TOLERANCE);
+                    assertEquals("α₁", expected[COLUMN_α1], c.getStartingAzimuth(),  azimuthTolerance / cosφ1);
+                    assertEquals("φ₂", expected[COLUMN_φ2], end.getOrdinate(0),      latitudeTolerance);
+                    assertEquals("λ₂", expected[COLUMN_λ2], end.getOrdinate(1),      longitudeTolerance);
+                    assertEquals("α₂", expected[COLUMN_α2], c.getEndingAzimuth(),    azimuthTolerance / cosφ2);
+                    assertEquals("∆s", expected[COLUMN_Δs], c.getGeodesicDistance(), linearTolerance);
+                } catch (GeodesicException | AssertionError e) {
+                    if (!isTestingInverse || e instanceof AssertionError || isFailure(expected)) {
+                        out.printf("Test failure at line %d: %s%n"
+                                + "The values provided in the test file are:%n"
+                                + "(φ₁,λ₁) = %16.12f %16.12f%n"
+                                + "(φ₂,λ₂) = %16.12f %16.12f%n"
+                                + "The values computed by the geodesic calculator are:%n",
+                                reader.getLineNumber(), e.getLocalizedMessage(),
+                                expected[0], expected[1], expected[3], expected[4]);
+                        out.println(c);
+                        throw e;
+                    }
+                    noConvergenceCount++;
+                }
             }
         }
+        assertTrue("Unexpected amount of no convergence errors.", noConvergenceCount <= 8);
     }
 
     /**
-     * Verifies that geodetic calculator results are equal to the given values.
-     * Order in the {@code data} array is as documented in {@link OptionalTestData#GEODESIC}.
+     * Tells whether failure to compute geodesic for the given data should cause the test case to fail.
+     * The default implementation always return {@code true}. Subclass can override if some points are
+     * known to fail.
+     *
+     * @param  expected  a row from the {@code $SIS_DATA/Tests/GeodTest.dat} file.
+     *         Use {@code COLUMN_*} constant for accessing values by column indices.
+     * @return whether the JUnit test should fail.
      */
-    private static void compareGeodeticData(final double[] expected, final GeodeticCalculator c,
-            final double latitudeTolerance, final double longitudeTolerance, final double azimuthTolerance,
-            final double linearTolerance) throws TransformException
-    {
-        final DirectPosition start = c.getStartPoint();
-        final DirectPosition end   = c.getEndPoint();
-        assertEquals("φ₁",  expected[0], start.getOrdinate(0),    Formulas.ANGULAR_TOLERANCE);
-        assertEquals("λ₁",  expected[1], start.getOrdinate(1),    Formulas.ANGULAR_TOLERANCE);
-        assertEquals("α₁",  expected[2], c.getStartingAzimuth(),  azimuthTolerance);
-        assertEquals("φ₂",  expected[3], end.getOrdinate(0),      latitudeTolerance);
-        assertEquals("λ₂",  expected[4], end.getOrdinate(1),      longitudeTolerance);
-        assertEquals("α₂",  expected[5], c.getEndingAzimuth(),    azimuthTolerance);
-        assertEquals("s₁₂", expected[6], c.getGeodesicDistance(), linearTolerance);
+    boolean isFailure(final double[] expected) {
+        return true;
     }
 }
diff --git a/core/sis-referencing/src/test/java/org/apache/sis/test/suite/ReferencingTestSuite.java b/core/sis-referencing/src/test/java/org/apache/sis/test/suite/ReferencingTestSuite.java
index 0c5c4e8..f46a2db 100644
--- a/core/sis-referencing/src/test/java/org/apache/sis/test/suite/ReferencingTestSuite.java
+++ b/core/sis-referencing/src/test/java/org/apache/sis/test/suite/ReferencingTestSuite.java
@@ -246,6 +246,7 @@ import org.junit.BeforeClass;
     org.apache.sis.referencing.cs.CodesTest.class,
     org.apache.sis.referencing.CRSTest.class,
     org.apache.sis.referencing.GeodeticCalculatorTest.class,
+    org.apache.sis.referencing.GeodesicsOnEllipsoidTest.class,
     org.apache.sis.internal.referencing.DefinitionVerifierTest.class,
     org.apache.sis.internal.referencing.CoordinateOperationsTest.class,
 
diff --git a/core/sis-utility/src/main/java/org/apache/sis/math/MathFunctions.java b/core/sis-utility/src/main/java/org/apache/sis/math/MathFunctions.java
index f0ee80c..48be82e 100644
--- a/core/sis-utility/src/main/java/org/apache/sis/math/MathFunctions.java
+++ b/core/sis-utility/src/main/java/org/apache/sis/math/MathFunctions.java
@@ -23,9 +23,13 @@ import org.apache.sis.util.ArgumentChecks;
 import org.apache.sis.util.resources.Errors;
 import org.apache.sis.internal.util.DoubleDouble;
 
+import static java.lang.Math.abs;
+import static java.lang.Math.sqrt;
 import static java.lang.Float.intBitsToFloat;
+import static java.lang.Float.floatToIntBits;
 import static java.lang.Float.floatToRawIntBits;
 import static java.lang.Double.longBitsToDouble;
+import static java.lang.Double.doubleToLongBits;
 import static java.lang.Double.doubleToRawLongBits;
 import static org.apache.sis.internal.util.Numerics.SIGN_BIT_MASK;
 import static org.apache.sis.internal.util.Numerics.SIGNIFICAND_SIZE;
@@ -223,7 +227,7 @@ public final class MathFunctions extends Static {
 
         // We have found a non-zero element. If it is the only one, returns its absolute value.
         double v2;
-        do if (i == 0) return Math.abs(v1);
+        do if (i == 0) return abs(v1);
         while ((v2 = vector[--i]) == 0);
 
         // If there is exactly 2 elements, use Math.hypot which is more robust than our algorithm.
@@ -406,7 +410,7 @@ public final class MathFunctions extends Static {
      * @since 0.6
      */
     public static double asinh(final double x) {
-        return Math.log(x + Math.sqrt(x*x + 1));
+        return Math.log(x + sqrt(x*x + 1));
     }
 
     /**
@@ -421,7 +425,7 @@ public final class MathFunctions extends Static {
      * @since 0.6
      */
     public static double acosh(final double x) {
-        return Math.log(x + Math.sqrt(x*x - 1));
+        return Math.log(x + sqrt(x*x - 1));
     }
 
     /**
@@ -601,7 +605,7 @@ public final class MathFunctions extends Static {
      * @return {@code true} if both values are equal given the tolerance threshold.
      */
     public static boolean epsilonEqual(final float v1, final float v2, final float ε) {
-        return (Math.abs(v1 - v2) <= ε) || Float.floatToIntBits(v1) == Float.floatToIntBits(v2);
+        return (abs(v1 - v2) <= ε) || floatToIntBits(v1) == floatToIntBits(v2);
     }
 
     /**
@@ -623,7 +627,7 @@ public final class MathFunctions extends Static {
      * @return {@code true} if both values are equal given the tolerance threshold.
      */
     public static boolean epsilonEqual(final double v1, final double v2, final double ε) {
-        return (Math.abs(v1 - v2) <= ε) || Double.doubleToLongBits(v1) == Double.doubleToLongBits(v2);
+        return (abs(v1 - v2) <= ε) || doubleToLongBits(v1) == doubleToLongBits(v2);
     }
 
     /**
@@ -777,7 +781,7 @@ public final class MathFunctions extends Static {
                     primes = Arrays.copyOf(primes, Math.min((index | 0xF) + 1, PRIMES_LENGTH_16_BITS));
                     do {
 testNextNumber:         while (true) {      // Simulate a "goto" statement (usually not recommanded...)
-                            final int stopAt = (int) Math.sqrt(n += 2);
+                            final int stopAt = (int) sqrt(n += 2);
                             int prime;
                             int j = 0;
                             do {
@@ -842,7 +846,7 @@ testNextNumber:         while (true) {      // Simulate a "goto" statement (usua
         if (number == 0) {
             return ArraysExt.EMPTY_INT;
         }
-        number = Math.abs(number);
+        number = abs(number);
         int[] divisors = new int[16];
         divisors[0] = 1;
         int count = 1;
@@ -852,7 +856,7 @@ testNextNumber:         while (true) {      // Simulate a "goto" statement (usua
          * values before that point, i.e. if n=p1*p2 and p2 is greater than 'sqrt', than p1
          * most be lower than 'sqrt'.
          */
-        final int sqrt = (int) Math.sqrt(number);               // Really want rounding toward 0.
+        final int sqrt = (int) sqrt(number);               // Really want rounding toward 0.
         for (int p,i=0; (p=primeNumberAt(i)) <= sqrt; i++) {
             if (number % p == 0) {
                 if (count == divisors.length) {
@@ -925,7 +929,7 @@ testNextNumber:         while (true) {      // Simulate a "goto" statement (usua
          */
         int minValue = Integer.MAX_VALUE;
         for (int i=0; i<numbers.length; i++) {
-            final int n = Math.abs(numbers[i]);
+            final int n = abs(numbers[i]);
             if (n <= minValue) {
                 minValue = n;
             }
@@ -937,7 +941,7 @@ testNextNumber:         while (true) {      // Simulate a "goto" statement (usua
          */
         int count = divisors.length;
         for (int i=0; i<numbers.length; i++) {
-            final int n = Math.abs(numbers[i]);
+            final int n = abs(numbers[i]);
             if (n != minValue) {
                 for (int j=count; --j>0;) {         // Do not test j==0, since divisors[0] ==  1.
                     if (n % divisors[j] != 0) {
@@ -948,4 +952,229 @@ testNextNumber:         while (true) {      // Simulate a "goto" statement (usua
         }
         return ArraysExt.resize(divisors, count);
     }
+
+    /**
+     * Returns the real (non-complex) roots of a polynomial equation having the given coefficients.
+     * This method returns the <var>x</var> values for which <var>y</var>=0 in the following equation:
+     *
+     * <blockquote><var>y</var> =
+     * <var>c<sub>0</sub></var> +
+     * <var>c<sub>1</sub></var>⋅<var>x</var> +
+     * <var>c<sub>2</sub></var>⋅<var>x</var><sup>2</sup> +
+     * <var>c<sub>3</sub></var>⋅<var>x</var><sup>3</sup> + … +
+     * <var>c<sub>n</sub></var>⋅<var>x</var><sup>n</sup>
+     * </blockquote>
+     *
+     * Current implementation can resolve polynomials described by a maximum of 5 coefficients, ignoring
+     * leading and trailing zeros. They correspond to linear, quadratic, cubic and quartic polynomials.
+     *
+     * @param  coefficients  the <var>c<sub>0</sub></var>, <var>c<sub>1</sub></var>, <var>c<sub>2</sub></var>, …
+     *                       <var>c<sub>n</sub></var> coefficients, in that order.
+     * @return the non-complex roots, or an empty array if none.
+     * @throws UnsupportedOperationException if given arguments contain more non-zero coefficients than this method can handle.
+     *
+     * @see java.awt.geom.QuadCurve2D#solveQuadratic(double[])
+     * @see java.awt.geom.CubicCurve2D#solveCubic(double[])
+     *
+     * @since 1.0
+     */
+    public static double[] polynomialRoots(final double... coefficients) {
+        int upper = coefficients.length;
+        while (upper > 0) {
+            double a = coefficients[--upper];                   // Coefficient of the term with highest exponent.
+            if (a == 0) {
+                continue;                                       // Search the last non-zero coefficient.
+            }
+            double c;                                           // Coefficient of the term with lowest exponent.
+            int lower = 0;
+            while ((c = coefficients[lower]) == 0) lower++;
+            switch (upper - lower) {
+                /*
+                 * c = 0
+                 *
+                 * Can not compute x. We could return an arbitrary value if c = 0, but we rather return
+                 * an empty array for keeping the number of root equals to the highest exponent.
+                 */
+                case 0: {
+                    break;
+                }
+                /*
+                 * ax + c = 0    →    x = -c/a
+                 */
+                case 1: {
+                    final double x = -c / a;
+                    if (Double.isNaN(x)) break;
+                    return new double[] {x};
+                }
+                /*
+                 * ax² + bx + c = 0    →    x = (-b ± √[b² - 4ac]) / (2a)
+                 *
+                 * Above equation is numerically unstable. More stable algorithm is given
+                 * by Numerical Recipes §5.6 (quadratic equation), which is applied below.
+                 */
+                case 2: {
+                    final double b = coefficients[lower + 1];
+                    final double q = -0.5 * (b + Math.copySign(sqrt(b*b - 4*a*c), b));
+                    final double x1 = q/a;
+                    final double x2 = c/q;
+                    if (Double.isNaN(x1) && Double.isNaN(x2)) break;
+                    return (x1 != x2) ? new double[] {x1, x2} : new double[] {x1};
+                }
+                /*
+                 * x³ + ax² + bx + c = 0
+                 *
+                 * Numerical Recipes §5.6 (cubic equation) is applied below.
+                 * Solution usually have either 1 or 3 roots.
+                 */
+                case 3: {
+                    return refineRoots(coefficients, solveCubic(
+                                       coefficients[lower + 2] / a,
+                                       coefficients[lower + 1] / a,
+                                       c / a, false));
+                }
+                /*
+                 * x⁴ + ax³ + bx² + cx + d = 0
+                 *
+                 * https://dlmf.nist.gov/1.11 in "Quartic equations" section.
+                 * This algorithm reduces the equation to a cubic equation.
+                 */
+                case 4: {
+                    double b,d;
+                    d = c / a;
+                    c = coefficients[lower + 1] / a;
+                    b = coefficients[lower + 2] / a;
+                    a = coefficients[lower + 3] / a;
+                    final double a2 = a*a;
+                    final double p = -3./8   * (a2)     +  b;
+                    final double q =  1./8   * (a2*a)   -  1./2  * (a*b)   +  c;
+                    final double r = -3./256 * (a2*a2)  +  1./16 * (a2*b)  -  1./4 * (a*c) + d;
+                    final double[] roots = solveCubic(-2*p, p*p-4*r, q*q, true);
+                    if (roots.length != 4) break;
+                    for (int i=0; i<3; i++) {
+                        roots[i] = sqrt(-roots[i]);
+                    }
+                    if (isPositive(q)) {
+                        for (int i=0; i<3; i++) {
+                            roots[i] = -roots[i];
+                        }
+                    }
+                    final double α = roots[0], β = roots[1], γ = roots[2];
+                    final double s = α + β + γ;
+                    if (Double.isNaN(s)) break;
+                    roots[0] = s/2 - (a /= 4);
+                    roots[1] = (+α - β - γ)/2 - a;
+                    roots[2] = (-α + β - γ)/2 - a;
+                    roots[3] = (-α - β + γ)/2 - a;
+                    return refineRoots(coefficients, removeDuplicated(roots));
+                }
+                default: {
+                    throw new UnsupportedOperationException();
+                }
+            }
+            break;
+        }
+        return ArraysExt.EMPTY_DOUBLE;
+    }
+
+    /**
+     * Solves cubic equation x³ + ax² + bx + c = 0. The solution before simplification has either 1 or 3 roots.
+     * The {@code quartic} argument specifies whether this cubic equation is used as a step for solving a quartic equation:
+     *
+     * <ul>
+     *   <li>If {@code true}, then we are interested only in the 3 roots solution and we do not check for duplicated values.
+     *       The length of returned array is 4 for allowing the caller to reuse the same array.</li>
+     *   <li>If {@code false}, then this method may simplify the 3 roots to 2 roots if two of them are equal,
+     *       or may return the 1 root solution. The length of returned array is the number of roots.</li>
+     * </ul>
+     */
+    private static double[] solveCubic(double a, double b, double c, final boolean quartic) {
+        final double Q = (a*a - 3*b) / 9;                           // Q from Numerical Recipes 5.6.10.
+        final double R = (a*(a*a - 4.5*b) + 13.5*c) / 27;           // R from Numerical Recipes 5.6.10.
+        final double Q3 = Q*Q*Q;
+        final double R2 = R*R;
+        a /= 3;                                                     // Last term of Numerical Recipes 5.6.12, 17 and 18.
+        if (R2 < Q3) {
+            /*
+             * Numerical Recipes 5.6.11 and 5.6.12 uses acos(R/sqrt(Q³)). It is possible to rewrite as
+             * atan2(sqrt(Q3 - R2), R) using the  cos(θ) = 1/√[1 + tan²θ]  trigonometric identity, but
+             * this substitution seems to decrease accuracy instead of increasing it in our tests.
+             */
+            b = Math.acos(R/sqrt(Q3)) / 3;                          // θ from Numerical recipes 5.6.11, then b = θ/3.
+            c = -2 * sqrt(Q);                                       // First part of Numerical Recipes 5.6.12.
+            double[] roots = new double[quartic ? 4 : 3];
+            roots[2] = c*Math.cos(b - 2*Math.PI/3) - a;             // TODO: try Math.fma with JDK9.
+            roots[1] = c*Math.cos(b + 2*Math.PI/3) - a;
+            roots[0] = c*Math.cos(b) - a;
+            if (!quartic) {
+                roots = removeDuplicated(roots);
+            }
+            return roots;
+        }
+        if (!quartic) {
+            b = -Math.copySign(Math.cbrt(abs(R) + sqrt(R2 - Q3)), R);       // A from Numerical Recipes 5.6.15.
+            final double x = (b == 0 ? 0 : b + Q/b) - a;
+            if (!Double.isNaN(x)) {
+                return new double[] {x};
+            }
+        }
+        return ArraysExt.EMPTY_DOUBLE;
+    }
+
+    /**
+     * Remove duplicated values in the given array. This method is okay only for very small arrays (3 or 4 elements).
+     * Duplicated values should be very rare and occur mostly as a consequence of rounding errors while computing the
+     * roots of polynomial equations. Because if the algebraic solution has less roots than what we would expect from
+     * the largest exponent (for example ax² + bx = 0 has only one root instead of two), then {@link #polynomialRoots}
+     * should have reduced the equation to a lower degrees (ax + b = 0 in above example), in which case there is no
+     * duplicated roots to remove.
+     */
+    private static double[] removeDuplicated(double[] roots) {
+        int i = 1;
+next:   while (i < roots.length) {
+            for (int j=i; --j >= 0;) {
+                if (roots[j] == roots[i]) {
+                    roots = ArraysExt.remove(roots, i, 1);
+                    continue next;
+                }
+            }
+            i++;
+        }
+        return roots;
+    }
+
+    /**
+     * Tries to improves accuracy of polynomial roots by applying small displacements
+     * to the <var>x</var> values using ∂y/∂x derivative around those values.
+     *
+     * <div class="note"><b>Purpose:</b>
+     * this refinement is significant in a {@link org.apache.sis.referencing.GeodesicsOnEllipsoid}
+     * test checking the value of an μ(x²,y²) function.</div>
+     *
+     * @param  coefficients  the user-specified coefficients.
+     * @param  roots         the roots. This array will be modified in place.
+     * @return {@code roots}.
+     */
+    private static double[] refineRoots(final double[] coefficients, final double[] roots) {
+        for (int i=0; i < roots.length; i++) {
+            double ymin = Double.POSITIVE_INFINITY;
+            double x = roots[i];
+            double dx;
+            do {
+                double px = 1;                              // Power of x: 1, x¹, x², x³, …
+                double dy = 0;                              // First derivative of polynomial at x.
+                double y  = coefficients[0];                // Value of polynomial at x.
+                double ey = 0, edy = 0;                     // Error terms for Kahan summation algorithm.
+                for (int j=1; j<coefficients.length; j++) {
+                    final double c = coefficients[j];
+                    double s;
+                    s = c * (px *  i) + edy; edy = s + (dy - (dy += s));        // Kahan summation of dy.
+                    s = c * (px *= x) +  ey;  ey = s + ( y - ( y += s));        // Kahan summation of y.
+                }
+                if (!(ymin > (ymin = abs(y)))) break;       // If result not better than previous result, stop.
+                roots[i] = x;
+                dx = y/dy;
+            } while (x != (x -= dx) && Double.isFinite(x));
+        }
+        return roots;
+    }
 }
diff --git a/core/sis-utility/src/test/java/org/apache/sis/math/MathFunctionsTest.java b/core/sis-utility/src/test/java/org/apache/sis/math/MathFunctionsTest.java
index 210e6f0..b899114 100644
--- a/core/sis-utility/src/test/java/org/apache/sis/math/MathFunctionsTest.java
+++ b/core/sis-utility/src/test/java/org/apache/sis/math/MathFunctionsTest.java
@@ -438,4 +438,47 @@ public final strictfp class MathFunctionsTest extends TestCase {
         assertArrayEquals(new int[] {1, 2}, commonDivisors(-4, 2));
         assertArrayEquals(new int[] {}, commonDivisors(0));
     }
+
+    /**
+     * Tests {@link MathFunctions#polynomialRoots(double...)} with values given in NIST example.
+     * Source: <a href="https://dlmf.nist.gov/1.11">NIST Digital Library of Mathematical Functions</a>.
+     */
+    @Test
+    public void testQuarticRoots() {
+        /*
+         * Given polynomial:  f(z) = z⁴ - 4z³ + 5z + 2
+         * Reduced:           g(w) = w⁴ - 6w² - 3w + 4
+         * Resolvent cubic:   z³ + 12z² + 20z + 9 = 0
+         */
+        assertArrayEquals(new double[] {
+             3.5615528128088303,            // ½(3 + √17)
+             1.618033988749895,             // ½(1 + √5)
+            -0.6180339887498949,            // ½(1 - √5)
+            -0.5615528128088303             // ½(3 - √17)
+        }, polynomialRoots(2, 5, 0, -4, 1), 1E-14);
+    }
+
+    /**
+     * Tests {@link MathFunctions#polynomialRoots(double...)} with random values.
+     */
+    @Test
+    public void testPolynomialRoots() {
+        final Random r = TestUtilities.createRandomNumberGenerator(62308330810132L);
+        final double[] coefficients = new double[5];
+        for (int n=2; n < coefficients.length; n++) {           // Number of coefficients to fill.
+            for (int t=n*10; --t >= 0;) {                       // Number of tests to execute.
+                for (int i=0; i<n; i++) {
+                    coefficients[i] = (r.nextInt(8) != 0) ? r.nextDouble() * 20 - 10 : 0;
+                }
+                for (final double x : polynomialRoots(coefficients)) {
+                    double p = 1;
+                    double y = coefficients[0];
+                    for (int i=1; i<coefficients.length; i++) {
+                        y += coefficients[i] * (p *= x);
+                    }
+                    assertEquals(0, y, 1E-12);
+                }
+            }
+        }
+    }
 }
diff --git a/core/sis-utility/src/test/java/org/apache/sis/test/TestCase.java b/core/sis-utility/src/test/java/org/apache/sis/test/TestCase.java
index da96ed4..b2ae8f1 100644
--- a/core/sis-utility/src/test/java/org/apache/sis/test/TestCase.java
+++ b/core/sis-utility/src/test/java/org/apache/sis/test/TestCase.java
@@ -126,6 +126,9 @@ public abstract strictfp class TestCase {
         out = new PrintWriter(buffer = new StringWriter());
         VERBOSE = Boolean.getBoolean(TestConfiguration.VERBOSE_OUTPUT_KEY);
         RUN_EXTENSIVE_TESTS = Boolean.getBoolean(TestConfiguration.EXTENSIVE_TESTS_KEY);
+        if (VERBOSE) {
+            System.setErr(System.out);      // For avoiding log records to be interleaved with block of text.
+        }
     }
 
     /**
diff --git a/storage/sis-storage/src/main/java/org/apache/sis/index/tree/QuadTree.java b/storage/sis-storage/src/main/java/org/apache/sis/index/tree/QuadTree.java
index 92fcdea..9bbce43 100644
--- a/storage/sis-storage/src/main/java/org/apache/sis/index/tree/QuadTree.java
+++ b/storage/sis-storage/src/main/java/org/apache/sis/index/tree/QuadTree.java
@@ -25,6 +25,7 @@ import org.apache.sis.geometry.Envelope2D;
 import org.apache.sis.distance.LatLonPointRadius;
 import org.apache.sis.referencing.CommonCRS;
 import org.apache.sis.referencing.GeodeticCalculator;
+import org.apache.sis.referencing.operation.GeodesicException;
 
 /**
  * Implementation of Quad Tree Index. Insertion algorithm implemented based on
@@ -342,9 +343,9 @@ public class QuadTree {
         if (node == null) {
             return matches;
         } else if (node.getNodeType() != NodeType.GRAY) {
-            if (node.getNodeType() == NodeType.WHITE)
+            if (node.getNodeType() == NodeType.WHITE) {
                 return matches;
-            else {
+            } else {
                 QuadTreeData[] data = node.getData();
                 for (int i = 0; i < node.getCount(); i++) {
                     DirectPosition2D latLon = data[i].getLatLon();
@@ -352,7 +353,11 @@ public class QuadTree {
                     synchronized (calculator) {
                         calculator.setStartGeographicPoint(latLon.y, latLon.x);
                         calculator.setEndGeographicPoint(point.y, point.x);
-                        distance = calculator.getGeodesicDistance();
+                        try {
+                            distance = calculator.getGeodesicDistance();
+                        } catch (GeodesicException e) {
+                            throw new RuntimeException(e);
+                        }
                     }
                     if (distance <= radiusKM) {
                         matches.add(data[i]);
@@ -360,7 +365,6 @@ public class QuadTree {
                 }
                 return matches;
             }
-
         } else {
             Rectangle2D swRectangle = new Rectangle2D.Double(nodeRegion.getX(), nodeRegion.getY(),
                     nodeRegion.getWidth() / 2, nodeRegion.getHeight() / 2);


Mime
View raw message