geometry/formulas/vertex_latitude.hpp

0001 // Boost.Geometry
0002
0003 // Copyright (c) 2016-2020 Oracle and/or its affiliates.
0004
0005 // Contributed and/or modified by Vissarion Fysikopoulos, on behalf of Oracle
0006 // Contributed and/or modified by Adam Wulkiewicz, on behalf of Oracle
0007
0008 // Use, modification and distribution is subject to the Boost Software License,
0009 // Version 1.0. (See accompanying file LICENSE_1_0.txt or copy at
0010 // http://www.boost.org/LICENSE_1_0.txt)
0011
0012 #ifndef BOOST_GEOMETRY_FORMULAS_MAXIMUM_LATITUDE_HPP
0013 #define BOOST_GEOMETRY_FORMULAS_MAXIMUM_LATITUDE_HPP
0014
0015
0016 #include <boost/geometry/core/static_assert.hpp>
0017 #include <boost/geometry/formulas/flattening.hpp>
0018 #include <boost/geometry/formulas/spherical.hpp>
0019
0020
0021 namespace boost { namespace geometry { namespace formula
0022 {
0023
0024 /*!
0025 \brief Algorithm to compute the vertex latitude of a geodesic segment. Vertex is
0026 a point on the geodesic that maximizes (or minimizes) the latitude.
0027 \author See
0028     [Wood96] Wood - Vertex Latitudes on Ellipsoid Geodesics, SIAM Rev., 38(4),
0029              637–644, 1996
0030 */
0031
0032 template <typename CT>
0033 class vertex_latitude_on_sphere
0034 {
0035
0036 public:
0037     template<typename T1, typename T2>
0038     static inline CT apply(T1 const& lat1,
0039                            T2 const& alp1)
0040     {
0041         return std::acos( math::abs(cos(lat1) * sin(alp1)) );
0042     }
0043 };
0044
0045 template <typename CT>
0046 class vertex_latitude_on_spheroid
0047 {
0048
0049 public:
0050 /*
0051  * formula based on paper
0052  *   [Wood96] Wood - Vertex Latitudes on Ellipsoid Geodesics, SIAM Rev., 38(4),
0053  *            637–644, 1996
0054     template <typename T1, typename T2, typename Spheroid>
0055     static inline CT apply(T1 const& lat1,
0056                            T2 const& alp1,
0057                            Spheroid const& spheroid)
0058     {
0059         CT const f = formula::flattening<CT>(spheroid);
0060
0061         CT const e2 = f * (CT(2) - f);
0062         CT const sin_alp1 = sin(alp1);
0063         CT const sin2_lat1 = math::sqr(sin(lat1));
0064         CT const cos2_lat1 = CT(1) - sin2_lat1;
0065
0066         CT const e2_sin2 = CT(1) - e2 * sin2_lat1;
0067         CT const cos2_sin2 = cos2_lat1 * math::sqr(sin_alp1);
0068         CT const vertex_lat = std::asin( math::sqrt((e2_sin2 - cos2_sin2)
0069                                                     / (e2_sin2 - e2 * cos2_sin2)));
0070         return vertex_lat;
0071     }
0072 */
0073
0074     // simpler formula based on Clairaut relation for spheroids
0075     template <typename T1, typename T2, typename Spheroid>
0076     static inline CT apply(T1 const& lat1,
0077                            T2 const& alp1,
0078                            Spheroid const& spheroid)
0079     {
0080         CT const f = formula::flattening<CT>(spheroid);
0081
0082         CT const one_minus_f = (CT(1) - f);
0083
0084         //get the reduced latitude
0085         CT const bet1 = atan( one_minus_f * tan(lat1) );
0086
0087         //apply Clairaut relation
0088         CT const betv =  vertex_latitude_on_sphere<CT>::apply(bet1, alp1);
0089
0090         //return the spheroid latitude
0091         return atan( tan(betv) / one_minus_f );
0092     }
0093
0094     /*
0095     template <typename T>
0096     inline static void sign_adjustment(CT lat1, CT lat2, CT vertex_lat, T& vrt_result)
0097     {
0098         // signbit returns a non-zero value (true) if the sign is negative;
0099         // and zero (false) otherwise.
0100         bool sign = std::signbit(std::abs(lat1) > std::abs(lat2) ? lat1 : lat2);
0101
0102         vrt_result.north = sign ? std::max(lat1, lat2) : vertex_lat;
0103         vrt_result.south = sign ? vertex_lat * CT(-1) : std::min(lat1, lat2);
0104     }
0105
0106     template <typename T>
0107     inline static bool vertex_on_segment(CT alp1, CT alp2, CT lat1, CT lat2, T& vrt_result)
0108     {
0109         CT const half_pi = math::pi<CT>() / CT(2);
0110
0111         // if the segment does not contain the vertex of the geodesic
0112         // then return the endpoint of max (min) latitude
0113         if ((alp1 < half_pi && alp2 < half_pi)
0114                 || (alp1 > half_pi && alp2 > half_pi))
0115         {
0116             vrt_result.north = std::max(lat1, lat2);
0117             vrt_result.south = std::min(lat1, lat2);
0118             return false;
0119         }
0120         return true;
0121     }
0122     */
0123 };
0124
0125
0126 template <typename CT, typename CS_Tag>
0127 struct vertex_latitude
0128 {
0129     BOOST_GEOMETRY_STATIC_ASSERT_FALSE(
0130         "Not implemented for this coordinate system.",
0131         CT, CS_Tag);
0132 };
0133
0134 template <typename CT>
0135 struct vertex_latitude<CT, spherical_equatorial_tag>
0136         : vertex_latitude_on_sphere<CT>
0137 {};
0138
0139 template <typename CT>
0140 struct vertex_latitude<CT, geographic_tag>
0141         : vertex_latitude_on_spheroid<CT>
0142 {};
0143
0144
0145 }}} // namespace boost::geometry::formula
0146
0147 #endif // BOOST_GEOMETRY_FORMULAS_MAXIMUM_LATITUDE_HPP