id,summary,reporter,owner,description,type,status,milestone,component,version,severity,resolution,keywords,cc
11930,huiller method is inaccurate,Charles Karney <charles@…>,Barend Gehrels,"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.

{{{#!c++
#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
  [https://en.wikipedia.org/wiki/Spherical_trigonometry#Area_and_spherical_excess 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
  [http://geographiclib.sourceforge.net/html/ GeographicLib] for
  computing geodesic areas.
",Bugs,closed,Boost 1.64.0,geometry,Boost 1.58.0,Problem,fixed,area,charles.karney@…