Karney gives often wrong results

[barendgehrels]

The Karney formula (inside Boost.Geometry) gives often very wrong results.

(In GeographicLib, the results are different - and correct)

It fails in several cases:

• Results: negative distances
• Output ≈ D_antipodal (~20003 km) for pairs that are nowhere near antipodal
• Off by exactly one quarter-meridian (~10001 km) for near-meridional polar-transit pairs
• ~25% relative error on short (~km-scale) pairs entirely above 45° latitude

The bug is not limited to edge cases — case 7 in the repro is a 6.5 km pair near Rotterdam (both endpoints at ~51° N) where our version of Karney returns 8009 m instead of 6520 m.

Minimal reproduction scenario:

#include <iostream>
#include <iomanip>

#include <boost/geometry/geometry.hpp>

namespace bg = boost::geometry;

using latlon_point = bg::model::point<double, 2, bg::cs::geographic<bg::degree>>;

struct distance_pair
{
    latlon_point first;
    latlon_point second;
    double geodesic_distance{0};
};

int main()
{
    using stype = bg::srs::spheroid<double>;

    std::vector<distance_pair> pairs;
    pairs.push_back({{0.0, 0.0}, {2.0, 0.0}, 222638.98}); // OK

    pairs.push_back({{0.0, 41.306061250198}, {100.88904290879072632, -52.570275453324655475}, 14174344.4466144}); // source: GeodTest.dat line 78768
    pairs.push_back({{0.0, 25.7158673231}, {32.9180106939834239, -66.243826410425021118}, 10569765.894952}); // source: GeodTest.dat line 78769
    pairs.push_back({{0.0, 58.256371249077}, {22.784150207906535006, 28.739702371452295937},  3721028.7556286}); // source: GeodTest.dat line 78770
    pairs.push_back({{0.0, 17.812015912354}, {179.998636794719653263, 45.432345524206019945}, 13000761.9422694}); // source: GeodTest.dat line 349760
    pairs.push_back({{0.0, 64.858589732434}, { 0.002464379937107867, -50.63203940688252426}, 12806725.1078506}); // source: GeodTest.dat line 349762

    pairs.push_back({{4.19132,51.1189}, {4.1005,51.106}, 6520.61}); // source: local random test

    bg::strategy::distance::karney<stype> karney{};
    bg::strategy::distance::vincenty<stype> vincenty{};

    for (auto const& item : pairs)
    {
        auto const k = bg::distance(item.first, item.second, karney);
        auto const v = bg::distance(item.first, item.second, vincenty);
        std::cout << "Karney: " << std::setprecision(10) << std::setw(12) << k
                  << "  Vincenty: " << std::setw(12) << v
                  << "  expected: " << std::setw(12) << item.geodesic_distance 
                  << std::endl;
    }

    return 0;
}

Giving:

Karney:  222638.9816  Vincenty:  222638.9816  expected:    222638.98
Karney: -5587.447352  Vincenty:  14174344.45  expected:  14174344.45
Karney: -3120.426059  Vincenty:  10569765.89  expected:  10569765.89
Karney:  19988085.62  Vincenty:  3721028.756  expected:  3721028.756
Karney:  20001692.52  Vincenty:  13000761.94  expected:  13000761.94
Karney:  2829433.304  Vincenty:  12806725.11  expected:  12806725.11
Karney:  8008.921612  Vincenty:  6519.839481  expected:      6520.61

Vincenty is not perfect. But Karney should be perfect and give identical results to the version in GeographicLib. However, it is so often off that it is a mystery what is the exact cause, and why it is not reported before...

This came from the investigation of #1461