title: Why is boost::geometry geographic Vincenty distance inaccurate around the Equator?
body: <p>I require a function to calculate the distance between a pair of <a href="http://www.icao.int/safety/pbn/documentation/eurocontrol/eurocontrol%20wgs%2084%20implementation%20manual.pdf" rel="noreferrer">WGS 84</a> positions to a high degree of accuracy and I was planning to use the <code>geographic</code> functions from <a href="http://www.boost.org/doc/libs/1_60_0/libs/geometry/doc/html/index.html" rel="noreferrer">boost geometry</a>.</p>
  
  <p>The <a href="http://www.boost.org/doc/libs/1_60_0/libs/geometry/doc/html/geometry/design.html" rel="noreferrer">boost geometry Design Rational</a> states:</p>
  
  <blockquote>
    <p>There is the Andoyer method, fast and precise, and there is the Vincenty method, slower and more precise..</p>
  </blockquote>
  
  <p>However, when testing the <code>boost::geometry::distance</code> function with both the <code>Andoyer</code> and <code>Vincenty</code> strategies, I got the following results:</p>
  
  <pre><code>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
  </code></pre>
  
  <p>The <code>Vincenty</code> distances along the Semimajor axis (i.e. around the Equator) are less than the distance around the Semiminor axis between the North and South Poles. That can't be correct.</p>
  
  <p>The Semiminor  and <code>Andoyer</code> distances look reasonable. Except when the points are on opposing side of the Earth, when the <code>boost</code> <code>Andoyer</code> function returns zero!</p>
  
  <p>Is the problem in: the <code>Vincenty</code> algorithm, the <code>boost geometry</code> implementation of it, or my test code?  </p>
  
  <p>Test code:</p>
  
  <pre><code>/// 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 &lt;boost/geometry.hpp&gt;
  #include &lt;cmath&gt;
  #include &lt;iostream&gt;
  #include &lt;ios&gt;
  
  // 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 &lt;&lt; "WGS 84 values (metres)\n";
    std::cout &lt;&lt; "\tSemimajor axis:\t\t"   &lt;&lt; a &lt;&lt; "\n";
    std::cout &lt;&lt; "\tFlattening:\t\t"       &lt;&lt; f &lt;&lt; "\n";
    std::cout &lt;&lt; "\tSemiminor axis:\t\t"   &lt;&lt; b &lt;&lt; "\n\n";
  
    std::cout &lt;&lt; "\tSemimajor distance:\t" &lt;&lt; M_PI * a &lt;&lt; "\n";
    std::cout &lt;&lt; "\tSemiminor distance:\t" &lt;&lt; M_PI * b &lt;&lt; "\n";
    std::cout &lt;&lt; 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&lt;boost::geometry::radian&gt; Wgs84Coords;
    typedef boost::geometry::model::point&lt;double, 2, Wgs84Coords&gt; GeographicPoint;
    // Note boost points are Long &amp; Lat NOT Lat &amp; 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 &lt;boost/geometry/core/srs.hpp&gt;
    typedef boost::geometry::srs::spheroid&lt;double&gt; SpheroidType;
    SpheroidType spheriod;
  
    //#include &lt;boost/geometry/strategies/geographic/distance_andoyer.hpp&gt;
    typedef boost::geometry::strategy::distance::andoyer&lt;SpheroidType&gt;
                                                                 AndoyerStrategy;
    AndoyerStrategy andoyer(spheriod);
  
    std::cout &lt;&lt; "Boost geometry near poles\n";
    std::cout &lt;&lt; "Andoyer function:\n";
    double andoyer_major(boost::geometry::distance(near_equator_east, near_equator_west, andoyer));
    std::cout &lt;&lt; "\tSemimajor distance:\t" &lt;&lt; andoyer_major &lt;&lt; "\n";
    double andoyer_minor(boost::geometry::distance(near_north_pole, near_south_pole, andoyer));
    std::cout &lt;&lt; "\tSemiminor distance:\t" &lt;&lt; andoyer_minor &lt;&lt; "\n";
  
    //#include &lt;boost/geometry/strategies/geographic/distance_vincenty.hpp&gt;
    typedef boost::geometry::strategy::distance::vincenty&lt;SpheroidType&gt;
                                                                 VincentyStrategy;
    VincentyStrategy vincenty(spheriod);
  
    std::cout &lt;&lt; "Vincenty function:\n";
    double vincenty_major(boost::geometry::distance(near_equator_east, near_equator_west, vincenty));
    std::cout &lt;&lt; "\tSemimajor distance:\t" &lt;&lt; vincenty_major &lt;&lt; "\n";
    double vincenty_minor(boost::geometry::distance(near_north_pole, near_south_pole, vincenty));
    std::cout &lt;&lt; "\tSemiminor distance:\t" &lt;&lt; vincenty_minor &lt;&lt; "\n\n";
  
    // Note boost points are Long &amp; Lat NOT Lat &amp; 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 &lt;&lt; "Boost geometry at poles\n";
    std::cout &lt;&lt; "Andoyer function:\n";
    andoyer_major = boost::geometry::distance(equator_east, equator_west, andoyer);
    std::cout &lt;&lt; "\tSemimajor distance:\t" &lt;&lt; andoyer_major &lt;&lt; "\n";
    andoyer_minor = boost::geometry::distance(north_pole, south_pole, andoyer);
    std::cout &lt;&lt; "\tSemiminor distance:\t" &lt;&lt; andoyer_minor &lt;&lt; "\n";
  
    std::cout &lt;&lt; "Vincenty function:\n";
    vincenty_major = boost::geometry::distance(equator_east, equator_west, vincenty);
    std::cout &lt;&lt; "\tSemimajor distance:\t" &lt;&lt; vincenty_major &lt;&lt; "\n";
    vincenty_minor = boost::geometry::distance(north_pole, south_pole, vincenty);
    std::cout &lt;&lt; "\tSemiminor distance:\t" &lt;&lt; vincenty_minor &lt;&lt; "\n";
  
    return 0;
  }
  </code></pre>