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>
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         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>
applyboost::geometry::formula::compute_vertex_lon240     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>
applyboost::geometry::formula::compute_vertex_lon261     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>
apply(CT & lon1,CT & lat1,CT & lon2,CT & lat2,CT const & vertex_lat,CT & alp1,Strategy const & azimuth_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 
340