]>
Commit | Line | Data |
---|---|---|
b32b8144 FG |
1 | // Boost.Geometry |
2 | ||
20effc67 | 3 | // Copyright (c) 2016-2020 Oracle and/or its affiliates. |
b32b8144 FG |
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 | ||
20effc67 TL |
15 | |
16 | #include <boost/geometry/core/static_assert.hpp> | |
b32b8144 FG |
17 | #include <boost/geometry/formulas/spherical.hpp> |
18 | #include <boost/geometry/formulas/flattening.hpp> | |
11fdf7f2 | 19 | |
b32b8144 FG |
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 | ||
11fdf7f2 TL |
212 | CT const sign = bet3 >= c0 |
213 | ? c1 | |
214 | : cminus1; | |
215 | ||
b32b8144 FG |
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 | { | |
20effc67 TL |
227 | BOOST_GEOMETRY_STATIC_ASSERT_FALSE( |
228 | "Not implemented for this coordinate system.", | |
229 | CT, CS_Tag); | |
b32b8144 FG |
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 | ||
92f5a8d4 TL |
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 | ||
b32b8144 FG |
350 | }}} // namespace boost::geometry::formula |
351 | #endif // BOOST_GEOMETRY_FORMULAS_MAXIMUM_LONGITUDE_HPP | |
352 |