https://svn.boost.org/trac10/ticket/11929 <p> The following code illustrates a problem with the "haversine" strategy for computing great circle distances. It computes the distance between two nearly antipodal points. The relative error in the result is 5e-9, much bigger than you should get with doubles. </p> <div class="wiki-code"><div class="code"><pre><span class="cp">#include</span> <span class="cpf">&lt;iostream&gt;</span><span class="cp"></span> <span class="cp">#include</span> <span class="cpf">&lt;cmath&gt;</span><span class="cp"></span> <span class="cp">#include</span> <span class="cpf">&lt;boost/geometry.hpp&gt;</span><span class="cp"></span> <span class="k">namespace</span> <span class="n">bg</span> <span class="o">=</span> <span class="n">boost</span><span class="o">::</span><span class="n">geometry</span><span class="p">;</span> <span class="k">typedef</span> <span class="n">bg</span><span class="o">::</span><span class="n">model</span><span class="o">::</span><span class="n">point</span><span class="o">&lt;</span><span class="kt">double</span><span class="p">,</span> <span class="mi">2</span><span class="p">,</span> <span class="n">bg</span><span class="o">::</span><span class="n">cs</span><span class="o">::</span><span class="n">spherical_equatorial</span><span class="o">&lt;</span><span class="n">bg</span><span class="o">::</span><span class="n">degree</span><span class="o">&gt;</span> <span class="o">&gt;</span> <span class="n">pt</span><span class="p">;</span> <span class="kt">int</span> <span class="nf">main</span><span class="p">()</span> <span class="p">{</span> <span class="kt">double</span> <span class="n">pi</span> <span class="o">=</span> <span class="n">std</span><span class="o">::</span><span class="n">atan2</span><span class="p">(</span><span class="mi">0</span><span class="p">,</span> <span class="o">-</span><span class="mi">1</span><span class="p">);</span> <span class="n">bg</span><span class="o">::</span><span class="n">model</span><span class="o">::</span><span class="n">segment</span><span class="o">&lt;</span><span class="n">pt</span><span class="o">&gt;</span> <span class="n">seg</span><span class="p">;</span> <span class="n">bg</span><span class="o">::</span><span class="n">read_wkt</span><span class="p">(</span><span class="s">&quot;SEGMENT(0.000001 0, 180 0)&quot;</span><span class="p">,</span> <span class="n">seg</span><span class="p">);</span> <span class="kt">double</span> <span class="c1">// arc length (converted to degrees</span> <span class="n">dist</span> <span class="o">=</span> <span class="n">bg</span><span class="o">::</span><span class="n">distance</span><span class="p">(</span><span class="n">seg</span><span class="p">.</span><span class="n">first</span><span class="p">,</span> <span class="n">seg</span><span class="p">.</span><span class="n">second</span><span class="p">)</span> <span class="o">*</span> <span class="p">(</span><span class="mi">180</span><span class="o">/</span><span class="n">pi</span><span class="p">),</span> <span class="c1">// the correct distance (in degrees)</span> <span class="n">dist0</span> <span class="o">=</span> <span class="mi">180</span> <span class="o">-</span> <span class="mf">0.000001</span><span class="p">,</span> <span class="n">err</span> <span class="o">=</span> <span class="p">(</span><span class="n">dist</span> <span class="o">-</span> <span class="n">dist0</span><span class="p">)</span> <span class="o">/</span> <span class="n">dist0</span><span class="p">;</span> <span class="n">std</span><span class="o">::</span><span class="n">cout</span> <span class="o">&lt;&lt;</span> <span class="s">&quot;relative error &quot;</span> <span class="o">&lt;&lt;</span> <span class="n">err</span> <span class="o">&lt;&lt;</span> <span class="s">&quot;</span><span class="se">\n</span><span class="s">&quot;</span><span class="p">;</span> <span class="p">}</span> </pre></div></div><p> Running this code gives </p> <pre class="wiki">relative error 5.55556e-09 </pre><p> Discussion: </p> <ul><li>Surely it's a bad idea to embed the algorithm name "haversine" into the API. Computing great circle distances is really straightforward. Users shouldn't need to worry about the underlying method. Maybe "great_circle" is a better name. </li><li>The use of haversines made sense when the log haversine function was tabulated and you needed trigonometric identities involving products to facilitate computations with log tables. Nowadays there are better formulas. </li><li>As this example illustrates, it's subject to a loss of accuracy with antipodal points (it takes the arcsine of a number close to 1.) </li><li>A good formula to use is given in the Wikipedia article on <a class="ext-link" href="https://en.wikipedia.org/wiki/Great-circle_navigation#cite_note-2"><span class="icon">​</span>great-circle navigation</a> (this is the formula used by Vincenty in the spherical limit -- but it surely doesn't originate with him). </li></ul> en-us /htdocs/site/boost.png https://svn.boost.org/trac10/ticket/11929 Trac 1.4.3 Charles Karney <charles@…> Thu, 21 Jan 2016 19:03:36 GMT https://svn.boost.org/trac10/ticket/11929#comment:1 https://svn.boost.org/trac10/ticket/11929#comment:1 <ul> <li><strong>summary</strong> <span class="trac-field-old">haversine method is accurate</span> → <span class="trac-field-new">haversine method is inaccurate</span> </li> </ul> Ticket