diff --git a/include/boost/geometry/formulas/thomas_direct.hpp b/include/boost/geometry/formulas/thomas_direct.hpp index 6a7ac3e41..acf363d75 100644 --- a/include/boost/geometry/formulas/thomas_direct.hpp +++ b/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)) diff --git a/include/boost/geometry/formulas/thomas_first_order_direct.hpp b/include/boost/geometry/formulas/thomas_first_order_direct.hpp new file mode 100644 index 000000000..1fa356787 --- /dev/null +++ b/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 + +#include + +#include +#include + +#include +#include +#include + + +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 result_type; + + template + 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 diff --git a/include/boost/geometry/strategies/geographic/parameters.hpp b/include/boost/geometry/strategies/geographic/parameters.hpp index 4896e42e8..8ac427526 100644 --- a/include/boost/geometry/strategies/geographic/parameters.hpp +++ b/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 #include #include #include @@ -36,7 +36,7 @@ struct andoyer bool EnableGeodesicScale = false > struct direct - : formula::thomas_direct + : formula::thomas_first_order_direct < CT, EnableCoordinates, EnableReverseAzimuth, EnableReducedLength, EnableGeodesicScale