geometry/formulas/vertex_longitude.hpp

0001 // Boost.Geometry
0002
0003 // Copyright (c) 2016-2020 Oracle and/or its affiliates.
0004
0005 // Contributed and/or modified by Vissarion Fysikopoulos, on behalf of Oracle
0006 // Contributed and/or modified by Adam Wulkiewicz, on behalf of Oracle
0007
0008 // Use, modification and distribution is subject to the Boost Software License,
0009 // Version 1.0. (See accompanying file LICENSE_1_0.txt or copy at
0010 // http://www.boost.org/LICENSE_1_0.txt)
0011
0012 #ifndef BOOST_GEOMETRY_FORMULAS_MAXIMUM_LONGITUDE_HPP
0013 #define BOOST_GEOMETRY_FORMULAS_MAXIMUM_LONGITUDE_HPP
0014
0015
0016 #include <boost/geometry/core/static_assert.hpp>
0017 #include <boost/geometry/formulas/spherical.hpp>
0018 #include <boost/geometry/formulas/flattening.hpp>
0019
0020 #include <boost/math/special_functions/hypot.hpp>
0021
0022 namespace boost { namespace geometry { namespace formula
0023 {
0024
0025 /*!
0026 \brief Algorithm to compute the vertex longitude of a geodesic segment. Vertex is
0027 a point on the geodesic that maximizes (or minimizes) the latitude. The algorithm
0028 is given the vertex latitude.
0029 */
0030
0031 //Classes for spesific CS
0032
0033 template <typename CT>
0034 class vertex_longitude_on_sphere
0035 {
0036
0037 public:
0038
0039     template <typename T>
0040     static inline CT apply(T const& lat1, //segment point 1
0041                            T const& lat2, //segment point 2
0042                            T const& lat3, //vertex latitude
0043                            T const& sin_l12,
0044                            T const& cos_l12) //lon1 -lon2
0045     {
0046         //https://en.wikipedia.org/wiki/Great-circle_navigation#Finding_way-points
0047         CT const A = sin(lat1) * cos(lat2) * cos(lat3) * sin_l12;
0048         CT const B = sin(lat1) * cos(lat2) * cos(lat3) * cos_l12
0049                 - cos(lat1) * sin(lat2) * cos(lat3);
0050         CT lon = atan2(B, A);
0051         return lon + math::pi<CT>();
0052     }
0053 };
0054
0055 template <typename CT>
0056 class vertex_longitude_on_spheroid
0057 {
0058     template<typename T>
0059     static inline void normalize(T& x, T& y)
0060     {
0061         T h = boost::math::hypot(x, y);
0062         x /= h;
0063         y /= h;
0064     }
0065
0066 public:
0067
0068     template <typename T, typename Spheroid>
0069     static inline CT apply(T const& lat1, //segment point 1
0070                            T const& lat2, //segment point 2
0071                            T const& lat3, //vertex latitude
0072                            T& alp1,
0073                            Spheroid const& spheroid)
0074     {
0075         // We assume that segment points lay on different side w.r.t.
0076         // the vertex
0077
0078         // Constants
0079         CT const c0 = 0;
0080         CT const c2 = 2;
0081         CT const half_pi = math::pi<CT>() / c2;
0082         if (math::equals(lat1, half_pi)
0083                 || math::equals(lat2, half_pi)
0084                 || math::equals(lat1, -half_pi)
0085                 || math::equals(lat2, -half_pi))
0086         {
0087             // one segment point is the pole
0088             return c0;
0089         }
0090
0091         // More constants
0092         CT const f = flattening<CT>(spheroid);
0093         CT const pi = math::pi<CT>();
0094         CT const c1 = 1;
0095         CT const cminus1 = -1;
0096
0097         // First, compute longitude on auxiliary sphere
0098
0099         CT const one_minus_f = c1 - f;
0100         CT const bet1 = atan(one_minus_f * tan(lat1));
0101         CT const bet2 = atan(one_minus_f * tan(lat2));
0102         CT const bet3 = atan(one_minus_f * tan(lat3));
0103
0104         CT cos_bet1 = cos(bet1);
0105         CT cos_bet2 = cos(bet2);
0106         CT const sin_bet1 = sin(bet1);
0107         CT const sin_bet2 = sin(bet2);
0108         CT const sin_bet3 = sin(bet3);
0109
0110         CT omg12 = 0;
0111
0112         if (bet1 < c0)
0113         {
0114             cos_bet1 *= cminus1;
0115             omg12 += pi;
0116         }
0117         if (bet2 < c0)
0118         {
0119             cos_bet2 *= cminus1;
0120             omg12 += pi;
0121         }
0122
0123         CT const sin_alp1 = sin(alp1);
0124         CT const cos_alp1 = math::sqrt(c1 - math::sqr(sin_alp1));
0125
0126         CT const norm = math::sqrt(math::sqr(cos_alp1) + math::sqr(sin_alp1 * sin_bet1));
0127         CT const sin_alp0 = sin(atan2(sin_alp1 * cos_bet1, norm));
0128
0129         BOOST_ASSERT(cos_bet2 != c0);
0130         CT const sin_alp2 = sin_alp1 * cos_bet1 / cos_bet2;
0131
0132         CT const cos_alp0 = math::sqrt(c1 - math::sqr(sin_alp0));
0133         CT const cos_alp2 = math::sqrt(c1 - math::sqr(sin_alp2));
0134
0135         CT const sig1 = atan2(sin_bet1, cos_alp1 * cos_bet1);
0136         CT const sig2 = atan2(sin_bet2, -cos_alp2 * cos_bet2); //lat3 is a vertex
0137
0138         CT const cos_sig1 = cos(sig1);
0139         CT const sin_sig1 = math::sqrt(c1 - math::sqr(cos_sig1));
0140
0141         CT const cos_sig2 = cos(sig2);
0142         CT const sin_sig2 = math::sqrt(c1 - math::sqr(cos_sig2));
0143
0144         CT const omg1 = atan2(sin_alp0 * sin_sig1, cos_sig1);
0145         CT const omg2 = atan2(sin_alp0 * sin_sig2, cos_sig2);
0146
0147         omg12 += omg1 - omg2;
0148
0149         CT const sin_omg12 = sin(omg12);
0150         CT const cos_omg12 = cos(omg12);
0151
0152         CT omg13 = geometry::formula::vertex_longitude_on_sphere<CT>
0153                 ::apply(bet1, bet2, bet3, sin_omg12, cos_omg12);
0154
0155         if (lat1 * lat2 < c0)//different hemispheres
0156         {
0157             if ((lat2 - lat1) * lat3  > c0)// ascending segment
0158             {
0159                 omg13 = pi - omg13;
0160             }
0161         }
0162
0163         // Second, compute the ellipsoidal longitude
0164
0165         CT const e2 = f * (c2 - f);
0166         CT const ep = math::sqrt(e2 / (c1 - e2));
0167         CT const k2 = math::sqr(ep * cos_alp0);
0168         CT const sqrt_k2_plus_one = math::sqrt(c1 + k2);
0169         CT const eps = (sqrt_k2_plus_one - c1) / (sqrt_k2_plus_one + c1);
0170         CT const eps2 = eps * eps;
0171         CT const n = f / (c2 - f);
0172
0173         // sig3 is the length from equator to the vertex
0174         CT sig3;
0175         if(sin_bet3 > c0)
0176         {
0177             sig3 = half_pi;
0178         } else {
0179             sig3 = -half_pi;
0180         }
0181         CT const cos_sig3 = 0;
0182         CT const sin_sig3 = 1;
0183
0184         CT sig13 = sig3 - sig1;
0185         if (sig13 > pi)
0186         {
0187             sig13 -= 2 * pi;
0188         }
0189
0190         // Order 2 approximation
0191         CT const c1over2 = 0.5;
0192         CT const c1over4 = 0.25;
0193         CT const c1over8 = 0.125;
0194         CT const c1over16 = 0.0625;
0195         CT const c4 = 4;
0196         CT const c8 = 8;
0197
0198         CT const A3 = 1 - (c1over2 - c1over2 * n) * eps - c1over4 * eps2;
0199         CT const C31 = (c1over4 - c1over4 * n) * eps + c1over8 * eps2;
0200         CT const C32 = c1over16 * eps2;
0201
0202         CT const sin2_sig3 = c2 * cos_sig3 * sin_sig3;
0203         CT const sin4_sig3 = sin_sig3 * (-c4 * cos_sig3
0204                                          + c8 * cos_sig3 * cos_sig3 * cos_sig3);
0205         CT const sin2_sig1 = c2 * cos_sig1 * sin_sig1;
0206         CT const sin4_sig1 = sin_sig1 * (-c4 * cos_sig1
0207                                          + c8 * cos_sig1 * cos_sig1 * cos_sig1);
0208         CT const I3 = A3 * (sig13
0209                             + C31 * (sin2_sig3 - sin2_sig1)
0210                             + C32 * (sin4_sig3 - sin4_sig1));
0211
0212         CT const sign = bet3 >= c0
0213                       ? c1
0214                       : cminus1;
0215
0216         CT const dlon_max = omg13 - sign * f * sin_alp0 * I3;
0217
0218         return dlon_max;
0219     }
0220 };
0221
0222 //CS_tag dispatching
0223
0224 template <typename CT, typename CS_Tag>
0225 struct compute_vertex_lon
0226 {
0227     BOOST_GEOMETRY_STATIC_ASSERT_FALSE(
0228         "Not implemented for this coordinate system.",
0229         CT, CS_Tag);
0230 };
0231
0232 template <typename CT>
0233 struct compute_vertex_lon<CT, spherical_equatorial_tag>
0234 {
0235     template <typename Strategy>
0236     static inline CT apply(CT const& lat1,
0237                            CT const& lat2,
0238                            CT const& vertex_lat,
0239                            CT const& sin_l12,
0240                            CT const& cos_l12,
0241                            CT,
0242                            Strategy)
0243     {
0244         return vertex_longitude_on_sphere<CT>
0245                 ::apply(lat1,
0246                         lat2,
0247                         vertex_lat,
0248                         sin_l12,
0249                         cos_l12);
0250     }
0251 };
0252
0253 template <typename CT>
0254 struct compute_vertex_lon<CT, geographic_tag>
0255 {
0256     template <typename Strategy>
0257     static inline CT apply(CT const& lat1,
0258                            CT const& lat2,
0259                            CT const& vertex_lat,
0260                            CT,
0261                            CT,
0262                            CT& alp1,
0263                            Strategy const& azimuth_strategy)
0264     {
0265         return vertex_longitude_on_spheroid<CT>
0266                 ::apply(lat1,
0267                         lat2,
0268                         vertex_lat,
0269                         alp1,
0270                         azimuth_strategy.model());
0271     }
0272 };
0273
0274 // Vertex longitude interface
0275 // Assume that lon1 < lon2 and vertex_lat is the latitude of the vertex
0276
0277 template <typename CT, typename CS_Tag>
0278 class vertex_longitude
0279 {
0280 public :
0281     template <typename Strategy>
0282     static inline CT apply(CT& lon1,
0283                            CT& lat1,
0284                            CT& lon2,
0285                            CT& lat2,
0286                            CT const& vertex_lat,
0287                            CT& alp1,
0288                            Strategy const& azimuth_strategy)
0289     {
0290         CT const c0 = 0;
0291         CT pi = math::pi<CT>();
0292
0293         //Vertex is a segment's point
0294         if (math::equals(vertex_lat, lat1))
0295         {
0296             return lon1;
0297         }
0298         if (math::equals(vertex_lat, lat2))
0299         {
0300             return lon2;
0301         }
0302
0303         //Segment lay on meridian
0304         if (math::equals(lon1, lon2))
0305         {
0306             return (std::max)(lat1, lat2);
0307         }
0308         BOOST_ASSERT(lon1 < lon2);
0309
0310         CT dlon = compute_vertex_lon<CT, CS_Tag>::apply(lat1, lat2,
0311                                                         vertex_lat,
0312                                                         sin(lon1 - lon2),
0313                                                         cos(lon1 - lon2),
0314                                                         alp1,
0315                                                         azimuth_strategy);
0316
0317         CT vertex_lon = std::fmod(lon1 + dlon, 2 * pi);
0318
0319         if (vertex_lat < c0)
0320         {
0321             vertex_lon -= pi;
0322         }
0323
0324         if (std::abs(lon1 - lon2) > pi)
0325         {
0326             vertex_lon -= pi;
0327         }
0328
0329         return vertex_lon;
0330     }
0331 };
0332
0333 template <typename CT>
0334 class vertex_longitude<CT, cartesian_tag>
0335 {
0336 public :
0337     template <typename Strategy>
0338     static inline CT apply(CT& /*lon1*/,
0339                            CT& /*lat1*/,
0340                            CT& lon2,
0341                            CT& /*lat2*/,
0342                            CT const& /*vertex_lat*/,
0343                            CT& /*alp1*/,
0344                            Strategy const& /*azimuth_strategy*/)
0345     {
0346         return lon2;
0347     }
0348 };
0349
0350 }}} // namespace boost::geometry::formula
0351 #endif // BOOST_GEOMETRY_FORMULAS_MAXIMUM_LONGITUDE_HPP
0352