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.
This archive was generated by hypermail 2.1.7 : 2017-02-16 18:50:19 UTC