diff options
Diffstat (limited to 'boost/geometry/formulas/meridian_direct.hpp')
| -rw-r--r-- | boost/geometry/formulas/meridian_direct.hpp | 168 |
1 files changed, 168 insertions, 0 deletions
diff --git a/boost/geometry/formulas/meridian_direct.hpp b/boost/geometry/formulas/meridian_direct.hpp new file mode 100644 index 0000000000..e55b35e88f --- /dev/null +++ b/boost/geometry/formulas/meridian_direct.hpp @@ -0,0 +1,168 @@ +// 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_MERIDIAN_DIRECT_HPP +#define BOOST_GEOMETRY_FORMULAS_MERIDIAN_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/meridian_inverse.hpp> +#include <boost/geometry/formulas/flattening.hpp> +#include <boost/geometry/formulas/quarter_meridian.hpp> +#include <boost/geometry/formulas/result_direct.hpp> + +namespace boost { namespace geometry { namespace formula +{ + +/*! +\brief Compute the direct geodesic problem on a meridian +*/ + +template < + typename CT, + bool EnableCoordinates = true, + bool EnableReverseAzimuth = false, + bool EnableReducedLength = false, + bool EnableGeodesicScale = false, + unsigned int Order = 4 +> +class meridian_direct +{ + static const bool CalcQuantities = EnableReducedLength || EnableGeodesicScale; + static const bool CalcRevAzimuth = EnableReverseAzimuth || CalcQuantities; + static const bool CalcCoordinates = EnableCoordinates || CalcRevAzimuth; + +public: + typedef result_direct<CT> result_type; + + template <typename T, typename Dist, typename Spheroid> + static inline result_type apply(T const& lo1, + T const& la1, + Dist const& distance, + bool north, + Spheroid const& spheroid) + { + result_type result; + + CT const half_pi = math::half_pi<CT>(); + CT const pi = math::pi<CT>(); + CT const one_and_a_half_pi = pi + half_pi; + CT const c0 = 0; + + CT azimuth = north ? c0 : pi; + + if (BOOST_GEOMETRY_CONDITION(CalcCoordinates)) + { + CT s0 = meridian_inverse<CT, Order>::apply(la1, spheroid); + int signed_distance = north ? distance : -distance; + result.lon2 = lo1; + result.lat2 = apply(s0 + signed_distance, spheroid); + } + + if (BOOST_GEOMETRY_CONDITION(CalcRevAzimuth)) + { + result.reverse_azimuth = azimuth; + + + if (result.lat2 > half_pi && + result.lat2 < one_and_a_half_pi) + { + result.reverse_azimuth = pi; + } + else if (result.lat2 < -half_pi && + result.lat2 > -one_and_a_half_pi) + { + result.reverse_azimuth = c0; + } + + } + + if (BOOST_GEOMETRY_CONDITION(CalcQuantities)) + { + CT const b = CT(get_radius<2>(spheroid)); + CT const f = formula::flattening<CT>(spheroid); + + boost::geometry::math::normalize_spheroidal_coordinates + < + boost::geometry::radian, + double + >(result.lon2, result.lat2); + + typedef differential_quantities + < + CT, + EnableReducedLength, + EnableGeodesicScale, + Order + > quantities; + quantities::apply(lo1, la1, result.lon2, result.lat2, + azimuth, result.reverse_azimuth, + b, f, + result.reduced_length, result.geodesic_scale); + } + return result; + } + + // 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_MERIDIAN_DIRECT_HPP |