[formulas] [strategies] Thomas first order direct formula

include/boost/geometry/formulas/thomas_direct.hpp

@@ -1,7 +1,8 @@
 // Boost.Geometry
 
-// Copyright (c) 2016-2017 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
 
 // Use, modification and distribution is subject to the Boost Software License,
@@ -30,13 +31,12 @@ namespace boost { namespace geometry { namespace formula
 
 /*!
 \brief The solution of the direct problem of geodesics on latlong coordinates,
-       Forsyth-Andoyer-Lambert type approximation with second order terms.
+       Forsyth-Andoyer-Lambert type approximation with first/second order terms.
 \author See
     - Technical Report: PAUL D. THOMAS, MATHEMATICAL MODELS FOR NAVIGATION SYSTEMS, 1965
       http://www.dtic.mil/docs/citations/AD0627893
     - Technical Report: PAUL D. THOMAS, SPHEROIDAL GEODESICS, REFERENCE SYSTEMS, AND LOCAL GEOMETRY, 1970
       http://www.dtic.mil/docs/citations/AD0703541
-
 */
 template <
     typename CT,
@@ -59,7 +59,8 @@ public:
                                     T const& la1,
                                     Dist const& distance,
                                     Azi const& azimuth12,
-                                    Spheroid const& spheroid)
+                                    Spheroid const& spheroid,
+                                    bool SecondOrder = true)
     {
         result_type result;
 
@@ -91,7 +92,7 @@ public:
         CT azi12_alt = azimuth12;
         CT lat1_alt = lat1;
         bool alter_result = vflip_if_south(lat1, azimuth12, lat1_alt, azi12_alt);
-        
+
         CT const theta1 = math::equals(lat1_alt, pi_half) ? lat1_alt :
                           math::equals(lat1_alt, -pi_half) ? lat1_alt :
                           atan(one_minus_f * tan(lat1_alt));
@@ -108,9 +109,18 @@ public:
         CT const N = cos_theta1 * cos_a12;
         CT const C1 = f * M; // lower-case c1 in the technical report
         CT const C2 = f * (c1 - math::sqr(M)) / c4; // lower-case c2 in the technical report
-        CT const D = (c1 - C2) * (c1 - C2 - C1 * M);
-        CT const P = C2 * (c1 + C1 * M / c2) / D;
-
+        CT D = 0;
+        CT P = 0;
+        if ( BOOST_GEOMETRY_CONDITION(SecondOrder) )
+        {
+            D = (c1 - C2) * (c1 - C2 - C1 * M);
+            P = C2 * (c1 + C1 * M / c2) / D;
+        }
+        else
+        {
+            D = c1 - c2 * C2 - C1 * M;
+            P = C2 / D;
+        }
         // special case for equator:
         // sin_theta0 = 0 <=> lat1 = 0 ^ |azimuth12| = pi/2
         // NOTE: in this case it doesn't matter what's the value of cos_sigma1 because
@@ -132,9 +142,14 @@ public:
 
         CT const W = c1 - c2 * P * cos_u;
         CT const V = cos_u * cos_d - sin_u * sin_d;
-        CT const X = math::sqr(C2) * sin_d * cos_d * (2 * math::sqr(V) - c1);
         CT const Y = c2 * P * V * W * sin_d;
-        CT const d_sigma = d + X - Y;
+        CT X = 0;
+        CT d_sigma = d - Y;
+        if ( BOOST_GEOMETRY_CONDITION(SecondOrder) )
+        {
+            X = math::sqr(C2) * sin_d * cos_d * (2 * math::sqr(V) - c1);
+            d_sigma += X;
+        }
         CT const sin_d_sigma = sin(d_sigma);
         CT const cos_d_sigma = cos(d_sigma);
 
@@ -151,11 +166,16 @@ public:
         if (BOOST_GEOMETRY_CONDITION(CalcCoordinates))
         {
             CT const S_sigma = c2 * sigma1 - d_sigma;
-            CT const cos_S_sigma = cos(S_sigma);
+            CT cos_S_sigma = 0;
+            CT H = C1 * d_sigma;
+            if ( BOOST_GEOMETRY_CONDITION(SecondOrder) )
+            {
+                cos_S_sigma = cos(S_sigma);
+                H = H * (c1 - C2) - C1 * C2 * sin_d_sigma * cos_S_sigma;
+            }
             CT const d_eta = atan2(sin_d_sigma * sin_a12, cos_theta1 * cos_d_sigma - sin_theta1 * sin_d_sigma * cos_a12);
-            CT const H = C1 * (c1 - C2) * d_sigma - C1 * C2 * sin_d_sigma * cos_S_sigma;
             CT const d_lambda = d_eta - H;
-            
+
             result.lon2 = lon1 + d_lambda;
 
             if (! math::equals(M, c0))

include/boost/geometry/formulas/thomas_first_order_direct.hpp

@@ -0,0 +1,72 @@
+// 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_THOMAS_FIRST_ORDER_DIRECT_HPP
+#define BOOST_GEOMETRY_FORMULAS_THOMAS_FIRST_ORDER_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/differential_quantities.hpp>
+#include <boost/geometry/formulas/flattening.hpp>
+#include <boost/geometry/formulas/result_direct.hpp>
+
+
+namespace boost { namespace geometry { namespace formula
+{
+
+
+/*!
+\brief The solution of the direct problem of geodesics on latlong coordinates,
+       Forsyth-Andoyer-Lambert type approximation with first order terms.
+\author See
+    - Technical Report: PAUL D. THOMAS, MATHEMATICAL MODELS FOR NAVIGATION SYSTEMS, 1965
+      http://www.dtic.mil/docs/citations/AD0627893
+    - Technical Report: PAUL D. THOMAS, SPHEROIDAL GEODESICS, REFERENCE SYSTEMS, AND LOCAL GEOMETRY, 1970
+      http://www.dtic.mil/docs/citations/AD0703541
+*/
+template <
+    typename CT,
+    bool EnableCoordinates = true,
+    bool EnableReverseAzimuth = false,
+    bool EnableReducedLength = false,
+    bool EnableGeodesicScale = false
+>
+class thomas_first_order_direct
+{
+
+public:
+    typedef result_direct<CT> result_type;
+
+    template <typename T, typename Dist, typename Azi, typename Spheroid>
+    static inline result_type apply(T const& lo1,
+                                    T const& la1,
+                                    Dist const& distance,
+                                    Azi const& azimuth12,
+                                    Spheroid const& spheroid)
+    {   return thomas_direct<
+                CT,
+                EnableCoordinates,
+                EnableReverseAzimuth,
+                EnableReducedLength,
+                EnableGeodesicScale
+        >::apply(lo1, la1, distance, azimuth12, spheroid, false);
+    }
+};
+
+}}} // namespace boost::geometry::formula
+
+
+#endif // BOOST_GEOMETRY_FORMULAS_THOMAS_FIRST_ORDER_DIRECT_HPP

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

@@ -10,7 +10,7 @@
 #ifndef BOOST_GEOMETRY_STRATEGIES_GEOGRAPHIC_PARAMETERS_HPP
 #define BOOST_GEOMETRY_STRATEGIES_GEOGRAPHIC_PARAMETERS_HPP
 
-
+#include <boost/geometry/formulas/thomas_first_order_direct.hpp>
 #include <boost/geometry/formulas/andoyer_inverse.hpp>
 #include <boost/geometry/formulas/thomas_direct.hpp>
 #include <boost/geometry/formulas/thomas_inverse.hpp>
@@ -36,7 +36,7 @@ struct andoyer
         bool EnableGeodesicScale = false
     >
     struct direct
-            : formula::thomas_direct
+            : formula::thomas_first_order_direct
               <
                   CT, EnableCoordinates, EnableReverseAzimuth,
                   EnableReducedLength, EnableGeodesicScale