1 // Boost.Geometry 2 3 // Copyright (c) 2016-2020 Oracle and/or its affiliates. 4 5 // Contributed and/or modified by Vissarion Fysikopoulos, on behalf of Oracle 6 // Contributed and/or modified by Adam Wulkiewicz, on behalf of Oracle 7 8 // Use, modification and distribution is subject to the Boost Software License, 9 // Version 1.0. (See accompanying file LICENSE_1_0.txt or copy at 10 // http://www.boost.org/LICENSE_1_0.txt) 11 12 #ifndef BOOST_GEOMETRY_FORMULAS_MAXIMUM_LATITUDE_HPP 13 #define BOOST_GEOMETRY_FORMULAS_MAXIMUM_LATITUDE_HPP 14 15 16 #include <boost/geometry/core/static_assert.hpp> 17 #include <boost/geometry/formulas/flattening.hpp> 18 #include <boost/geometry/formulas/spherical.hpp> 19 20 21 namespace boost { namespace geometry { namespace formula 22 { 23 24 /*! 25 \brief Algorithm to compute the vertex latitude of a geodesic segment. Vertex is 26 a point on the geodesic that maximizes (or minimizes) the latitude. 27 \author See 28 [Wood96] Wood - Vertex Latitudes on Ellipsoid Geodesics, SIAM Rev., 38(4), 29 637–644, 1996 30 */ 31 32 template <typename CT> 33 class vertex_latitude_on_sphere 34 { 35 36 public: 37 template<typename T1, typename T2> apply(T1 const & lat1,T2 const & alp1)38 static inline CT apply(T1 const& lat1, 39 T2 const& alp1) 40 { 41 return std::acos( math::abs(cos(lat1) * sin(alp1)) ); 42 } 43 }; 44 45 template <typename CT> 46 class vertex_latitude_on_spheroid 47 { 48 49 public: 50 /* 51 * formula based on paper 52 * [Wood96] Wood - Vertex Latitudes on Ellipsoid Geodesics, SIAM Rev., 38(4), 53 * 637–644, 1996 54 template <typename T1, typename T2, typename Spheroid> 55 static inline CT apply(T1 const& lat1, 56 T2 const& alp1, 57 Spheroid const& spheroid) 58 { 59 CT const f = formula::flattening<CT>(spheroid); 60 61 CT const e2 = f * (CT(2) - f); 62 CT const sin_alp1 = sin(alp1); 63 CT const sin2_lat1 = math::sqr(sin(lat1)); 64 CT const cos2_lat1 = CT(1) - sin2_lat1; 65 66 CT const e2_sin2 = CT(1) - e2 * sin2_lat1; 67 CT const cos2_sin2 = cos2_lat1 * math::sqr(sin_alp1); 68 CT const vertex_lat = std::asin( math::sqrt((e2_sin2 - cos2_sin2) 69 / (e2_sin2 - e2 * cos2_sin2))); 70 return vertex_lat; 71 } 72 */ 73 74 // simpler formula based on Clairaut relation for spheroids 75 template <typename T1, typename T2, typename Spheroid> apply(T1 const & lat1,T2 const & alp1,Spheroid const & spheroid)76 static inline CT apply(T1 const& lat1, 77 T2 const& alp1, 78 Spheroid const& spheroid) 79 { 80 CT const f = formula::flattening<CT>(spheroid); 81 82 CT const one_minus_f = (CT(1) - f); 83 84 //get the reduced latitude 85 CT const bet1 = atan( one_minus_f * tan(lat1) ); 86 87 //apply Clairaut relation 88 CT const betv = vertex_latitude_on_sphere<CT>::apply(bet1, alp1); 89 90 //return the spheroid latitude 91 return atan( tan(betv) / one_minus_f ); 92 } 93 94 /* 95 template <typename T> 96 inline static void sign_adjustment(CT lat1, CT lat2, CT vertex_lat, T& vrt_result) 97 { 98 // signbit returns a non-zero value (true) if the sign is negative; 99 // and zero (false) otherwise. 100 bool sign = std::signbit(std::abs(lat1) > std::abs(lat2) ? lat1 : lat2); 101 102 vrt_result.north = sign ? std::max(lat1, lat2) : vertex_lat; 103 vrt_result.south = sign ? vertex_lat * CT(-1) : std::min(lat1, lat2); 104 } 105 106 template <typename T> 107 inline static bool vertex_on_segment(CT alp1, CT alp2, CT lat1, CT lat2, T& vrt_result) 108 { 109 CT const half_pi = math::pi<CT>() / CT(2); 110 111 // if the segment does not contain the vertex of the geodesic 112 // then return the endpoint of max (min) latitude 113 if ((alp1 < half_pi && alp2 < half_pi) 114 || (alp1 > half_pi && alp2 > half_pi)) 115 { 116 vrt_result.north = std::max(lat1, lat2); 117 vrt_result.south = std::min(lat1, lat2); 118 return false; 119 } 120 return true; 121 } 122 */ 123 }; 124 125 126 template <typename CT, typename CS_Tag> 127 struct vertex_latitude 128 { 129 BOOST_GEOMETRY_STATIC_ASSERT_FALSE( 130 "Not implemented for this coordinate system.", 131 CT, CS_Tag); 132 }; 133 134 template <typename CT> 135 struct vertex_latitude<CT, spherical_equatorial_tag> 136 : vertex_latitude_on_sphere<CT> 137 {}; 138 139 template <typename CT> 140 struct vertex_latitude<CT, geographic_tag> 141 : vertex_latitude_on_spheroid<CT> 142 {}; 143 144 145 }}} // namespace boost::geometry::formula 146 147 #endif // BOOST_GEOMETRY_FORMULAS_MAXIMUM_LATITUDE_HPP 148