]>
Commit | Line | Data |
---|---|---|
b32b8144 FG |
1 | // Boost.Geometry |
2 | ||
20effc67 | 3 | // Copyright (c) 2016-2020, Oracle and/or its affiliates. |
b32b8144 FG |
4 | // Contributed and/or modified by Vissarion Fysikopoulos, on behalf of Oracle |
5 | // Contributed and/or modified by Adam Wulkiewicz, on behalf of Oracle | |
6 | ||
7 | // Use, modification and distribution is subject to the Boost Software License, | |
8 | // Version 1.0. (See accompanying file LICENSE_1_0.txt or copy at | |
9 | // http://www.boost.org/LICENSE_1_0.txt) | |
10 | ||
11 | #ifndef BOOST_GEOMETRY_FORMULAS_SPHERICAL_HPP | |
12 | #define BOOST_GEOMETRY_FORMULAS_SPHERICAL_HPP | |
13 | ||
14 | #include <boost/geometry/core/coordinate_system.hpp> | |
15 | #include <boost/geometry/core/coordinate_type.hpp> | |
92f5a8d4 | 16 | #include <boost/geometry/core/cs.hpp> |
b32b8144 FG |
17 | #include <boost/geometry/core/access.hpp> |
18 | #include <boost/geometry/core/radian_access.hpp> | |
92f5a8d4 | 19 | #include <boost/geometry/core/radius.hpp> |
b32b8144 FG |
20 | |
21 | //#include <boost/geometry/arithmetic/arithmetic.hpp> | |
22 | #include <boost/geometry/arithmetic/cross_product.hpp> | |
23 | #include <boost/geometry/arithmetic/dot_product.hpp> | |
24 | ||
25 | #include <boost/geometry/util/math.hpp> | |
26 | #include <boost/geometry/util/normalize_spheroidal_coordinates.hpp> | |
27 | #include <boost/geometry/util/select_coordinate_type.hpp> | |
28 | ||
92f5a8d4 TL |
29 | #include <boost/geometry/formulas/result_direct.hpp> |
30 | ||
b32b8144 FG |
31 | namespace boost { namespace geometry { |
32 | ||
33 | namespace formula { | |
34 | ||
35 | template <typename T> | |
36 | struct result_spherical | |
37 | { | |
38 | result_spherical() | |
39 | : azimuth(0) | |
40 | , reverse_azimuth(0) | |
41 | {} | |
42 | ||
43 | T azimuth; | |
44 | T reverse_azimuth; | |
45 | }; | |
46 | ||
47 | template <typename T> | |
48 | static inline void sph_to_cart3d(T const& lon, T const& lat, T & x, T & y, T & z) | |
49 | { | |
50 | T const cos_lat = cos(lat); | |
51 | x = cos_lat * cos(lon); | |
52 | y = cos_lat * sin(lon); | |
53 | z = sin(lat); | |
54 | } | |
55 | ||
56 | template <typename Point3d, typename PointSph> | |
57 | static inline Point3d sph_to_cart3d(PointSph const& point_sph) | |
58 | { | |
59 | typedef typename coordinate_type<Point3d>::type calc_t; | |
60 | ||
61 | calc_t const lon = get_as_radian<0>(point_sph); | |
62 | calc_t const lat = get_as_radian<1>(point_sph); | |
63 | calc_t x, y, z; | |
64 | sph_to_cart3d(lon, lat, x, y, z); | |
65 | ||
66 | Point3d res; | |
67 | set<0>(res, x); | |
68 | set<1>(res, y); | |
69 | set<2>(res, z); | |
70 | ||
71 | return res; | |
72 | } | |
73 | ||
74 | template <typename T> | |
75 | static inline void cart3d_to_sph(T const& x, T const& y, T const& z, T & lon, T & lat) | |
76 | { | |
77 | lon = atan2(y, x); | |
78 | lat = asin(z); | |
79 | } | |
80 | ||
81 | template <typename PointSph, typename Point3d> | |
82 | static inline PointSph cart3d_to_sph(Point3d const& point_3d) | |
83 | { | |
84 | typedef typename coordinate_type<PointSph>::type coord_t; | |
85 | typedef typename coordinate_type<Point3d>::type calc_t; | |
86 | ||
87 | calc_t const x = get<0>(point_3d); | |
88 | calc_t const y = get<1>(point_3d); | |
89 | calc_t const z = get<2>(point_3d); | |
90 | calc_t lonr, latr; | |
91 | cart3d_to_sph(x, y, z, lonr, latr); | |
92 | ||
93 | PointSph res; | |
94 | set_from_radian<0>(res, lonr); | |
95 | set_from_radian<1>(res, latr); | |
96 | ||
97 | coord_t lon = get<0>(res); | |
98 | coord_t lat = get<1>(res); | |
99 | ||
100 | math::normalize_spheroidal_coordinates | |
101 | < | |
1e59de90 | 102 | typename geometry::detail::cs_angular_units<PointSph>::type, |
b32b8144 FG |
103 | coord_t |
104 | >(lon, lat); | |
105 | ||
106 | set<0>(res, lon); | |
107 | set<1>(res, lat); | |
108 | ||
109 | return res; | |
110 | } | |
111 | ||
112 | // -1 right | |
113 | // 1 left | |
114 | // 0 on | |
115 | template <typename Point3d1, typename Point3d2> | |
116 | static inline int sph_side_value(Point3d1 const& norm, Point3d2 const& pt) | |
117 | { | |
118 | typedef typename select_coordinate_type<Point3d1, Point3d2>::type calc_t; | |
119 | calc_t c0 = 0; | |
120 | calc_t d = dot_product(norm, pt); | |
121 | return math::equals(d, c0) ? 0 | |
122 | : d > c0 ? 1 | |
123 | : -1; // d < 0 | |
124 | } | |
125 | ||
126 | template <typename CT, bool ReverseAzimuth, typename T1, typename T2> | |
127 | static inline result_spherical<CT> spherical_azimuth(T1 const& lon1, | |
128 | T1 const& lat1, | |
129 | T2 const& lon2, | |
130 | T2 const& lat2) | |
131 | { | |
132 | typedef result_spherical<CT> result_type; | |
133 | result_type result; | |
134 | ||
135 | // http://williams.best.vwh.net/avform.htm#Crs | |
136 | // https://en.wikipedia.org/wiki/Great-circle_navigation | |
137 | CT dlon = lon2 - lon1; | |
138 | ||
139 | // An optimization which should kick in often for Boxes | |
140 | //if ( math::equals(dlon, ReturnType(0)) ) | |
141 | //if ( get<0>(p1) == get<0>(p2) ) | |
142 | //{ | |
143 | // return - sin(get_as_radian<1>(p1)) * cos_p2lat); | |
144 | //} | |
145 | ||
146 | CT const cos_dlon = cos(dlon); | |
147 | CT const sin_dlon = sin(dlon); | |
148 | CT const cos_lat1 = cos(lat1); | |
149 | CT const cos_lat2 = cos(lat2); | |
150 | CT const sin_lat1 = sin(lat1); | |
151 | CT const sin_lat2 = sin(lat2); | |
152 | ||
153 | { | |
154 | // "An alternative formula, not requiring the pre-computation of d" | |
155 | // In the formula below dlon is used as "d" | |
156 | CT const y = sin_dlon * cos_lat2; | |
157 | CT const x = cos_lat1 * sin_lat2 - sin_lat1 * cos_lat2 * cos_dlon; | |
158 | result.azimuth = atan2(y, x); | |
159 | } | |
160 | ||
161 | if (ReverseAzimuth) | |
162 | { | |
163 | CT const y = sin_dlon * cos_lat1; | |
164 | CT const x = sin_lat2 * cos_lat1 * cos_dlon - cos_lat2 * sin_lat1; | |
165 | result.reverse_azimuth = atan2(y, x); | |
166 | } | |
167 | ||
168 | return result; | |
169 | } | |
170 | ||
171 | template <typename ReturnType, typename T1, typename T2> | |
172 | inline ReturnType spherical_azimuth(T1 const& lon1, T1 const& lat1, | |
173 | T2 const& lon2, T2 const& lat2) | |
174 | { | |
175 | return spherical_azimuth<ReturnType, false>(lon1, lat1, lon2, lat2).azimuth; | |
176 | } | |
177 | ||
178 | template <typename T> | |
179 | inline T spherical_azimuth(T const& lon1, T const& lat1, T const& lon2, T const& lat2) | |
180 | { | |
181 | return spherical_azimuth<T, false>(lon1, lat1, lon2, lat2).azimuth; | |
182 | } | |
183 | ||
184 | template <typename T> | |
185 | inline int azimuth_side_value(T const& azi_a1_p, T const& azi_a1_a2) | |
186 | { | |
92f5a8d4 | 187 | T const c0 = 0; |
b32b8144 | 188 | T const pi = math::pi<T>(); |
b32b8144 FG |
189 | |
190 | // instead of the formula from XTD | |
191 | //calc_t a_diff = asin(sin(azi_a1_p - azi_a1_a2)); | |
192 | ||
193 | T a_diff = azi_a1_p - azi_a1_a2; | |
20effc67 TL |
194 | // normalize, angle in (-pi, pi] |
195 | math::detail::normalize_angle_loop<radian>(a_diff); | |
b32b8144 FG |
196 | |
197 | // NOTE: in general it shouldn't be required to support the pi/-pi case | |
198 | // because in non-cartesian systems it makes sense to check the side | |
199 | // only "between" the endpoints. | |
200 | // However currently the winding strategy calls the side strategy | |
201 | // for vertical segments to check if the point is "between the endpoints. | |
202 | // This could be avoided since the side strategy is not required for that | |
203 | // because meridian is the shortest path. So a difference of | |
20effc67 | 204 | // longitudes would be sufficient (of course normalized to (-pi, pi]). |
b32b8144 FG |
205 | |
206 | // NOTE: with the above said, the pi/-pi check is temporary | |
207 | // however in case if this was required | |
208 | // the geodesics on ellipsoid aren't "symmetrical" | |
209 | // therefore instead of comparing a_diff to pi and -pi | |
210 | // one should probably use inverse azimuths and compare | |
211 | // the difference to 0 as well | |
212 | ||
213 | // positive azimuth is on the right side | |
92f5a8d4 | 214 | return math::equals(a_diff, c0) |
b32b8144 FG |
215 | || math::equals(a_diff, pi) |
216 | || math::equals(a_diff, -pi) ? 0 | |
217 | : a_diff > 0 ? -1 // right | |
218 | : 1; // left | |
219 | } | |
220 | ||
92f5a8d4 TL |
221 | template |
222 | < | |
223 | bool Coordinates, | |
224 | bool ReverseAzimuth, | |
225 | typename CT, | |
226 | typename Sphere | |
227 | > | |
228 | inline result_direct<CT> spherical_direct(CT const& lon1, | |
229 | CT const& lat1, | |
230 | CT const& sig12, | |
231 | CT const& alp1, | |
232 | Sphere const& sphere) | |
233 | { | |
234 | result_direct<CT> result; | |
235 | ||
236 | CT const sin_alp1 = sin(alp1); | |
237 | CT const sin_lat1 = sin(lat1); | |
238 | CT const cos_alp1 = cos(alp1); | |
239 | CT const cos_lat1 = cos(lat1); | |
240 | ||
241 | CT const norm = math::sqrt(cos_alp1 * cos_alp1 + sin_alp1 * sin_alp1 | |
242 | * sin_lat1 * sin_lat1); | |
243 | CT const alp0 = atan2(sin_alp1 * cos_lat1, norm); | |
244 | CT const sig1 = atan2(sin_lat1, cos_alp1 * cos_lat1); | |
245 | CT const sig2 = sig1 + sig12 / get_radius<0>(sphere); | |
246 | ||
247 | CT const cos_sig2 = cos(sig2); | |
248 | CT const sin_alp0 = sin(alp0); | |
249 | CT const cos_alp0 = cos(alp0); | |
250 | ||
251 | if (Coordinates) | |
252 | { | |
253 | CT const sin_sig2 = sin(sig2); | |
254 | CT const sin_sig1 = sin(sig1); | |
255 | CT const cos_sig1 = cos(sig1); | |
256 | ||
257 | CT const norm2 = math::sqrt(cos_alp0 * cos_alp0 * cos_sig2 * cos_sig2 | |
258 | + sin_alp0 * sin_alp0); | |
259 | CT const lat2 = atan2(cos_alp0 * sin_sig2, norm2); | |
260 | ||
261 | CT const omg1 = atan2(sin_alp0 * sin_sig1, cos_sig1); | |
262 | CT const lon2 = atan2(sin_alp0 * sin_sig2, cos_sig2); | |
263 | ||
264 | result.lon2 = lon1 + lon2 - omg1; | |
265 | result.lat2 = lat2; | |
20effc67 TL |
266 | |
267 | // For longitudes close to the antimeridian the result can be out | |
268 | // of range. Therefore normalize. | |
269 | math::detail::normalize_angle_cond<radian>(result.lon2); | |
92f5a8d4 TL |
270 | } |
271 | ||
272 | if (ReverseAzimuth) | |
273 | { | |
274 | CT const alp2 = atan2(sin_alp0, cos_alp0 * cos_sig2); | |
275 | result.reverse_azimuth = alp2; | |
276 | } | |
277 | ||
278 | return result; | |
279 | } | |
280 | ||
b32b8144 FG |
281 | } // namespace formula |
282 | ||
283 | }} // namespace boost::geometry | |
284 | ||
285 | #endif // BOOST_GEOMETRY_FORMULAS_SPHERICAL_HPP |