diff options
Diffstat (limited to 'boost/geometry/formulas/vertex_latitude.hpp')
-rw-r--r-- | boost/geometry/formulas/vertex_latitude.hpp | 148 |
1 files changed, 148 insertions, 0 deletions
diff --git a/boost/geometry/formulas/vertex_latitude.hpp b/boost/geometry/formulas/vertex_latitude.hpp new file mode 100644 index 0000000000..755067b08d --- /dev/null +++ b/boost/geometry/formulas/vertex_latitude.hpp @@ -0,0 +1,148 @@ +// 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_LATITUDE_HPP +#define BOOST_GEOMETRY_FORMULAS_MAXIMUM_LATITUDE_HPP + +#include <boost/geometry/core/srs.hpp> +#include <boost/geometry/formulas/flattening.hpp> +#include <boost/geometry/formulas/spherical.hpp> +#include <boost/mpl/assert.hpp> + +namespace boost { namespace geometry { namespace formula +{ + +/*! +\brief Algorithm to compute the vertex latitude of a geodesic segment. Vertex is +a point on the geodesic that maximizes (or minimizes) the latitude. +\author See + [Wood96] Wood - Vertex Latitudes on Ellipsoid Geodesics, SIAM Rev., 38(4), + 637–644, 1996 +*/ + +template <typename CT> +class vertex_latitude_on_sphere +{ + +public: + template<typename T1, typename T2> + static inline CT apply(T1 const& lat1, + T2 const& alp1) + { + return std::acos( math::abs(cos(lat1) * sin(alp1)) ); + } +}; + +template <typename CT> +class vertex_latitude_on_spheroid +{ + +public: +/* + * formula based on paper + * [Wood96] Wood - Vertex Latitudes on Ellipsoid Geodesics, SIAM Rev., 38(4), + * 637–644, 1996 + template <typename T1, typename T2, typename Spheroid> + static inline CT apply(T1 const& lat1, + T2 const& alp1, + Spheroid const& spheroid) + { + CT const f = formula::flattening<CT>(spheroid); + + CT const e2 = f * (CT(2) - f); + CT const sin_alp1 = sin(alp1); + CT const sin2_lat1 = math::sqr(sin(lat1)); + CT const cos2_lat1 = CT(1) - sin2_lat1; + + CT const e2_sin2 = CT(1) - e2 * sin2_lat1; + CT const cos2_sin2 = cos2_lat1 * math::sqr(sin_alp1); + CT const vertex_lat = std::asin( math::sqrt((e2_sin2 - cos2_sin2) + / (e2_sin2 - e2 * cos2_sin2))); + return vertex_lat; + } +*/ + + // simpler formula based on Clairaut relation for spheroids + template <typename T1, typename T2, typename Spheroid> + static inline CT apply(T1 const& lat1, + T2 const& alp1, + Spheroid const& spheroid) + { + CT const f = formula::flattening<CT>(spheroid); + + CT const one_minus_f = (CT(1) - f); + + //get the reduced latitude + CT const bet1 = atan( one_minus_f * tan(lat1) ); + + //apply Clairaut relation + CT const betv = vertex_latitude_on_sphere<CT>::apply(bet1, alp1); + + //return the spheroid latitude + return atan( tan(betv) / one_minus_f ); + } + + /* + template <typename T> + inline static void sign_adjustment(CT lat1, CT lat2, CT vertex_lat, T& vrt_result) + { + // signbit returns a non-zero value (true) if the sign is negative; + // and zero (false) otherwise. + bool sign = std::signbit(std::abs(lat1) > std::abs(lat2) ? lat1 : lat2); + + vrt_result.north = sign ? std::max(lat1, lat2) : vertex_lat; + vrt_result.south = sign ? vertex_lat * CT(-1) : std::min(lat1, lat2); + } + + template <typename T> + inline static bool vertex_on_segment(CT alp1, CT alp2, CT lat1, CT lat2, T& vrt_result) + { + CT const half_pi = math::pi<CT>() / CT(2); + + // if the segment does not contain the vertex of the geodesic + // then return the endpoint of max (min) latitude + if ((alp1 < half_pi && alp2 < half_pi) + || (alp1 > half_pi && alp2 > half_pi)) + { + vrt_result.north = std::max(lat1, lat2); + vrt_result.south = std::min(lat1, lat2); + return false; + } + return true; + } + */ +}; + + +template <typename CT, typename CS_Tag> +struct vertex_latitude +{ + BOOST_MPL_ASSERT_MSG + ( + false, NOT_IMPLEMENTED_FOR_THIS_COORDINATE_SYSTEM, (types<CS_Tag>) + ); + +}; + +template <typename CT> +struct vertex_latitude<CT, spherical_equatorial_tag> + : vertex_latitude_on_sphere<CT> +{}; + +template <typename CT> +struct vertex_latitude<CT, geographic_tag> + : vertex_latitude_on_spheroid<CT> +{}; + + +}}} // namespace boost::geometry::formula + +#endif // BOOST_GEOMETRY_FORMULAS_MAXIMUM_LATITUDE_HPP |