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

   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_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 #include <boost/geometry/core/srs.hpp>
  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>
  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         int sign = c1;
 213         if (bet3 < c0)
 214         {
 215             sign = cminus1;
 216         }
 217
 218         CT const dlon_max = omg13 - sign * f * sin_alp0 * I3;
 219
 220         return dlon_max;
 221     }
 222 };
 223
 224 //CS_tag dispatching
 225
 226 template <typename CT, typename CS_Tag>
 227 struct compute_vertex_lon
 228 {
 229     BOOST_MPL_ASSERT_MSG
 230     (
 231         false, NOT_IMPLEMENTED_FOR_THIS_COORDINATE_SYSTEM, (types<CS_Tag>)
 232     );
 233
 234 };
 235
 236 template <typename CT>
 237 struct compute_vertex_lon<CT, spherical_equatorial_tag>
 238 {
 239     template <typename Strategy>
 240     static inline CT apply(CT const& lat1,
 241                            CT const& lat2,
 242                            CT const& vertex_lat,
 243                            CT const& sin_l12,
 244                            CT const& cos_l12,
 245                            CT,
 246                            Strategy)
 247     {
 248         return vertex_longitude_on_sphere<CT>
 249                 ::apply(lat1,
 250                         lat2,
 251                         vertex_lat,
 252                         sin_l12,
 253                         cos_l12);
 254     }
 255 };
 256
 257 template <typename CT>
 258 struct compute_vertex_lon<CT, geographic_tag>
 259 {
 260     template <typename Strategy>
 261     static inline CT apply(CT const& lat1,
 262                            CT const& lat2,
 263                            CT const& vertex_lat,
 264                            CT,
 265                            CT,
 266                            CT& alp1,
 267                            Strategy const& azimuth_strategy)
 268     {
 269         return vertex_longitude_on_spheroid<CT>
 270                 ::apply(lat1,
 271                         lat2,
 272                         vertex_lat,
 273                         alp1,
 274                         azimuth_strategy.model());
 275     }
 276 };
 277
 278 // Vertex longitude interface
 279 // Assume that lon1 < lon2 and vertex_lat is the latitude of the vertex
 280
 281 template <typename CT, typename CS_Tag>
 282 class vertex_longitude
 283 {
 284 public :
 285     template <typename Strategy>
 286     static inline CT apply(CT& lon1,
 287                            CT& lat1,
 288                            CT& lon2,
 289                            CT& lat2,
 290                            CT const& vertex_lat,
 291                            CT& alp1,
 292                            Strategy const& azimuth_strategy)
 293     {
 294         CT const c0 = 0;
 295         CT pi = math::pi<CT>();
 296
 297         //Vertex is a segment's point
 298         if (math::equals(vertex_lat, lat1))
 299         {
 300             return lon1;
 301         }
 302         if (math::equals(vertex_lat, lat2))
 303         {
 304             return lon2;
 305         }
 306
 307         //Segment lay on meridian
 308         if (math::equals(lon1, lon2))
 309         {
 310             return (std::max)(lat1, lat2);
 311         }
 312         BOOST_ASSERT(lon1 < lon2);
 313
 314         CT dlon = compute_vertex_lon<CT, CS_Tag>::apply(lat1, lat2,
 315                                                         vertex_lat,
 316                                                         sin(lon1 - lon2),
 317                                                         cos(lon1 - lon2),
 318                                                         alp1,
 319                                                         azimuth_strategy);
 320
 321         CT vertex_lon = std::fmod(lon1 + dlon, 2 * pi);
 322
 323         if (vertex_lat < c0)
 324         {
 325             vertex_lon -= pi;
 326         }
 327
 328         if (std::abs(lon1 - lon2) > pi)
 329         {
 330             vertex_lon -= pi;
 331         }
 332
 333         return vertex_lon;
 334     }
 335 };
 336
 337 }}} // namespace boost::geometry::formula
 338 #endif // BOOST_GEOMETRY_FORMULAS_MAXIMUM_LONGITUDE_HPP
 339