]>
Commit | Line | Data |
---|---|---|
b32b8144 FG |
1 | // Boost.Geometry |
2 | ||
92f5a8d4 | 3 | // Copyright (c) 2016-2019 Oracle and/or its affiliates. |
b32b8144 FG |
4 | |
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_INVERSE_DIFFERENTIAL_QUANTITIES_HPP | |
12 | #define BOOST_GEOMETRY_FORMULAS_INVERSE_DIFFERENTIAL_QUANTITIES_HPP | |
13 | ||
92f5a8d4 | 14 | #include <boost/geometry/core/assert.hpp> |
b32b8144 FG |
15 | |
16 | #include <boost/geometry/util/condition.hpp> | |
17 | #include <boost/geometry/util/math.hpp> | |
18 | ||
19 | ||
20 | namespace boost { namespace geometry { namespace formula | |
21 | { | |
22 | ||
23 | /*! | |
24 | \brief The solution of a part of the inverse problem - differential quantities. | |
25 | \author See | |
26 | - Charles F.F Karney, Algorithms for geodesics, 2011 | |
27 | https://arxiv.org/pdf/1109.4448.pdf | |
28 | */ | |
29 | template < | |
30 | typename CT, | |
31 | bool EnableReducedLength, | |
32 | bool EnableGeodesicScale, | |
33 | unsigned int Order = 2, | |
34 | bool ApproxF = true | |
35 | > | |
36 | class differential_quantities | |
37 | { | |
38 | public: | |
39 | static inline void apply(CT const& lon1, CT const& lat1, | |
40 | CT const& lon2, CT const& lat2, | |
41 | CT const& azimuth, CT const& reverse_azimuth, | |
42 | CT const& b, CT const& f, | |
43 | CT & reduced_length, CT & geodesic_scale) | |
44 | { | |
45 | CT const dlon = lon2 - lon1; | |
46 | CT const sin_lat1 = sin(lat1); | |
47 | CT const cos_lat1 = cos(lat1); | |
48 | CT const sin_lat2 = sin(lat2); | |
49 | CT const cos_lat2 = cos(lat2); | |
50 | ||
51 | apply(dlon, sin_lat1, cos_lat1, sin_lat2, cos_lat2, | |
52 | azimuth, reverse_azimuth, | |
53 | b, f, | |
54 | reduced_length, geodesic_scale); | |
55 | } | |
56 | ||
57 | static inline void apply(CT const& dlon, | |
58 | CT const& sin_lat1, CT const& cos_lat1, | |
59 | CT const& sin_lat2, CT const& cos_lat2, | |
60 | CT const& azimuth, CT const& reverse_azimuth, | |
61 | CT const& b, CT const& f, | |
62 | CT & reduced_length, CT & geodesic_scale) | |
63 | { | |
64 | CT const c0 = 0; | |
65 | CT const c1 = 1; | |
66 | CT const one_minus_f = c1 - f; | |
67 | ||
68 | CT sin_bet1 = one_minus_f * sin_lat1; | |
69 | CT sin_bet2 = one_minus_f * sin_lat2; | |
70 | ||
71 | // equator | |
72 | if (math::equals(sin_bet1, c0) && math::equals(sin_bet2, c0)) | |
73 | { | |
92f5a8d4 | 74 | CT const sig_12 = dlon / one_minus_f; |
b32b8144 FG |
75 | if (BOOST_GEOMETRY_CONDITION(EnableReducedLength)) |
76 | { | |
92f5a8d4 TL |
77 | BOOST_GEOMETRY_ASSERT((-math::pi<CT>() <= azimuth && azimuth <= math::pi<CT>())); |
78 | ||
79 | int azi_sign = math::sign(azimuth) >= 0 ? 1 : -1; // for antipodal | |
80 | CT m12 = azi_sign * sin(sig_12) * b; | |
b32b8144 FG |
81 | reduced_length = m12; |
82 | } | |
83 | ||
84 | if (BOOST_GEOMETRY_CONDITION(EnableGeodesicScale)) | |
85 | { | |
86 | CT M12 = cos(sig_12); | |
87 | geodesic_scale = M12; | |
88 | } | |
89 | } | |
90 | else | |
91 | { | |
92 | CT const c2 = 2; | |
93 | CT const e2 = f * (c2 - f); | |
94 | CT const ep2 = e2 / math::sqr(one_minus_f); | |
95 | ||
96 | CT const sin_alp1 = sin(azimuth); | |
97 | CT const cos_alp1 = cos(azimuth); | |
98 | //CT const sin_alp2 = sin(reverse_azimuth); | |
99 | CT const cos_alp2 = cos(reverse_azimuth); | |
100 | ||
101 | CT cos_bet1 = cos_lat1; | |
102 | CT cos_bet2 = cos_lat2; | |
103 | ||
104 | normalize(sin_bet1, cos_bet1); | |
105 | normalize(sin_bet2, cos_bet2); | |
106 | ||
107 | CT sin_sig1 = sin_bet1; | |
108 | CT cos_sig1 = cos_alp1 * cos_bet1; | |
109 | CT sin_sig2 = sin_bet2; | |
110 | CT cos_sig2 = cos_alp2 * cos_bet2; | |
111 | ||
112 | normalize(sin_sig1, cos_sig1); | |
113 | normalize(sin_sig2, cos_sig2); | |
114 | ||
115 | CT const sin_alp0 = sin_alp1 * cos_bet1; | |
116 | CT const cos_alp0_sqr = c1 - math::sqr(sin_alp0); | |
117 | ||
118 | CT const J12 = BOOST_GEOMETRY_CONDITION(ApproxF) ? | |
119 | J12_f(sin_sig1, cos_sig1, sin_sig2, cos_sig2, cos_alp0_sqr, f) : | |
120 | J12_ep_sqr(sin_sig1, cos_sig1, sin_sig2, cos_sig2, cos_alp0_sqr, ep2) ; | |
121 | ||
122 | CT const dn1 = math::sqrt(c1 + ep2 * math::sqr(sin_bet1)); | |
123 | CT const dn2 = math::sqrt(c1 + ep2 * math::sqr(sin_bet2)); | |
124 | ||
125 | if (BOOST_GEOMETRY_CONDITION(EnableReducedLength)) | |
126 | { | |
127 | CT const m12_b = dn2 * (cos_sig1 * sin_sig2) | |
128 | - dn1 * (sin_sig1 * cos_sig2) | |
129 | - cos_sig1 * cos_sig2 * J12; | |
130 | CT const m12 = m12_b * b; | |
131 | ||
132 | reduced_length = m12; | |
133 | } | |
134 | ||
135 | if (BOOST_GEOMETRY_CONDITION(EnableGeodesicScale)) | |
136 | { | |
137 | CT const cos_sig12 = cos_sig1 * cos_sig2 + sin_sig1 * sin_sig2; | |
138 | CT const t = ep2 * (cos_bet1 - cos_bet2) * (cos_bet1 + cos_bet2) / (dn1 + dn2); | |
139 | CT const M12 = cos_sig12 + (t * sin_sig2 - cos_sig2 * J12) * sin_sig1 / dn1; | |
140 | ||
141 | geodesic_scale = M12; | |
142 | } | |
143 | } | |
144 | } | |
145 | ||
146 | private: | |
147 | /*! Approximation of J12, expanded into taylor series in f | |
148 | Maxima script: | |
149 | ep2: f * (2-f) / ((1-f)^2); | |
150 | k2: ca02 * ep2; | |
151 | assume(f < 1); | |
152 | assume(sig > 0); | |
153 | I1(sig):= integrate(sqrt(1 + k2 * sin(s)^2), s, 0, sig); | |
154 | I2(sig):= integrate(1/sqrt(1 + k2 * sin(s)^2), s, 0, sig); | |
155 | J(sig):= I1(sig) - I2(sig); | |
156 | S: taylor(J(sig), f, 0, 3); | |
157 | S1: factor( 2*integrate(sin(s)^2,s,0,sig)*ca02*f ); | |
158 | S2: factor( ((integrate(-6*ca02^2*sin(s)^4+6*ca02*sin(s)^2,s,0,sig)+integrate(-2*ca02^2*sin(s)^4+6*ca02*sin(s)^2,s,0,sig))*f^2)/4 ); | |
159 | S3: factor( ((integrate(30*ca02^3*sin(s)^6-54*ca02^2*sin(s)^4+24*ca02*sin(s)^2,s,0,sig)+integrate(6*ca02^3*sin(s)^6-18*ca02^2*sin(s)^4+24*ca02*sin(s)^2,s,0,sig))*f^3)/12 ); | |
160 | */ | |
161 | static inline CT J12_f(CT const& sin_sig1, CT const& cos_sig1, | |
162 | CT const& sin_sig2, CT const& cos_sig2, | |
163 | CT const& cos_alp0_sqr, CT const& f) | |
164 | { | |
165 | if (Order == 0) | |
166 | { | |
167 | return 0; | |
168 | } | |
169 | ||
170 | CT const c2 = 2; | |
171 | ||
172 | CT const sig_12 = atan2(cos_sig1 * sin_sig2 - sin_sig1 * cos_sig2, | |
173 | cos_sig1 * cos_sig2 + sin_sig1 * sin_sig2); | |
174 | CT const sin_2sig1 = c2 * cos_sig1 * sin_sig1; // sin(2sig1) | |
175 | CT const sin_2sig2 = c2 * cos_sig2 * sin_sig2; // sin(2sig2) | |
176 | CT const sin_2sig_12 = sin_2sig2 - sin_2sig1; | |
177 | CT const L1 = sig_12 - sin_2sig_12 / c2; | |
178 | ||
179 | if (Order == 1) | |
180 | { | |
181 | return cos_alp0_sqr * f * L1; | |
182 | } | |
183 | ||
184 | CT const sin_4sig1 = c2 * sin_2sig1 * (math::sqr(cos_sig1) - math::sqr(sin_sig1)); // sin(4sig1) | |
185 | CT const sin_4sig2 = c2 * sin_2sig2 * (math::sqr(cos_sig2) - math::sqr(sin_sig2)); // sin(4sig2) | |
186 | CT const sin_4sig_12 = sin_4sig2 - sin_4sig1; | |
187 | ||
188 | CT const c8 = 8; | |
189 | CT const c12 = 12; | |
190 | CT const c16 = 16; | |
191 | CT const c24 = 24; | |
192 | ||
193 | CT const L2 = -( cos_alp0_sqr * sin_4sig_12 | |
194 | + (-c8 * cos_alp0_sqr + c12) * sin_2sig_12 | |
195 | + (c12 * cos_alp0_sqr - c24) * sig_12) | |
196 | / c16; | |
197 | ||
198 | if (Order == 2) | |
199 | { | |
200 | return cos_alp0_sqr * f * (L1 + f * L2); | |
201 | } | |
202 | ||
203 | CT const c4 = 4; | |
204 | CT const c9 = 9; | |
205 | CT const c48 = 48; | |
206 | CT const c60 = 60; | |
207 | CT const c64 = 64; | |
208 | CT const c96 = 96; | |
209 | CT const c128 = 128; | |
210 | CT const c144 = 144; | |
211 | ||
212 | CT const cos_alp0_quad = math::sqr(cos_alp0_sqr); | |
213 | CT const sin3_2sig1 = math::sqr(sin_2sig1) * sin_2sig1; | |
214 | CT const sin3_2sig2 = math::sqr(sin_2sig2) * sin_2sig2; | |
215 | CT const sin3_2sig_12 = sin3_2sig2 - sin3_2sig1; | |
216 | ||
217 | CT const A = (c9 * cos_alp0_quad - c12 * cos_alp0_sqr) * sin_4sig_12; | |
218 | CT const B = c4 * cos_alp0_quad * sin3_2sig_12; | |
219 | CT const C = (-c48 * cos_alp0_quad + c96 * cos_alp0_sqr - c64) * sin_2sig_12; | |
220 | CT const D = (c60 * cos_alp0_quad - c144 * cos_alp0_sqr + c128) * sig_12; | |
221 | ||
222 | CT const L3 = (A + B + C + D) / c64; | |
223 | ||
224 | // Order 3 and higher | |
225 | return cos_alp0_sqr * f * (L1 + f * (L2 + f * L3)); | |
226 | } | |
227 | ||
228 | /*! Approximation of J12, expanded into taylor series in e'^2 | |
229 | Maxima script: | |
230 | k2: ca02 * ep2; | |
231 | assume(sig > 0); | |
232 | I1(sig):= integrate(sqrt(1 + k2 * sin(s)^2), s, 0, sig); | |
233 | I2(sig):= integrate(1/sqrt(1 + k2 * sin(s)^2), s, 0, sig); | |
234 | J(sig):= I1(sig) - I2(sig); | |
235 | S: taylor(J(sig), ep2, 0, 3); | |
236 | S1: factor( integrate(sin(s)^2,s,0,sig)*ca02*ep2 ); | |
237 | S2: factor( (integrate(sin(s)^4,s,0,sig)*ca02^2*ep2^2)/2 ); | |
238 | S3: factor( (3*integrate(sin(s)^6,s,0,sig)*ca02^3*ep2^3)/8 ); | |
239 | */ | |
240 | static inline CT J12_ep_sqr(CT const& sin_sig1, CT const& cos_sig1, | |
241 | CT const& sin_sig2, CT const& cos_sig2, | |
242 | CT const& cos_alp0_sqr, CT const& ep_sqr) | |
243 | { | |
244 | if (Order == 0) | |
245 | { | |
246 | return 0; | |
247 | } | |
248 | ||
249 | CT const c2 = 2; | |
250 | CT const c4 = 4; | |
251 | ||
252 | CT const c2a0ep2 = cos_alp0_sqr * ep_sqr; | |
253 | ||
254 | CT const sig_12 = atan2(cos_sig1 * sin_sig2 - sin_sig1 * cos_sig2, | |
255 | cos_sig1 * cos_sig2 + sin_sig1 * sin_sig2); // sig2 - sig1 | |
256 | CT const sin_2sig1 = c2 * cos_sig1 * sin_sig1; // sin(2sig1) | |
257 | CT const sin_2sig2 = c2 * cos_sig2 * sin_sig2; // sin(2sig2) | |
258 | CT const sin_2sig_12 = sin_2sig2 - sin_2sig1; | |
259 | ||
260 | CT const L1 = (c2 * sig_12 - sin_2sig_12) / c4; | |
261 | ||
262 | if (Order == 1) | |
263 | { | |
264 | return c2a0ep2 * L1; | |
265 | } | |
266 | ||
267 | CT const c8 = 8; | |
268 | CT const c64 = 64; | |
269 | ||
270 | CT const sin_4sig1 = c2 * sin_2sig1 * (math::sqr(cos_sig1) - math::sqr(sin_sig1)); // sin(4sig1) | |
271 | CT const sin_4sig2 = c2 * sin_2sig2 * (math::sqr(cos_sig2) - math::sqr(sin_sig2)); // sin(4sig2) | |
272 | CT const sin_4sig_12 = sin_4sig2 - sin_4sig1; | |
273 | ||
274 | CT const L2 = (sin_4sig_12 - c8 * sin_2sig_12 + 12 * sig_12) / c64; | |
275 | ||
276 | if (Order == 2) | |
277 | { | |
278 | return c2a0ep2 * (L1 + c2a0ep2 * L2); | |
279 | } | |
280 | ||
281 | CT const sin3_2sig1 = math::sqr(sin_2sig1) * sin_2sig1; | |
282 | CT const sin3_2sig2 = math::sqr(sin_2sig2) * sin_2sig2; | |
283 | CT const sin3_2sig_12 = sin3_2sig2 - sin3_2sig1; | |
284 | ||
285 | CT const c9 = 9; | |
286 | CT const c48 = 48; | |
287 | CT const c60 = 60; | |
288 | CT const c512 = 512; | |
289 | ||
290 | CT const L3 = (c9 * sin_4sig_12 + c4 * sin3_2sig_12 - c48 * sin_2sig_12 + c60 * sig_12) / c512; | |
291 | ||
292 | // Order 3 and higher | |
293 | return c2a0ep2 * (L1 + c2a0ep2 * (L2 + c2a0ep2 * L3)); | |
294 | } | |
295 | ||
296 | static inline void normalize(CT & x, CT & y) | |
297 | { | |
298 | CT const len = math::sqrt(math::sqr(x) + math::sqr(y)); | |
299 | x /= len; | |
300 | y /= len; | |
301 | } | |
302 | }; | |
303 | ||
304 | }}} // namespace boost::geometry::formula | |
305 | ||
306 | ||
307 | #endif // BOOST_GEOMETRY_FORMULAS_INVERSE_DIFFERENTIAL_QUANTITIES_HPP |