Add direct formulas

include/boost/geometry/formulas/elliptic_meridian_arc_direct.hpp

@@ -0,0 +1,87 @@
+// Boost.Geometry
+
+// Copyright (c) 2018 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_MERIDIAN_ARC_DIRECT_HPP
+#define BOOST_GEOMETRY_FORMULAS_ELLIPTIC_MERIDIAN_ARC_DIRECT_HPP
+
+#include <boost/math/constants/constants.hpp>
+
+#include <boost/geometry/core/radius.hpp>
+
+#include <boost/geometry/util/condition.hpp>
+#include <boost/geometry/util/math.hpp>
+
+#include <boost/geometry/formulas/flattening.hpp>
+#include <boost/geometry/formulas/quarter_meridian.hpp>
+
+namespace boost { namespace geometry { namespace formula
+{
+
+/*!
+\brief Compute the direct geodesic problem on a meridian
+*/
+
+template <typename CT, unsigned int Order = 1>
+class elliptic_meridian_arc_direct
+{
+
+public :
+
+    // https://en.wikipedia.org/wiki/Meridian_arc#The_inverse_meridian_problem_for_the_ellipsoid
+    // latitudes are assumed to be in radians and in [-pi/2,pi/2]
+    template <typename T, typename Spheroid>
+    static CT apply(T m, Spheroid const& spheroid)
+    {
+        CT const f = formula::flattening<CT>(spheroid);
+        CT n = f / (CT(2) - f);
+        CT mp = formula::quarter_meridian<CT>(spheroid);
+        CT mu = geometry::math::pi<CT>()/CT(2) * m / mp;
+
+        if (Order == 0)
+        {
+            return mu;
+        }
+
+        CT H2 = 1.5 * n;
+
+        if (Order == 1)
+        {
+            return mu + H2 * sin(2*mu);
+        }
+
+        CT n2 = n * n;
+        CT H4 = 1.3125 * n2;
+
+        if (Order == 2)
+        {
+            return mu + H2 * sin(2*mu) + H4 * sin(4*mu);
+        }
+
+        CT n3 = n2 * n;
+        H2 -= 0.84375 * n3;
+        CT H6 = 1.572916667 * n3;
+
+        if (Order == 3)
+        {
+            return mu + H2 * sin(2*mu) + H4 * sin(4*mu) + H6 * sin(6*mu);
+        }
+
+        CT n4 = n2 * n2;
+        H4 -= 1.71875 * n4;
+        CT H8 = 2.142578125 * n4;
+
+        // Order 4 or higher
+        return mu + H2 * sin(2*mu) + H4 * sin(4*mu) + H6 * sin(6*mu) + H8 * sin(8*mu);
+    }
+};
+
+}}} // namespace boost::geometry::formula
+
+#endif // BOOST_GEOMETRY_FORMULAS_ELLIPTIC_MERIDIAN_ARC_DIRECT_HPP

include/boost/geometry/formulas/elliptic_arc_length.hpp → include/boost/geometry/formulas/elliptic_meridian_arc_inverse.hpp

@@ -1,15 +1,15 @@
 // Boost.Geometry
 
-// Copyright (c) 2017 Oracle and/or its affiliates.
+// Copyright (c) 2017-2018 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
+#ifndef BOOST_GEOMETRY_FORMULAS_ELLIPTIC_MERIDIAN_ARC_INVERSE_HPP
+#define BOOST_GEOMETRY_FORMULAS_ELLIPTIC_MERIDIAN_ARC_INVERSE_HPP
 
 #include <boost/math/constants/constants.hpp>
 
@@ -29,7 +29,7 @@ namespace boost { namespace geometry { namespace formula
 */
 
 template <typename CT, unsigned int Order = 1>
-class elliptic_arc_length
+class elliptic_meridian_arc_inverse
 {
 
 public :
@@ -167,68 +167,9 @@ public :
                   + 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
+#endif // BOOST_GEOMETRY_FORMULAS_ELLIPTIC_MERIDIAN_ARC_INVERSE_HPP

include/boost/geometry/formulas/quarter_meridian.hpp

@@ -0,0 +1,102 @@
+// Boost.Geometry
+
+// Copyright (c) 2018 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_QUARTER_MERIDIAN_HPP
+#define BOOST_GEOMETRY_FORMULAS_QUARTER_MERIDIAN_HPP
+
+#include <boost/geometry/core/radius.hpp>
+#include <boost/geometry/core/tag.hpp>
+#include <boost/geometry/core/tags.hpp>
+
+#include <boost/geometry/algorithms/not_implemented.hpp>
+
+namespace boost { namespace geometry
+{
+
+#ifndef DOXYGEN_NO_DISPATCH
+namespace formula_dispatch
+{
+
+template <typename ResultType, typename Geometry, typename Tag = typename tag<Geometry>::type>
+struct quarter_meridian
+    : not_implemented<Tag>
+{};
+
+template <typename ResultType, typename Geometry>
+struct quarter_meridian<ResultType, Geometry, srs_spheroid_tag>
+{
+    //https://en.wikipedia.org/wiki/Meridian_arc#Generalized_series
+    static inline ResultType apply(Geometry const& geometry)
+    {
+        //order 8 expansion
+        ResultType const C[] =
+        {
+            1073741824,
+            268435456,
+            16777216,
+            4194304,
+            1638400,
+            802816,
+            451584,
+            278784,
+            184041
+        };
+
+        ResultType const c2 = 2;
+        ResultType const c4 = 4;
+        ResultType const f = formula::flattening<ResultType>(geometry);
+        ResultType const n = f / (c2 - f);
+        ResultType const ab4 = (get_radius<0>(geometry)
+                                + get_radius<2>(geometry)) / c4;
+        return geometry::math::pi<ResultType>() * ab4 *
+                 horner_evaluate(n*n, C, C+8) / C[0];
+    }
+
+private :
+    //TODO: move the following to a more general space to be used by other
+    //      classes as well
+    /*
+        Evaluate the polynomial in x using Horner's method.
+    */
+    template <typename NT, typename IteratorType>
+    static inline NT horner_evaluate(NT x,
+                                     IteratorType begin,
+                                     IteratorType end)
+    {
+        NT result(0);
+        IteratorType it = end;
+        do
+        {
+            result = result * x + *--it;
+        }
+        while (it != begin);
+        return result;
+    }
+};
+
+} // namespace formula_dispatch
+#endif // DOXYGEN_NO_DISPATCH
+
+#ifndef DOXYGEN_NO_DETAIL
+namespace formula
+{
+
+template <typename ResultType, typename Geometry>
+ResultType quarter_meridian(Geometry const& geometry)
+{
+    return formula_dispatch::quarter_meridian<ResultType, Geometry>::apply(geometry);
+}
+
+} // namespace formula
+#endif // DOXYGEN_NO_DETAIL
+
+}} // namespace boost::geometry
+
+#endif // BOOST_GEOMETRY_FORMULAS_QUARTER_MERIDIAN_HPP

include/boost/geometry/formulas/spherical.hpp

@@ -1,6 +1,6 @@
 // Boost.Geometry
 
-// Copyright (c) 2016, Oracle and/or its affiliates.
+// Copyright (c) 2016-2018, 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
 
@@ -24,6 +24,8 @@
 #include <boost/geometry/util/normalize_spheroidal_coordinates.hpp>
 #include <boost/geometry/util/select_coordinate_type.hpp>
 
+#include <boost/geometry/formulas/result_direct.hpp>
+
 namespace boost { namespace geometry {
 
 namespace formula {
@@ -217,6 +219,47 @@ inline int azimuth_side_value(T const& azi_a1_p, T const& azi_a1_a2)
         : 1; // left
 }
 
+template <typename CT, typename Sphere>
+inline result_direct<CT> spherical_direct(CT const& lon1,
+                                          CT const& lat1,
+                                          CT const& sig12,
+                                          CT const& alp1,
+                                          Sphere const& sphere)
+{
+    result_direct<CT> result;
+
+    CT const sin_alp1 = sin(alp1);
+    CT const sin_lat1 = sin(lat1);
+    CT const cos_alp1 = cos(alp1);
+    CT const cos_lat1 = cos(lat1);
+
+    CT const norm = math::sqrt(cos_alp1 * cos_alp1 + sin_alp1 * sin_alp1
+                                                   * sin_lat1 * sin_lat1);
+    CT const alp0 = atan2(sin_alp1 * cos_lat1, norm);
+    CT const sig1 = atan2(sin_lat1, cos_alp1 * cos_lat1);
+    CT const sig2 = sig1 + sig12 / get_radius<0>(sphere);
+
+    CT const sin_sig2 = sin(sig2);
+    CT const cos_sig2 = cos(sig2);
+    CT const sin_sig1 = sin(sig1);
+    CT const cos_sig1 = cos(sig1);
+    CT const sin_alp0 = sin(alp0);
+    CT const cos_alp0 = cos(alp0);
+
+    CT const norm2 = math::sqrt(cos_alp0 * cos_alp0 * cos_sig2 * cos_sig2
+                                + sin_alp0 * sin_alp0);
+    CT const lat2 = atan2(cos_alp0 * sin_sig2, norm2);
+    CT const alp2 = atan2(sin_alp0, cos_alp0 * cos_sig2);
+    CT const omg1 = atan2(sin_alp0 * sin_sig1, cos_sig1);
+    CT const lon2 = atan2(sin_alp0 * sin_sig2, cos_sig2);
+
+    result.lon2 = lon1 + lon2 - omg1;
+    result.lat2 = lat2;
+    result.reverse_azimuth = alp2;
+
+    return result;
+}
+
 } // namespace formula
 
 }} // namespace boost::geometry