1 // Boost.Geometry 2 3 // Copyright (c) 2016-2018 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_LONGITUDE_HPP 13 #define BOOST_GEOMETRY_FORMULAS_MAXIMUM_LONGITUDE_HPP 14 15 #include <boost/geometry/formulas/spherical.hpp> 16 #include <boost/geometry/formulas/flattening.hpp> 17 18 #include <boost/mpl/assert.hpp> 19 20 #include <boost/math/special_functions/hypot.hpp> 21 22 namespace boost { namespace geometry { namespace formula 23 { 24 25 /*! 26 \brief Algorithm to compute the vertex longitude of a geodesic segment. Vertex is 27 a point on the geodesic that maximizes (or minimizes) the latitude. The algorithm 28 is given the vertex latitude. 29 */ 30 31 //Classes for spesific CS 32 33 template <typename CT> 34 class vertex_longitude_on_sphere 35 { 36 37 public: 38 39 template <typename T> apply(T const & lat1,T const & lat2,T const & lat3,T const & sin_l12,T const & cos_l12)40 static inline CT apply(T const& lat1, //segment point 1 41 T const& lat2, //segment point 2 42 T const& lat3, //vertex latitude 43 T const& sin_l12, 44 T const& cos_l12) //lon1 -lon2 45 { 46 //https://en.wikipedia.org/wiki/Great-circle_navigation#Finding_way-points 47 CT const A = sin(lat1) * cos(lat2) * cos(lat3) * sin_l12; 48 CT const B = sin(lat1) * cos(lat2) * cos(lat3) * cos_l12 49 - cos(lat1) * sin(lat2) * cos(lat3); 50 CT lon = atan2(B, A); 51 return lon + math::pi<CT>(); 52 } 53 }; 54 55 template <typename CT> 56 class vertex_longitude_on_spheroid 57 { 58 template<typename T> normalize(T & x,T & y)59 static inline void normalize(T& x, T& y) 60 { 61 T h = boost::math::hypot(x, y); 62 x /= h; 63 y /= h; 64 } 65 66 public: 67 68 template <typename T, typename Spheroid> apply(T const & lat1,T const & lat2,T const & lat3,T & alp1,Spheroid const & spheroid)69 static inline CT apply(T const& lat1, //segment point 1 70 T const& lat2, //segment point 2 71 T const& lat3, //vertex latitude 72 T& alp1, 73 Spheroid const& spheroid) 74 { 75 // We assume that segment points lay on different side w.r.t. 76 // the vertex 77 78 // Constants 79 CT const c0 = 0; 80 CT const c2 = 2; 81 CT const half_pi = math::pi<CT>() / c2; 82 if (math::equals(lat1, half_pi) 83 || math::equals(lat2, half_pi) 84 || math::equals(lat1, -half_pi) 85 || math::equals(lat2, -half_pi)) 86 { 87 // one segment point is the pole 88 return c0; 89 } 90 91 // More constants 92 CT const f = flattening<CT>(spheroid); 93 CT const pi = math::pi<CT>(); 94 CT const c1 = 1; 95 CT const cminus1 = -1; 96 97 // First, compute longitude on auxiliary sphere 98 99 CT const one_minus_f = c1 - f; 100 CT const bet1 = atan(one_minus_f * tan(lat1)); 101 CT const bet2 = atan(one_minus_f * tan(lat2)); 102 CT const bet3 = atan(one_minus_f * tan(lat3)); 103 104 CT cos_bet1 = cos(bet1); 105 CT cos_bet2 = cos(bet2); 106 CT const sin_bet1 = sin(bet1); 107 CT const sin_bet2 = sin(bet2); 108 CT const sin_bet3 = sin(bet3); 109 110 CT omg12 = 0; 111 112 if (bet1 < c0) 113 { 114 cos_bet1 *= cminus1; 115 omg12 += pi; 116 } 117 if (bet2 < c0) 118 { 119 cos_bet2 *= cminus1; 120 omg12 += pi; 121 } 122 123 CT const sin_alp1 = sin(alp1); 124 CT const cos_alp1 = math::sqrt(c1 - math::sqr(sin_alp1)); 125 126 CT const norm = math::sqrt(math::sqr(cos_alp1) + math::sqr(sin_alp1 * sin_bet1)); 127 CT const sin_alp0 = sin(atan2(sin_alp1 * cos_bet1, norm)); 128 129 BOOST_ASSERT(cos_bet2 != c0); 130 CT const sin_alp2 = sin_alp1 * cos_bet1 / cos_bet2; 131 132 CT const cos_alp0 = math::sqrt(c1 - math::sqr(sin_alp0)); 133 CT const cos_alp2 = math::sqrt(c1 - math::sqr(sin_alp2)); 134 135 CT const sig1 = atan2(sin_bet1, cos_alp1 * cos_bet1); 136 CT const sig2 = atan2(sin_bet2, -cos_alp2 * cos_bet2); //lat3 is a vertex 137 138 CT const cos_sig1 = cos(sig1); 139 CT const sin_sig1 = math::sqrt(c1 - math::sqr(cos_sig1)); 140 141 CT const cos_sig2 = cos(sig2); 142 CT const sin_sig2 = math::sqrt(c1 - math::sqr(cos_sig2)); 143 144 CT const omg1 = atan2(sin_alp0 * sin_sig1, cos_sig1); 145 CT const omg2 = atan2(sin_alp0 * sin_sig2, cos_sig2); 146 147 omg12 += omg1 - omg2; 148 149 CT const sin_omg12 = sin(omg12); 150 CT const cos_omg12 = cos(omg12); 151 152 CT omg13 = geometry::formula::vertex_longitude_on_sphere<CT> 153 ::apply(bet1, bet2, bet3, sin_omg12, cos_omg12); 154 155 if (lat1 * lat2 < c0)//different hemispheres 156 { 157 if ((lat2 - lat1) * lat3 > c0)// ascending segment 158 { 159 omg13 = pi - omg13; 160 } 161 } 162 163 // Second, compute the ellipsoidal longitude 164 165 CT const e2 = f * (c2 - f); 166 CT const ep = math::sqrt(e2 / (c1 - e2)); 167 CT const k2 = math::sqr(ep * cos_alp0); 168 CT const sqrt_k2_plus_one = math::sqrt(c1 + k2); 169 CT const eps = (sqrt_k2_plus_one - c1) / (sqrt_k2_plus_one + c1); 170 CT const eps2 = eps * eps; 171 CT const n = f / (c2 - f); 172 173 // sig3 is the length from equator to the vertex 174 CT sig3; 175 if(sin_bet3 > c0) 176 { 177 sig3 = half_pi; 178 } else { 179 sig3 = -half_pi; 180 } 181 CT const cos_sig3 = 0; 182 CT const sin_sig3 = 1; 183 184 CT sig13 = sig3 - sig1; 185 if (sig13 > pi) 186 { 187 sig13 -= 2 * pi; 188 } 189 190 // Order 2 approximation 191 CT const c1over2 = 0.5; 192 CT const c1over4 = 0.25; 193 CT const c1over8 = 0.125; 194 CT const c1over16 = 0.0625; 195 CT const c4 = 4; 196 CT const c8 = 8; 197 198 CT const A3 = 1 - (c1over2 - c1over2 * n) * eps - c1over4 * eps2; 199 CT const C31 = (c1over4 - c1over4 * n) * eps + c1over8 * eps2; 200 CT const C32 = c1over16 * eps2; 201 202 CT const sin2_sig3 = c2 * cos_sig3 * sin_sig3; 203 CT const sin4_sig3 = sin_sig3 * (-c4 * cos_sig3 204 + c8 * cos_sig3 * cos_sig3 * cos_sig3); 205 CT const sin2_sig1 = c2 * cos_sig1 * sin_sig1; 206 CT const sin4_sig1 = sin_sig1 * (-c4 * cos_sig1 207 + c8 * cos_sig1 * cos_sig1 * cos_sig1); 208 CT const I3 = A3 * (sig13 209 + C31 * (sin2_sig3 - sin2_sig1) 210 + C32 * (sin4_sig3 - sin4_sig1)); 211 212 CT const sign = bet3 >= c0 213 ? c1 214 : cminus1; 215 216 CT const dlon_max = omg13 - sign * f * sin_alp0 * I3; 217 218 return dlon_max; 219 } 220 }; 221 222 //CS_tag dispatching 223 224 template <typename CT, typename CS_Tag> 225 struct compute_vertex_lon 226 { 227 BOOST_MPL_ASSERT_MSG 228 ( 229 false, NOT_IMPLEMENTED_FOR_THIS_COORDINATE_SYSTEM, (types<CS_Tag>) 230 ); 231 232 }; 233 234 template <typename CT> 235 struct compute_vertex_lon<CT, spherical_equatorial_tag> 236 { 237 template <typename Strategy> applyboost::geometry::formula::compute_vertex_lon238 static inline CT apply(CT const& lat1, 239 CT const& lat2, 240 CT const& vertex_lat, 241 CT const& sin_l12, 242 CT const& cos_l12, 243 CT, 244 Strategy) 245 { 246 return vertex_longitude_on_sphere<CT> 247 ::apply(lat1, 248 lat2, 249 vertex_lat, 250 sin_l12, 251 cos_l12); 252 } 253 }; 254 255 template <typename CT> 256 struct compute_vertex_lon<CT, geographic_tag> 257 { 258 template <typename Strategy> applyboost::geometry::formula::compute_vertex_lon259 static inline CT apply(CT const& lat1, 260 CT const& lat2, 261 CT const& vertex_lat, 262 CT, 263 CT, 264 CT& alp1, 265 Strategy const& azimuth_strategy) 266 { 267 return vertex_longitude_on_spheroid<CT> 268 ::apply(lat1, 269 lat2, 270 vertex_lat, 271 alp1, 272 azimuth_strategy.model()); 273 } 274 }; 275 276 // Vertex longitude interface 277 // Assume that lon1 < lon2 and vertex_lat is the latitude of the vertex 278 279 template <typename CT, typename CS_Tag> 280 class vertex_longitude 281 { 282 public : 283 template <typename Strategy> apply(CT & lon1,CT & lat1,CT & lon2,CT & lat2,CT const & vertex_lat,CT & alp1,Strategy const & azimuth_strategy)284 static inline CT apply(CT& lon1, 285 CT& lat1, 286 CT& lon2, 287 CT& lat2, 288 CT const& vertex_lat, 289 CT& alp1, 290 Strategy const& azimuth_strategy) 291 { 292 CT const c0 = 0; 293 CT pi = math::pi<CT>(); 294 295 //Vertex is a segment's point 296 if (math::equals(vertex_lat, lat1)) 297 { 298 return lon1; 299 } 300 if (math::equals(vertex_lat, lat2)) 301 { 302 return lon2; 303 } 304 305 //Segment lay on meridian 306 if (math::equals(lon1, lon2)) 307 { 308 return (std::max)(lat1, lat2); 309 } 310 BOOST_ASSERT(lon1 < lon2); 311 312 CT dlon = compute_vertex_lon<CT, CS_Tag>::apply(lat1, lat2, 313 vertex_lat, 314 sin(lon1 - lon2), 315 cos(lon1 - lon2), 316 alp1, 317 azimuth_strategy); 318 319 CT vertex_lon = std::fmod(lon1 + dlon, 2 * pi); 320 321 if (vertex_lat < c0) 322 { 323 vertex_lon -= pi; 324 } 325 326 if (std::abs(lon1 - lon2) > pi) 327 { 328 vertex_lon -= pi; 329 } 330 331 return vertex_lon; 332 } 333 }; 334 335 template <typename CT> 336 class vertex_longitude<CT, cartesian_tag> 337 { 338 public : 339 template <typename Strategy> apply(CT &,CT &,CT & lon2,CT &,CT const &,CT &,Strategy const &)340 static inline CT apply(CT& /*lon1*/, 341 CT& /*lat1*/, 342 CT& lon2, 343 CT& /*lat2*/, 344 CT const& /*vertex_lat*/, 345 CT& /*alp1*/, 346 Strategy const& /*azimuth_strategy*/) 347 { 348 return lon2; 349 } 350 }; 351 352 }}} // namespace boost::geometry::formula 353 #endif // BOOST_GEOMETRY_FORMULAS_MAXIMUM_LONGITUDE_HPP 354 355