]>
Commit | Line | Data |
---|---|---|
b32b8144 FG |
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 |