Почему буст :: геометрия географического расстояния Винсенти неточна вокруг экватора?

Мне нужна функция для расчета расстояния между парой WGS 84 позиции с высокой степенью точности, и я планировал использовать geographic функции от повысить геометрию.

повысить геометрию Дизайн Rational состояния:

Есть метод Андойера, быстрый и точный, и метод Винсенти, более медленный и более точный.

Однако при тестировании boost::geometry::distance функция как с Andoyer а также Vincenty Стратегии, я получил следующие результаты:

WGS 84 values (metres)
Semimajor axis:         6378137.000000
Flattening:             0.003353
Semiminor axis:         6356752.314245

Semimajor distance:     20037508.342789
Semiminor distance:     19970326.371123

Boost geometry near poles
Andoyer function:
Semimajor distance:     20037508.151445
Semiminor distance:     20003917.164970
Vincenty function:
Semimajor distance:     **19970326.180419**
Semiminor distance:     20003931.266635

Boost geometry at poles
Andoyer function:
Semimajor distance:     0.000000
Semiminor distance:     0.000000
Vincenty function:
Semimajor distance:     **19970326.371122**
Semiminor distance:     20003931.458623

Vincenty расстояния вдоль большой оси (то есть вокруг экватора) меньше, чем расстояние вокруг оси полуминора между северным и южным полюсами. Это не может быть правильно.

Полуминор и Andoyer расстояния выглядят разумно. За исключением случаев, когда точки находятся на противоположной стороне Земли, когда boost Andoyer функция возвращает ноль!

Проблема в: Vincenty алгоритм, boost geometry реализация этого или мой тестовый код?

Тестовый код:

/// boost geometry WGS84 distance issue

// Note: M_PI is not part of the C or C++ standards, _USE_MATH_DEFINES enables it
#define _USE_MATH_DEFINES
#include <boost/geometry.hpp>
#include <cmath>
#include <iostream>
#include <ios>

// WGS 84 parameters from: Eurocontrol WGS 84 Implementation Manual
// Version 2.4 Chapter 3, page 14

/// The Semimajor axis measured in metres.
/// This is the radius at the equator.
constexpr double a = 6378137.0;

/// Flattening, a ratio.
/// This is the flattening of the ellipse at the poles
constexpr double f = 1.0/298.257223563;

/// The Semiminor axis measured in metres.
/// This is the radius at the poles.
/// Note: this is derived from the Semimajor axis and the flattening.
/// See WGS 84 Implementation Manual equation B-2, page 69.
constexpr double b = a * (1.0 - f);

int main(int /*argc*/, char ** /*argv*/)
{
std::cout.setf(std::ios::fixed);

std::cout << "WGS 84 values (metres)\n";
std::cout << "\tSemimajor axis:\t\t"   << a << "\n";
std::cout << "\tFlattening:\t\t"       << f << "\n";
std::cout << "\tSemiminor axis:\t\t"   << b << "\n\n";

std::cout << "\tSemimajor distance:\t" << M_PI * a << "\n";
std::cout << "\tSemiminor distance:\t" << M_PI * b << "\n";
std::cout << std::endl;

// Min value for delta. 0.000000014 causes Andoyer to fail.
const double DELTA(0.000000015);

// For boost::geometry:
typedef boost::geometry::cs::geographic<boost::geometry::radian> Wgs84Coords;
typedef boost::geometry::model::point<double, 2, Wgs84Coords> GeographicPoint;
// Note boost points are Long & Lat NOT Lat & Long
GeographicPoint near_north_pole   (0.0,  M_PI_2 - DELTA);
GeographicPoint near_south_pole   (0.0, -M_PI_2 + DELTA);

GeographicPoint near_equator_east ( M_PI_2 - DELTA, 0.0);
GeographicPoint near_equator_west (-M_PI_2 + DELTA, 0.0);

// Note: the default boost geometry spheroid is WGS84
// #include <boost/geometry/core/srs.hpp>
typedef boost::geometry::srs::spheroid<double> SpheroidType;
SpheroidType spheriod;

//#include <boost/geometry/strategies/geographic/distance_andoyer.hpp>
typedef boost::geometry::strategy::distance::andoyer<SpheroidType>
AndoyerStrategy;
AndoyerStrategy andoyer(spheriod);

std::cout << "Boost geometry near poles\n";
std::cout << "Andoyer function:\n";
double andoyer_major(boost::geometry::distance(near_equator_east, near_equator_west, andoyer));
std::cout << "\tSemimajor distance:\t" << andoyer_major << "\n";
double andoyer_minor(boost::geometry::distance(near_north_pole, near_south_pole, andoyer));
std::cout << "\tSemiminor distance:\t" << andoyer_minor << "\n";

//#include <boost/geometry/strategies/geographic/distance_vincenty.hpp>
typedef boost::geometry::strategy::distance::vincenty<SpheroidType>
VincentyStrategy;
VincentyStrategy vincenty(spheriod);

std::cout << "Vincenty function:\n";
double vincenty_major(boost::geometry::distance(near_equator_east, near_equator_west, vincenty));
std::cout << "\tSemimajor distance:\t" << vincenty_major << "\n";
double vincenty_minor(boost::geometry::distance(near_north_pole, near_south_pole, vincenty));
std::cout << "\tSemiminor distance:\t" << vincenty_minor << "\n\n";

// Note boost points are Long & Lat NOT Lat & Long
GeographicPoint north_pole   (0.0,  M_PI_2);
GeographicPoint south_pole   (0.0, -M_PI_2);

GeographicPoint equator_east ( M_PI_2, 0.0);
GeographicPoint equator_west (-M_PI_2, 0.0);

std::cout << "Boost geometry at poles\n";
std::cout << "Andoyer function:\n";
andoyer_major = boost::geometry::distance(equator_east, equator_west, andoyer);
std::cout << "\tSemimajor distance:\t" << andoyer_major << "\n";
andoyer_minor = boost::geometry::distance(north_pole, south_pole, andoyer);
std::cout << "\tSemiminor distance:\t" << andoyer_minor << "\n";

std::cout << "Vincenty function:\n";
vincenty_major = boost::geometry::distance(equator_east, equator_west, vincenty);
std::cout << "\tSemimajor distance:\t" << vincenty_major << "\n";
vincenty_minor = boost::geometry::distance(north_pole, south_pole, vincenty);
std::cout << "\tSemiminor distance:\t" << vincenty_minor << "\n";

return 0;
}