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