Formula for meridian distance in geography

include/boost/geometry/formulas/andoyer_inverse.hpp

@@ -214,7 +214,7 @@ class andoyer_inverse
     static inline void normalize_azimuth(CT & azimuth, CT const& A, CT const& dA)
     {
         CT const c0 = 0;
-        
+
         if (A >= c0) // A indicates Eastern hemisphere
         {
             if (dA >= c0) // A altered towards 0

include/boost/geometry/formulas/elliptic_arc_length.hpp

@@ -0,0 +1,166 @@
+// Boost.Geometry
+
+// Copyright (c) 2017 Oracle and/or its affiliates.
+
+// Contributed and/or modified by Vissarion Fysikopoulos, on behalf of Oracle
+
+// Use, modification and distribution is subject to the Boost Software License,
+// Version 1.0. (See accompanying file LICENSE_1_0.txt or copy at
+// http://www.boost.org/LICENSE_1_0.txt)
+
+#ifndef BOOST_GEOMETRY_FORMULAS_ELLIPTIC_ARC_LENGTH_HPP
+#define BOOST_GEOMETRY_FORMULAS_ELLIPTIC_ARC_LENGTH_HPP
+
+
+#include <boost/math/constants/constants.hpp>
+
+#include <boost/geometry/core/radius.hpp>
+#include <boost/geometry/core/srs.hpp>
+
+#include <boost/geometry/util/condition.hpp>
+#include <boost/geometry/util/math.hpp>
+
+#include <boost/geometry/formulas/flattening.hpp>
+
+
+namespace boost { namespace geometry { namespace formula
+{
+
+/*!
+\brief Compute the arc length of an ellipse.
+*/
+
+template <typename CT, unsigned int Order = 1>
+class elliptic_arc_length
+{
+
+public :
+
+    // Distance computation on meridians using series approximations
+    // to elliptic integrals. Formula to compute distance from lattitude 0 to lat
+    // https://en.wikipedia.org/wiki/Meridian_arc
+    // latitudes are assumed to be in radians and in [-pi/2,pi/2]
+    template <typename T, typename Spheroid>
+    static CT apply(T lat, Spheroid const& spheroid)
+    {
+        CT const a = get_radius<0>(spheroid);
+        CT const f = formula::flattening<CT>(spheroid);
+        CT n = f / (CT(2) - f);
+        CT M = a/(1+n);
+        CT C0 = 1;
+
+        if (Order == 0)
+        {
+           return M * C0 * lat;
+        }
+
+        CT C2 = -1.5 * n;
+
+        if (Order == 1)
+        {
+            return M * (C0 * lat + C2 * sin(2*lat));
+        }
+
+        CT n2 = n * n;
+        C0 += .25 * n2;
+        CT C4 = 0.9375 * n2;
+
+        if (Order == 2)
+        {
+            return M * (C0 * lat + C2 * sin(2*lat) + C4 * sin(4*lat));
+        }
+
+        CT n3 = n2 * n;
+        C2 += 0.1875 * n3;
+        CT C6 = -0.729166667 * n3;
+
+        if (Order == 3)
+        {
+            return M * (C0 * lat + C2 * sin(2*lat) + C4 * sin(4*lat)
+                      + C6 * sin(6*lat));
+        }
+
+        CT n4 = n2 * n2;
+        C4 -= 0.234375 * n4;
+        CT C8 = 0.615234375 * n4;
+
+        if (Order == 4)
+        {
+            return M * (C0 * lat + C2 * sin(2*lat) + C4 * sin(4*lat)
+                      + C6 * sin(6*lat) + C8 * sin(8*lat));
+        }
+
+        CT n5 = n4 * n;
+        C6 += 0.227864583 * n5;
+        CT C10 = -0.54140625 * n5;
+
+        // Order 5 or higher
+        return M * (C0 * lat + C2 * sin(2*lat) + C4 * sin(4*lat)
+                  + C6 * sin(6*lat) + C8 * sin(8*lat) + C10 * sin(10*lat));
+
+    }
+
+    // Iterative method to elliptic arc length based on
+    // http://www.codeguru.com/cpp/cpp/algorithms/article.php/c5115/
+    // Geographic-Distance-and-Azimuth-Calculations.htm
+    // latitudes are assumed to be in radians and in [-pi/2,pi/2]
+    template <typename T1, typename T2, typename Spheroid>
+    CT interative_method(T1 lat1,
+                         T2 lat2,
+                         Spheroid const& spheroid)
+    {
+        CT result = 0;
+        CT const zero = 0;
+        CT const one = 1;
+        CT const c1 = 2;
+        CT const c2 = 0.5;
+        CT const c3 = 4000;
+
+        CT const a = get_radius<0>(spheroid);
+        CT const f = formula::flattening<CT>(spheroid);
+
+        // how many steps to use
+
+        CT lat1_deg = lat1 * geometry::math::r2d<CT>();
+        CT lat2_deg = lat2 * geometry::math::r2d<CT>();
+
+        int steps = c1 + (c2 + (lat2_deg > lat1_deg) ? CT(lat2_deg - lat1_deg)
+                                                     : CT(lat1_deg - lat2_deg));
+        steps = (steps > c3) ? c3 : steps;
+
+        //std::cout << "Steps=" << steps << std::endl;
+
+        CT snLat1 = sin(lat1);
+        CT snLat2 = sin(lat2);
+        CT twoF   = 2 * f - f * f;
+
+        // limits of integration
+        CT x1 = a * cos(lat1) /
+                sqrt(1 - twoF * snLat1 * snLat1);
+        CT x2 = a * cos(lat2) /
+                sqrt(1 - twoF * snLat2 * snLat2);
+
+        CT dx = (x2 - x1) / (steps - one);
+        CT x, y1, y2, dy, dydx;
+        CT adx = (dx < zero) ? -dx : dx;    // absolute value of dx
+
+        CT a2 = a * a;
+        CT oneF = 1 - f;
+
+        // now loop through each step adding up all the little
+        // hypotenuses
+        for (int i = 0; i < (steps - 1); i++){
+            x = x1 + dx * i;
+            dydx = ((a * oneF * sqrt((one - ((x+dx)*(x+dx))/a2))) -
+                    (a * oneF * sqrt((one - (x*x)/a2)))) / dx;
+            result += adx * sqrt(one + dydx*dydx);
+        }
+
+        return result;
+    }
+};
+
+}}} // namespace boost::geometry::formula
+
+
+#endif // BOOST_GEOMETRY_FORMULAS_ELLIPTIC_ARC_LENGTH_HPP

include/boost/geometry/strategies/geographic/distance.hpp

@@ -5,6 +5,7 @@
 // This file was modified by Oracle on 2014-2017.
 // Modifications copyright (c) 2014-2017 Oracle and/or its affiliates.
 
+// Contributed and/or modified by Vissarion Fysikopoulos, on behalf of Oracle
 // Contributed and/or modified by Adam Wulkiewicz, on behalf of Oracle
 
 // Use, modification and distribution is subject to the Boost Software License,
@@ -28,6 +29,7 @@
 #include <boost/geometry/strategies/geographic/parameters.hpp>
 
 #include <boost/geometry/util/math.hpp>
+#include <boost/geometry/util/normalize_spheroidal_coordinates.hpp>
 #include <boost/geometry/util/promote_floating_point.hpp>
 #include <boost/geometry/util/select_calculation_type.hpp>
 
@@ -74,13 +76,52 @@ public :
     inline typename calculation_type<Point1, Point2>::type
     apply(Point1 const& point1, Point2 const& point2) const
     {
+        typedef typename calculation_type<Point1, Point2>::type CT;
+
+        CT lon1 = get_as_radian<0>(point1);
+        CT lat1 = get_as_radian<1>(point1);
+        CT lon2 = get_as_radian<0>(point2);
+        CT lat2 = get_as_radian<1>(point2);
+
+        CT c0 = 0;
+        CT pi = math::pi<CT>();
+        CT half_pi = pi/CT(2);
+        CT diff = math::longitude_distance_signed<geometry::radian>(lon1, lon2);
+
+        typedef typename formula::elliptic_arc_length
+        <
+            CT, strategy::default_order<FormulaPolicy>::value
+        > elliptic_arc_length;
+
+        if (math::equals(diff, c0))
+        {
+            // single meridian not crossing pole
+            if (lat1 > lat2)
+            {
+                std::swap(lat1, lat2);
+            }
+            return elliptic_arc_length::apply(lat2, m_spheroid)
+                  - elliptic_arc_length::apply(lat1, m_spheroid);
+        }
+
+        if (math::equals(math::abs(diff), pi))
+        {
+            // meridian crosses pole
+            CT lat_sign = 1;
+            if (lat1+lat2 < c0)
+            {
+                lat_sign = CT(-1);
+            }
+            return math::abs(lat_sign * CT(2) *
+                             elliptic_arc_length::apply(half_pi, m_spheroid)
+                             - elliptic_arc_length::apply(lat1, m_spheroid)
+                             - elliptic_arc_length::apply(lat2, m_spheroid));
+        }
+
         return FormulaPolicy::template inverse
-            <
-                typename calculation_type<Point1, Point2>::type,
-                true, false, false, false, false
-            >::apply(get_as_radian<0>(point1), get_as_radian<1>(point1),
-                     get_as_radian<0>(point2), get_as_radian<1>(point2),
-                     m_spheroid).distance;
+               <
+                   CT, true, false, false, false, false
+               >::apply(lon1, lat1, lon2, lat2, m_spheroid).distance;
     }
 
     inline Spheroid const& model() const

include/boost/geometry/strategies/geographic/parameters.hpp

@@ -12,10 +12,10 @@
 
 
 #include <boost/geometry/formulas/andoyer_inverse.hpp>
-#include <boost/geometry/formulas/thomas_inverse.hpp>
-#include <boost/geometry/formulas/vincenty_inverse.hpp>
 #include <boost/geometry/formulas/thomas_direct.hpp>
+#include <boost/geometry/formulas/thomas_inverse.hpp>
 #include <boost/geometry/formulas/vincenty_direct.hpp>
+#include <boost/geometry/formulas/vincenty_inverse.hpp>
 
 #include <boost/mpl/assert.hpp>
 #include <boost/mpl/integral_c.hpp>

test/strategies/andoyer.cpp

@@ -5,9 +5,10 @@
 // Copyright (c) 2008-2016 Bruno Lalande, Paris, France.
 // Copyright (c) 2009-2016 Mateusz Loskot, London, UK.
 
-// This file was modified by Oracle on 2014, 2015, 2016.
+// This file was modified by Oracle on 2014-2017.
 // Modifications copyright (c) 2014-2016 Oracle and/or its affiliates.
 
+// Contributed and/or modified by Vissarion Fysikopoulos, on behalf of Oracle
 // Contributed and/or modified by Adam Wulkiewicz, on behalf of Oracle
 
 // Parts of Boost.Geometry are redesigned from Geodan's Geographic Library
@@ -94,7 +95,24 @@ void test_distance(double lon1, double lat1, double lon2, double lat2, double ex
 
     return_type d_strategy = andoyer.apply(p1, p2);
     return_type d_function = bg::distance(p1, p2, andoyer);
-    return_type d_formula = andoyer_inverse_type::apply(to_rad(lon1), to_rad(lat1), to_rad(lon2), to_rad(lat2), stype()).distance;
+
+    double diff = bg::math::longitude_distance_signed<bg::degree>(lon1, lon2);
+    return_type d_formula;
+
+    // if the points lay on a meridian, distance strategy calls the special formula
+    // for meridian distance that returns different result than andoyer formula
+    // for nearly antipodal points
+    if (bg::math::equals(diff, 0.0)
+       || bg::math::equals(bg::math::abs(diff), 180.0))
+    {
+        d_formula = d_strategy;
+    }
+    else
+    {
+        d_formula = andoyer_inverse_type::apply(to_rad(lon1), to_rad(lat1),
+                                                to_rad(lon2), to_rad(lat2),
+                                                stype()).distance;
+    }
 
     BOOST_CHECK_CLOSE(d_strategy / 1000.0, expected_km, 0.001);
     BOOST_CHECK_CLOSE(d_function / 1000.0, expected_km, 0.001);
@@ -214,10 +232,10 @@ void test_all()
     // antipodal
     // ok? in those cases shorter path would pass through a pole
     // but 90 or -90 would be consistent with distance?
-    test_distazi<P1, P2>(0, 0,  180, 0, 20037.5, 0.0);
-    test_distazi<P1, P2>(0, 0, -180, 0, 20037.5, 0.0);
-    test_distazi<P1, P2>(-90, 0, 90, 0, 20037.5, 0.0);
-    test_distazi<P1, P2>(90, 0, -90, 0, 20037.5, 0.0);
+    test_distazi<P1, P2>(0, 0,  180, 0, 20003.9, 0.0);
+    test_distazi<P1, P2>(0, 0, -180, 0, 20003.9, 0.0);
+    test_distazi<P1, P2>(-90, 0, 90, 0, 20003.9, 0.0);
+    test_distazi<P1, P2>(90, 0, -90, 0, 20003.9, 0.0);
 
     // 0, 45, 90 ...
     for (int i = 0 ; i < 360 ; i += 45)
@@ -264,7 +282,7 @@ void test_all()
         test_distazi_symm<P1, P2>(normlized_deg(l-44.99), -44.99, normlized_deg(l+135), 45, 20008.1, 0.0);
         test_distazi_symm<P1, P2>(normlized_deg(l-44.999), -44.999, normlized_deg(l+135), 45, 20009.4, 0.0);
         // antipodal
-        test_distazi_symm<P1, P2>(normlized_deg(l-45), -45, normlized_deg(l+135), 45, 20020.7, 0.0);
+        test_distazi_symm<P1, P2>(normlized_deg(l-45), -45, normlized_deg(l+135), 45, 20003.92, 0.0);
     }
 
     /* SQL Server gives: