huiller method is inaccurate

[Charles Karney <charles@…>]

The following code illustrates a problem with the "huiller" strategy for computing spherical areas. It computes the area of a sequence of progressively skinny triangles. The relative errors in the computed areas are much larger than they need be.

#include <iostream>
#include <cmath>
#include <boost/geometry.hpp>

namespace bg = boost::geometry;
typedef bg::model::point<double, 2,
                         bg::cs::spherical_equatorial<bg::degree> > pt;

int main() {
  bg::model::polygon<pt, false, false> poly;
  bg::read_wkt("POLYGON((0 45,0 -45,0 0))", poly);
  double c = 0.014458780939650811;
  std::cout << "  area        relerr" << "\n";
  for (int i = 3; i <= 6; ++i) {
    double h = std::pow(10.0, -i);
    poly.outer()[2].set<0>(h);
    double
      area = bg::area(poly),
      // Relative error in this estimate is 2e-11 for h = 0.001
      area0 = h * c;
    std::cout << area << " " << (area - area0)/(area0) << "\n";
  }
}

Running this code gives

  area        relerr
1.44588e-05 -2.0871e-06
1.44594e-06 4.55941e-05
1.43961e-07 -0.00433739
0 -1

Discussion:

• Surely it's a bad idea to embed the algorithm name "huiller" into the API. Computing spherical areas is really straightforward. Users shouldn't need to worry about the underlying method. Maybe "spherical_area" is a better name.
• In any case, you've misspelled the guy's name. It's either "l'Huilier" or "l'Huillier" with two "i"s.
• Furthermore, you're needlessly blackening l'Huilier's name by applying his formula for the area of a spherical triangle in terms of its sides.
• The three-side formula is a terrible choice because if a + b is nearly equal to c (as is the case in my example), the triangle is badly condititioned.
• Much better is to use the formula of the area of a triangle in terms of two sides and the included angle. See the last two formulas in the Wikipedia article on ​spherical excess. This nicely reduces to the familar trapezium rule in the plane limit.
• This formula comes with a sign attached, so there's no need to kludge a sign on afterwards as there is with the current implementation.
• This formula together with an ellipsoidal correction is used by ​GeographicLib for computing geodesic areas.

[awulkiew]

Thanks!

Fixed by this PR: ​https://github.com/boostorg/geometry/pull/359

The name is changed to area::spherical.

The new results are:

  area        relerr
1.44588e-05 2.10295e-11
1.44588e-06 2.10165e-13
1.44588e-07 2.01378e-15
1.44588e-08 0