[Boost-bugs] [Boost C++ Libraries] #11930: huiller method is inaccurate

Subject: [Boost-bugs] [Boost C++ Libraries] #11930: huiller method is inaccurate
From: Boost C++ Libraries (noreply_at_[hidden])
Date: 2016-01-21 19:05:41

#11930: huiller method is inaccurate
----------------------------------------+---------------------------
 Reporter:  Charles Karney <charles@…>  |      Owner:  barendgehrels
     Type:  Bugs                        |     Status:  new
Milestone:  To Be Determined            |  Component:  geometry
  Version:  Boost 1.58.0                |   Severity:  Problem
 Keywords:  area                        |
----------------------------------------+---------------------------
 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.

-- 
Ticket URL: <https://svn.boost.org/trac/boost/ticket/11930>
Boost C++ Libraries <http://www.boost.org/>
Boost provides free peer-reviewed portable C++ source libraries.