diff options
Diffstat (limited to 'boost/geometry/formulas')
| -rw-r--r-- | boost/geometry/formulas/vertex_longitude.hpp | 339 |
1 files changed, 339 insertions, 0 deletions
diff --git a/boost/geometry/formulas/vertex_longitude.hpp b/boost/geometry/formulas/vertex_longitude.hpp new file mode 100644 index 0000000000..00f2fd4e7a --- /dev/null +++ b/boost/geometry/formulas/vertex_longitude.hpp @@ -0,0 +1,339 @@ +// Boost.Geometry + +// Copyright (c) 2016-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, +// 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_MAXIMUM_LONGITUDE_HPP +#define BOOST_GEOMETRY_FORMULAS_MAXIMUM_LONGITUDE_HPP + +#include <boost/geometry/formulas/spherical.hpp> +#include <boost/geometry/formulas/flattening.hpp> +#include <boost/geometry/core/srs.hpp> +#include <boost/mpl/assert.hpp> + +#include <boost/math/special_functions/hypot.hpp> + +namespace boost { namespace geometry { namespace formula +{ + +/*! +\brief Algorithm to compute the vertex longitude of a geodesic segment. Vertex is +a point on the geodesic that maximizes (or minimizes) the latitude. The algorithm +is given the vertex latitude. +*/ + +//Classes for spesific CS + +template <typename CT> +class vertex_longitude_on_sphere +{ + +public: + + template <typename T> + static inline CT apply(T const& lat1, //segment point 1 + T const& lat2, //segment point 2 + T const& lat3, //vertex latitude + T const& sin_l12, + T const& cos_l12) //lon1 -lon2 + { + //https://en.wikipedia.org/wiki/Great-circle_navigation#Finding_way-points + CT const A = sin(lat1) * cos(lat2) * cos(lat3) * sin_l12; + CT const B = sin(lat1) * cos(lat2) * cos(lat3) * cos_l12 + - cos(lat1) * sin(lat2) * cos(lat3); + CT lon = atan2(B, A); + return lon + math::pi<CT>(); + } +}; + +template <typename CT> +class vertex_longitude_on_spheroid +{ + template<typename T> + static inline void normalize(T& x, T& y) + { + T h = boost::math::hypot(x, y); + x /= h; + y /= h; + } + +public: + + template <typename T, typename Spheroid> + static inline CT apply(T const& lat1, //segment point 1 + T const& lat2, //segment point 2 + T const& lat3, //vertex latitude + T& alp1, + Spheroid const& spheroid) + { + // We assume that segment points lay on different side w.r.t. + // the vertex + + // Constants + CT const c0 = 0; + CT const c2 = 2; + CT const half_pi = math::pi<CT>() / c2; + if (math::equals(lat1, half_pi) + || math::equals(lat2, half_pi) + || math::equals(lat1, -half_pi) + || math::equals(lat2, -half_pi)) + { + // one segment point is the pole + return c0; + } + + // More constants + CT const f = flattening<CT>(spheroid); + CT const pi = math::pi<CT>(); + CT const c1 = 1; + CT const cminus1 = -1; + + // First, compute longitude on auxiliary sphere + + CT const one_minus_f = c1 - f; + CT const bet1 = atan(one_minus_f * tan(lat1)); + CT const bet2 = atan(one_minus_f * tan(lat2)); + CT const bet3 = atan(one_minus_f * tan(lat3)); + + CT cos_bet1 = cos(bet1); + CT cos_bet2 = cos(bet2); + CT const sin_bet1 = sin(bet1); + CT const sin_bet2 = sin(bet2); + CT const sin_bet3 = sin(bet3); + + CT omg12 = 0; + + if (bet1 < c0) + { + cos_bet1 *= cminus1; + omg12 += pi; + } + if (bet2 < c0) + { + cos_bet2 *= cminus1; + omg12 += pi; + } + + CT const sin_alp1 = sin(alp1); + CT const cos_alp1 = math::sqrt(c1 - math::sqr(sin_alp1)); + + CT const norm = math::sqrt(math::sqr(cos_alp1) + math::sqr(sin_alp1 * sin_bet1)); + CT const sin_alp0 = sin(atan2(sin_alp1 * cos_bet1, norm)); + + BOOST_ASSERT(cos_bet2 != c0); + CT const sin_alp2 = sin_alp1 * cos_bet1 / cos_bet2; + + CT const cos_alp0 = math::sqrt(c1 - math::sqr(sin_alp0)); + CT const cos_alp2 = math::sqrt(c1 - math::sqr(sin_alp2)); + + CT const sig1 = atan2(sin_bet1, cos_alp1 * cos_bet1); + CT const sig2 = atan2(sin_bet2, -cos_alp2 * cos_bet2); //lat3 is a vertex + + CT const cos_sig1 = cos(sig1); + CT const sin_sig1 = math::sqrt(c1 - math::sqr(cos_sig1)); + + CT const cos_sig2 = cos(sig2); + CT const sin_sig2 = math::sqrt(c1 - math::sqr(cos_sig2)); + + CT const omg1 = atan2(sin_alp0 * sin_sig1, cos_sig1); + CT const omg2 = atan2(sin_alp0 * sin_sig2, cos_sig2); + + omg12 += omg1 - omg2; + + CT const sin_omg12 = sin(omg12); + CT const cos_omg12 = cos(omg12); + + CT omg13 = geometry::formula::vertex_longitude_on_sphere<CT> + ::apply(bet1, bet2, bet3, sin_omg12, cos_omg12); + + if (lat1 * lat2 < c0)//different hemispheres + { + if ((lat2 - lat1) * lat3 > c0)// ascending segment + { + omg13 = pi - omg13; + } + } + + // Second, compute the ellipsoidal longitude + + CT const e2 = f * (c2 - f); + CT const ep = math::sqrt(e2 / (c1 - e2)); + CT const k2 = math::sqr(ep * cos_alp0); + CT const sqrt_k2_plus_one = math::sqrt(c1 + k2); + CT const eps = (sqrt_k2_plus_one - c1) / (sqrt_k2_plus_one + c1); + CT const eps2 = eps * eps; + CT const n = f / (c2 - f); + + // sig3 is the length from equator to the vertex + CT sig3; + if(sin_bet3 > c0) + { + sig3 = half_pi; + } else { + sig3 = -half_pi; + } + CT const cos_sig3 = 0; + CT const sin_sig3 = 1; + + CT sig13 = sig3 - sig1; + if (sig13 > pi) + { + sig13 -= 2 * pi; + } + + // Order 2 approximation + CT const c1over2 = 0.5; + CT const c1over4 = 0.25; + CT const c1over8 = 0.125; + CT const c1over16 = 0.0625; + CT const c4 = 4; + CT const c8 = 8; + + CT const A3 = 1 - (c1over2 - c1over2 * n) * eps - c1over4 * eps2; + CT const C31 = (c1over4 - c1over4 * n) * eps + c1over8 * eps2; + CT const C32 = c1over16 * eps2; + + CT const sin2_sig3 = c2 * cos_sig3 * sin_sig3; + CT const sin4_sig3 = sin_sig3 * (-c4 * cos_sig3 + + c8 * cos_sig3 * cos_sig3 * cos_sig3); + CT const sin2_sig1 = c2 * cos_sig1 * sin_sig1; + CT const sin4_sig1 = sin_sig1 * (-c4 * cos_sig1 + + c8 * cos_sig1 * cos_sig1 * cos_sig1); + CT const I3 = A3 * (sig13 + + C31 * (sin2_sig3 - sin2_sig1) + + C32 * (sin4_sig3 - sin4_sig1)); + + int sign = c1; + if (bet3 < c0) + { + sign = cminus1; + } + + CT const dlon_max = omg13 - sign * f * sin_alp0 * I3; + + return dlon_max; + } +}; + +//CS_tag dispatching + +template <typename CT, typename CS_Tag> +struct compute_vertex_lon +{ + BOOST_MPL_ASSERT_MSG + ( + false, NOT_IMPLEMENTED_FOR_THIS_COORDINATE_SYSTEM, (types<CS_Tag>) + ); + +}; + +template <typename CT> +struct compute_vertex_lon<CT, spherical_equatorial_tag> +{ + template <typename Strategy> + static inline CT apply(CT const& lat1, + CT const& lat2, + CT const& vertex_lat, + CT const& sin_l12, + CT const& cos_l12, + CT, + Strategy) + { + return vertex_longitude_on_sphere<CT> + ::apply(lat1, + lat2, + vertex_lat, + sin_l12, + cos_l12); + } +}; + +template <typename CT> +struct compute_vertex_lon<CT, geographic_tag> +{ + template <typename Strategy> + static inline CT apply(CT const& lat1, + CT const& lat2, + CT const& vertex_lat, + CT, + CT, + CT& alp1, + Strategy const& azimuth_strategy) + { + return vertex_longitude_on_spheroid<CT> + ::apply(lat1, + lat2, + vertex_lat, + alp1, + azimuth_strategy.model()); + } +}; + +// Vertex longitude interface +// Assume that lon1 < lon2 and vertex_lat is the latitude of the vertex + +template <typename CT, typename CS_Tag> +class vertex_longitude +{ +public : + template <typename Strategy> + static inline CT apply(CT& lon1, + CT& lat1, + CT& lon2, + CT& lat2, + CT const& vertex_lat, + CT& alp1, + Strategy const& azimuth_strategy) + { + CT const c0 = 0; + CT pi = math::pi<CT>(); + + //Vertex is a segment's point + if (math::equals(vertex_lat, lat1)) + { + return lon1; + } + if (math::equals(vertex_lat, lat2)) + { + return lon2; + } + + //Segment lay on meridian + if (math::equals(lon1, lon2)) + { + return (std::max)(lat1, lat2); + } + BOOST_ASSERT(lon1 < lon2); + + CT dlon = compute_vertex_lon<CT, CS_Tag>::apply(lat1, lat2, + vertex_lat, + sin(lon1 - lon2), + cos(lon1 - lon2), + alp1, + azimuth_strategy); + + CT vertex_lon = std::fmod(lon1 + dlon, 2 * pi); + + if (vertex_lat < c0) + { + vertex_lon -= pi; + } + + if (std::abs(lon1 - lon2) > pi) + { + vertex_lon -= pi; + } + + return vertex_lon; + } +}; + +}}} // namespace boost::geometry::formula +#endif // BOOST_GEOMETRY_FORMULAS_MAXIMUM_LONGITUDE_HPP + |