ceph/src/boost/boost/geometry/formulas/vertex_longitude.hpp

   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_LONGITUDE_HPP
  13 #define BOOST_GEOMETRY_FORMULAS_MAXIMUM_LONGITUDE_HPP
  14
  15
  16 #include <boost/geometry/core/static_assert.hpp>
  17 #include <boost/geometry/formulas/spherical.hpp>
  18 #include <boost/geometry/formulas/flattening.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>
  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>
  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>
  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_GEOMETRY_STATIC_ASSERT_FALSE(
 228         "Not implemented for this coordinate system.",
 229         CT, CS_Tag);
 230 };
 231
 232 template <typename CT>
 233 struct compute_vertex_lon<CT, spherical_equatorial_tag>
 234 {
 235     template <typename Strategy>
 236     static inline CT apply(CT const& lat1,
 237                            CT const& lat2,
 238                            CT const& vertex_lat,
 239                            CT const& sin_l12,
 240                            CT const& cos_l12,
 241                            CT,
 242                            Strategy)
 243     {
 244         return vertex_longitude_on_sphere<CT>
 245                 ::apply(lat1,
 246                         lat2,
 247                         vertex_lat,
 248                         sin_l12,
 249                         cos_l12);
 250     }
 251 };
 252
 253 template <typename CT>
 254 struct compute_vertex_lon<CT, geographic_tag>
 255 {
 256     template <typename Strategy>
 257     static inline CT apply(CT const& lat1,
 258                            CT const& lat2,
 259                            CT const& vertex_lat,
 260                            CT,
 261                            CT,
 262                            CT& alp1,
 263                            Strategy const& azimuth_strategy)
 264     {
 265         return vertex_longitude_on_spheroid<CT>
 266                 ::apply(lat1,
 267                         lat2,
 268                         vertex_lat,
 269                         alp1,
 270                         azimuth_strategy.model());
 271     }
 272 };
 273
 274 // Vertex longitude interface
 275 // Assume that lon1 < lon2 and vertex_lat is the latitude of the vertex
 276
 277 template <typename CT, typename CS_Tag>
 278 class vertex_longitude
 279 {
 280 public :
 281     template <typename Strategy>
 282     static inline CT apply(CT& lon1,
 283                            CT& lat1,
 284                            CT& lon2,
 285                            CT& lat2,
 286                            CT const& vertex_lat,
 287                            CT& alp1,
 288                            Strategy const& azimuth_strategy)
 289     {
 290         CT const c0 = 0;
 291         CT pi = math::pi<CT>();
 292
 293         //Vertex is a segment's point
 294         if (math::equals(vertex_lat, lat1))
 295         {
 296             return lon1;
 297         }
 298         if (math::equals(vertex_lat, lat2))
 299         {
 300             return lon2;
 301         }
 302
 303         //Segment lay on meridian
 304         if (math::equals(lon1, lon2))
 305         {
 306             return (std::max)(lat1, lat2);
 307         }
 308         BOOST_ASSERT(lon1 < lon2);
 309
 310         CT dlon = compute_vertex_lon<CT, CS_Tag>::apply(lat1, lat2,
 311                                                         vertex_lat,
 312                                                         sin(lon1 - lon2),
 313                                                         cos(lon1 - lon2),
 314                                                         alp1,
 315                                                         azimuth_strategy);
 316
 317         CT vertex_lon = std::fmod(lon1 + dlon, 2 * pi);
 318
 319         if (vertex_lat < c0)
 320         {
 321             vertex_lon -= pi;
 322         }
 323
 324         if (std::abs(lon1 - lon2) > pi)
 325         {
 326             vertex_lon -= pi;
 327         }
 328
 329         return vertex_lon;
 330     }
 331 };
 332
 333 template <typename CT>
 334 class vertex_longitude<CT, cartesian_tag>
 335 {
 336 public :
 337     template <typename Strategy>
 338     static inline CT apply(CT& /*lon1*/,
 339                            CT& /*lat1*/,
 340                            CT& lon2,
 341                            CT& /*lat2*/,
 342                            CT const& /*vertex_lat*/,
 343                            CT& /*alp1*/,
 344                            Strategy const& /*azimuth_strategy*/)
 345     {
 346         return lon2;
 347     }
 348 };
 349
 350 }}} // namespace boost::geometry::formula
 351 #endif // BOOST_GEOMETRY_FORMULAS_MAXIMUM_LONGITUDE_HPP
 352