/*
 * 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.awt.Shape;
import java.awt.geom.Path2D;
import java.awt.geom.Point2D;
import java.awt.geom.Rectangle2D;
import java.awt.geom.PathIterator;
import java.util.Arrays;
import java.util.Random;
import java.io.IOException;
import java.io.LineNumberReader;
import org.opengis.geometry.DirectPosition;
import org.opengis.referencing.cs.AxisDirection;
import org.apache.sis.internal.referencing.j2d.ShapeUtilitiesExt;
import org.apache.sis.internal.referencing.Formulas;
import org.apache.sis.referencing.crs.HardCodedCRS;
import org.apache.sis.geometry.DirectPosition2D;
import org.apache.sis.util.CharSequences;
import org.apache.sis.math.StatisticsFormat;
import org.apache.sis.math.Statistics;
import org.apache.sis.measure.Units;
import org.apache.sis.test.widget.VisualCheck;
import org.apache.sis.test.OptionalTestData;
import org.apache.sis.test.DependsOnMethod;
import org.apache.sis.test.TestUtilities;
import org.apache.sis.test.TestCase;
import net.sf.geographiclib.Geodesic;
import net.sf.geographiclib.GeodesicData;
import org.junit.Test;

import static java.lang.StrictMath.*;
import static org.opengis.test.Assert.*;
import static org.apache.sis.internal.metadata.ReferencingServices.AUTHALIC_RADIUS;


/**
 * Tests {@link GeodeticCalculator}. Test values come from the following sources:
 *
 * <ul>
 *   <li><a href="https://en.wikipedia.org/wiki/Great-circle_navigation#Example">Great-circle navigation on Wikipedia.</a></li>
 *   <li><a href="http://doi.org/10.5281/zenodo.32156">Karney, C. F. F. (2010). Test set for geodesics [Data set]. Zenodo.</a></li>
 *   <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
 * @module
 */
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.
     *
     * @param φ  the expected latitude value, in degrees.
     * @param λ  the expected longitude value, in degrees.
     * @param p  the actual position to verify.
     * @param ε  the tolerance threshold.
     */
    static void assertPositionEquals(final double φ, final double λ, final DirectPosition p, final double ε) {
        assertEquals("φ", φ, p.getOrdinate(0), ε);
        assertEquals("λ", λ, p.getOrdinate(1), ε);
    }

    /**
     * Asserts that a Java2D point is equal to the expected value. Used for verifying geodesic paths.
     *
     * @param x  the expected <var>x</var> coordinates.
     * @param y  the expected <var>y</var> coordinates.
     * @param p  the actual position to verify.
     * @param ε  the tolerance threshold.
     */
    private static void assertPointEquals(final double x, final double y, final Point2D p, final double ε) {
        assertEquals("x", x, p.getX(), ε);
        assertEquals("y", y, p.getY(), ε);
    }

    /**
     * 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.
     */
    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.semiMajorAxis, 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.
     */
    @Test
    public void testCardinalAzimuths() {
        final GeodeticCalculator c = create(false);
        final double tolerance = 0.2;
        c.setStartGeographicPoint(20, 12);
        c.setEndGeographicPoint(20, 13);  assertEquals("East",   90, c.getStartingAzimuth(), tolerance);
        c.setEndGeographicPoint(21, 12);  assertEquals("North",   0, c.getStartingAzimuth(), tolerance);
        c.setEndGeographicPoint(20, 11);  assertEquals("West",  -90, c.getStartingAzimuth(), tolerance);
        c.setEndGeographicPoint(19, 12);  assertEquals("South", 180, c.getStartingAzimuth(), tolerance);
    }

    /**
     * Tests azimuths at poles.
     */
    @Test
    @DependsOnMethod("testCardinalAzimuths")
    public void testAzimuthAtPoles() {
        final GeodeticCalculator c = create(false);
        final double tolerance = 0.2;
        c.setStartGeographicPoint( 90,  30);
        c.setEndGeographicPoint  ( 20,  20);  assertEquals(-170, c.getStartingAzimuth(), tolerance);
        c.setEndGeographicPoint  ( 20,  40);  assertEquals( 170, c.getStartingAzimuth(), tolerance);
        c.setEndGeographicPoint  ( 20,  30);  assertEquals( 180, c.getStartingAzimuth(), tolerance);
        c.setEndGeographicPoint  (-20,  30);  assertEquals( 180, c.getStartingAzimuth(), tolerance);
        c.setEndGeographicPoint  (-90,  30);  assertEquals( 180, c.getStartingAzimuth(), tolerance);

        c.setStartGeographicPoint( 90,   0);
        c.setEndGeographicPoint  ( 20,  20);  assertEquals( 160, c.getStartingAzimuth(), tolerance);
        c.setEndGeographicPoint  ( 20, -20);  assertEquals(-160, c.getStartingAzimuth(), tolerance);
        c.setEndGeographicPoint  ( 20,   0);  assertEquals( 180, c.getStartingAzimuth(), tolerance);
        c.setEndGeographicPoint  (-90,   0);  assertEquals( 180, c.getStartingAzimuth(), tolerance);
    }

    /**
     * Tests geodesic distances and rhumb line length on the equator.
     */
    @Test
    public void testDistanceAtEquator() {
        final Random random = TestUtilities.createRandomNumberGenerator(86909845084528L);
        final GeodeticCalculator c = create(false);
        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 = (2*λmax) * random.nextDouble() - λmax;
            c.setEndGeographicPoint(0, x);
            final double expected = abs(x) * r;
            assertEquals("Geodesic",   expected, c.getGeodesicDistance(), Formulas.LINEAR_TOLERANCE);
            assertEquals("Rhumb line", expected, c.getRhumblineLength(),  Formulas.LINEAR_TOLERANCE);
        }
    }

    /**
     * Tests {@link GeodeticCalculator#getGeodesicDistance()} and azimuths with the example given in Wikipedia.
     * This computes the great circle route from Valparaíso (33°N 71.6W) to Shanghai (31.4°N 121.8°E).
     *
     * @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() {
        final GeodeticCalculator c = create(false);
        c.setStartGeographicPoint(-33.0, -71.6);            // Valparaíso
        c.setEndGeographicPoint  ( 31.4, 121.8);            // Shanghai
        /*
         * Wikipedia example gives:
         *
         *     Δλ = −166.6°
         *     α₁ = −94.41°
         *     α₂ = −78.42°
         *     Δσ = 168.56°   →    taking R = 6371 km, the distance is 18743 km.
         */
        assertEquals(Units.METRE,         c.getDistanceUnit());
        assertEquals("α₁",        -94.41, c.getStartingAzimuth(), 0.005);
        assertEquals("α₂",        -78.42, c.getEndingAzimuth(),   0.005);
        assertEquals("distance",   18743, c.getGeodesicDistance() / 1000, 0.5);
        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.41);
        c.setGeodesicDistance(18743000);
        assertEquals("α₁",        -94.41, c.getStartingAzimuth(), 1E-12);           // Should be the specified value.
        assertEquals("α₂",        -78.42, c.getEndingAzimuth(),   0.01);
        assertEquals("distance",   18743, c.getGeodesicDistance() / 1000, STRICT);  // Should be the specified value.
        assertPositionEquals(31.4, 121.8, c.getEndPoint(), 0.01);
    }

    /**
     * Tests geodetic calculator involving a coordinate operation.
     * This test uses a simple CRS with only the axis order interchanged.
     * The coordinates are the same than {@link #testGeodesicDistanceAndAzimuths()}.
     */
    @Test
    @DependsOnMethod("testGeodesicDistanceAndAzimuths")
    public void testUsingTransform() {
        final GeodeticCalculator c = create(true);
        assertAxisDirectionsEqual("GeographicCRS", c.getGeographicCRS().getCoordinateSystem(), AxisDirection.NORTH, AxisDirection.EAST);
        assertAxisDirectionsEqual("PositionCRS",     c.getPositionCRS().getCoordinateSystem(), AxisDirection.EAST, AxisDirection.NORTH);
        final double φ = -33.0;
        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.2);
    }

    /**
     * Tests {@link GeodeticCalculator#moveToEndPoint()}.
     */
    @Test
    public void testMoveToEndPoint() {
        // 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);
    }

    /**
     * Tests {@link GeodeticCalculator#createGeodesicCircle2D(double)}.
     */
    @Test
    @DependsOnMethod("testUsingTransform")
    public void testGeodesicCircle2D() {
        final GeodeticCalculator c = create(true);
        c.setStartGeographicPoint(-33.0, -71.6);                // Valparaíso
        c.setGeodesicDistance(100000);                          // 100 km
        Shape region = c.createGeodesicCircle2D(10000);
        if (VisualCheck.SHOW_WIDGET) {
            VisualCheck.show(region);
        }
        final Rectangle2D bounds = region.getBounds2D();
        assertEquals("x",    -71.6, bounds.getCenterX(), 1E-3);
        assertEquals("y",    -33.0, bounds.getCenterY(), 1E-3);
        assertEquals("width",  2.1, bounds.getWidth(),   0.1);
        assertEquals("height", 1.8, bounds.getHeight(),  0.1);
    }

    /**
     * Tests {@link GeodeticCalculator#createGeodesicPath2D(double)}. This method uses a CRS that swap axis order
     * as a way to verify that user-specified CRS is taken in account. The start point and end point are the same
     * than in {@link #testGeodesicDistanceAndAzimuths()}. Note that this path crosses the anti-meridian,
     * so the end point needs to be shifted by 360°.
     */
    @Test
    @DependsOnMethod("testUsingTransform")
    public void testGeodesicPath2D() {
        final GeodeticCalculator c = create(true);
        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),   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.
         */
        assertInstanceOf("Multicurves", Path2D.class, multiCurves);
        if (VisualCheck.SHOW_WIDGET) {
            VisualCheck.show(singleCurve, multiCurves);
        }
    }

    /**
     * Verifies that all <var>y</var> coordinates are zero for a geodesic path on equator.
     */
    @Test
    public void testGeodesicPathOnEquator() {
        final GeodeticCalculator c = create(false);
        final double tolerance = 1E-12;
        c.setStartGeographicPoint(0, 20);
        c.setEndGeographicPoint  (0, 12);
        assertEquals(-90, c.getStartingAzimuth(), tolerance);
        assertEquals(-90, c.getEndingAzimuth(),   tolerance);
        final Shape geodeticCurve = c.createGeodesicPath2D(1);
        final double[] coords = new double[2];
        for (final PathIterator it = geodeticCurve.getPathIterator(null, 1); !it.isDone(); it.next()) {
            it.currentSegment(coords);
            assertEquals ("φ",  0, coords[0], tolerance);
            assertBetween("λ", 12, 20, coords[1]);
        }
    }

    /**
     * Tests geodesic path between random points. The coordinates are compared with values computed by
     * <a href="https://geographiclib.sourceforge.io/">GeographicLib</a>, taken as reference implementation.
     */
    @Test
    public void testBetweenRandomPoints() {
        final Random random = TestUtilities.createRandomNumberGenerator();
        final GeodeticCalculator c = create(false);
        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;
            final double φ2 = random.nextDouble() * 180 -  90;
            final double Δλ = random.nextDouble() * 360 - 180;
            final double λ2 = IEEEremainder(λ1 + Δλ, 360);
            c.setStartGeographicPoint(φ1, λ1);
            c.setEndGeographicPoint  (φ2, λ2);
            final double geodesic  = c.getGeodesicDistance();
            final double rhumbLine = c.getRhumblineLength();
            final GeodesicData expected = reference.Inverse(φ1, λ1, φ2, λ2);
            assertEquals("Geodesic distance", expected.s12,  geodesic,               Formulas.LINEAR_TOLERANCE);
            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 (VERBOSE) {
                // Checks the geodesic path on only 10% of test data, because this computation is expensive.
                if ((i % 10) == 0) {
                    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 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) {
        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 double       azimuth  = c.getStartingAzimuth();
        final double       toMetres = (PI/180) * c.semiMajorAxis;
        final double[]     buffer   = new double[2];
        while (!iterator.isDone()) {
            switch (iterator.currentSegment(buffer)) {
                default: fail("Unexpected segment"); break;
                case PathIterator.SEG_MOVETO: break;
                case PathIterator.SEG_LINETO: {
                    c.setEndGeographicPoint(buffer[0], buffer[1]);
                    aErrors.accept(abs(c.getStartingAzimuth() - azimuth));
                    c.setStartingAzimuth(azimuth);
                    DirectPosition endPoint = c.getEndPoint();
                    final double φ = endPoint.getOrdinate(0);
                    final double λ = endPoint.getOrdinate(1);
                    double dy =              (buffer[0] - φ)      * toMetres;
                    double dx = IEEEremainder(buffer[1] - λ, 360) * toMetres * cos(toRadians(φ));
                    yError.accept(abs(dy) / resolution);
                    xError.accept(abs(dx) / resolution);
                }
            }
            iterator.next();
        }
        return new Statistics[] {xError, yError, aErrors};
    }

    /**
     * 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.
     *
     * @throws IOException if an error occurred while reading the test file.
     */
    @Test
    public void compareAgainstDataset() throws IOException {
        int noConvergenceCount = 0;
        try (LineNumberReader reader = OptionalTestData.GEODESIC.reader()) {
            final Random random = TestUtilities.createRandomNumberGenerator();
            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.semiMajorAxis);
                    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 {
                    /*
                     * When ellipsoidal formulas are used, we aim for an 1 mm accuracy in coordinate values.
                     * We also aim for azimuths 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.
                     */
                    linearTolerance    = 2 * GeodesicsOnEllipsoid.ITERATION_TOLERANCE * AUTHALIC_RADIUS;
                    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.
                        }
                    }
                }
                /*
                 * 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 (GeodeticException | 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);
    }

    /**
     * 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.
     */
    boolean isFailure(final double[] expected) {
        return true;
    }

    /**
     * Tests {@link GeodesicsOnEllipsoid#getRhumblineLength()} using points given by Bennett (1996).
     * This is an anti-regression test since the result was computed by SIS for the spherical case.
     */
    @Test
    public void testRhumblineLength() {
        final GeodeticCalculator c = create(false);
        c.setStartGeographicPoint(10+18.4/60,  37+41.7/60);
        c.setEndGeographicPoint  (53+29.5/60, 113+17.1/60);
        assertEquals(8344561, c.getRhumblineLength(), 1);
        assertEquals(54.8682, c.getConstantAzimuth(), 1E-4);
    }

    /**
     * Tests {@link GeodesicsOnEllipsoid#getRhumblineLength()} using points given by Bennett (1996).
     * This is an anti-regression test since the result was computed by SIS for the spherical case.
     */
    @Test
    public void testRhumblineNearlyEquatorial() {
        final GeodeticCalculator c = create(false);
        c.setStartGeographicPoint(-52-47.8/60, -97-31.6/60);
        c.setEndGeographicPoint  (-53-10.8/60, -41-34.6/60);
        assertEquals(3745332, c.getRhumblineLength(), 1);
        assertEquals(90.6521, c.getConstantAzimuth(), 1E-4);
    }

    /**
     * Tests {@link GeodesicsOnEllipsoid#getRhumblineLength()} using points given by Bennett (1996).
     * This is an anti-regression test since the result was computed by SIS for the spherical case.
     */
    @Test
    public void testRhumblineEquatorial() {
        final GeodeticCalculator c = create(false);
        c.setStartGeographicPoint(48+45.0/60, -61-31.1/60);
        c.setEndGeographicPoint  (48+45.0/60,   5+13.2/60);
        assertEquals(4892987, c.getRhumblineLength(), 1);
        assertEquals(90, c.getConstantAzimuth(), STRICT);
    }
}