]>
Commit | Line | Data |
---|---|---|
b32b8144 FG |
1 | // Boost.Geometry |
2 | ||
3 | // Copyright (c) 2015-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_AREA_FORMULAS_HPP | |
13 | #define BOOST_GEOMETRY_FORMULAS_AREA_FORMULAS_HPP | |
14 | ||
15 | #include <boost/geometry/formulas/flattening.hpp> | |
16 | #include <boost/math/special_functions/hypot.hpp> | |
17 | ||
18 | namespace boost { namespace geometry { namespace formula | |
19 | { | |
20 | ||
21 | /*! | |
22 | \brief Formulas for computing spherical and ellipsoidal polygon area. | |
23 | The current class computes the area of the trapezoid defined by a segment | |
24 | the two meridians passing by the endpoints and the equator. | |
25 | \author See | |
26 | - Danielsen JS, The area under the geodesic. Surv Rev 30(232): | |
27 | 61–66, 1989 | |
28 | - Charles F.F Karney, Algorithms for geodesics, 2011 | |
29 | https://arxiv.org/pdf/1109.4448.pdf | |
30 | */ | |
31 | ||
32 | template < | |
33 | typename CT, | |
34 | std::size_t SeriesOrder = 2, | |
35 | bool ExpandEpsN = true | |
36 | > | |
37 | class area_formulas | |
38 | { | |
39 | ||
40 | public: | |
41 | ||
42 | //TODO: move the following to a more general space to be used by other | |
43 | // classes as well | |
44 | /* | |
45 | Evaluate the polynomial in x using Horner's method. | |
46 | */ | |
47 | template <typename NT, typename IteratorType> | |
48 | static inline NT horner_evaluate(NT x, | |
49 | IteratorType begin, | |
50 | IteratorType end) | |
51 | { | |
52 | NT result(0); | |
53 | IteratorType it = end; | |
54 | do | |
55 | { | |
56 | result = result * x + *--it; | |
57 | } | |
58 | while (it != begin); | |
59 | return result; | |
60 | } | |
61 | ||
62 | /* | |
63 | Clenshaw algorithm for summing trigonometric series | |
64 | https://en.wikipedia.org/wiki/Clenshaw_algorithm | |
65 | */ | |
66 | template <typename NT, typename IteratorType> | |
67 | static inline NT clenshaw_sum(NT cosx, | |
68 | IteratorType begin, | |
69 | IteratorType end) | |
70 | { | |
71 | IteratorType it = end; | |
72 | bool odd = true; | |
73 | CT b_k, b_k1(0), b_k2(0); | |
74 | do | |
75 | { | |
76 | CT c_k = odd ? *--it : NT(0); | |
77 | b_k = c_k + NT(2) * cosx * b_k1 - b_k2; | |
78 | b_k2 = b_k1; | |
79 | b_k1 = b_k; | |
80 | odd = !odd; | |
81 | } | |
82 | while (it != begin); | |
83 | ||
84 | return *begin + b_k1 * cosx - b_k2; | |
85 | } | |
86 | ||
87 | template<typename T> | |
88 | static inline void normalize(T& x, T& y) | |
89 | { | |
90 | T h = boost::math::hypot(x, y); | |
91 | x /= h; | |
92 | y /= h; | |
93 | } | |
94 | ||
95 | /* | |
96 | Generate and evaluate the series expansion of the following integral | |
97 | ||
98 | I4 = -integrate( (t(ep2) - t(k2*sin(sigma1)^2)) / (ep2 - k2*sin(sigma1)^2) | |
99 | * sin(sigma1)/2, sigma1, pi/2, sigma ) | |
100 | where | |
101 | ||
102 | t(x) = sqrt(1+1/x)*asinh(sqrt(x)) + x | |
103 | ||
104 | valid for ep2 and k2 small. We substitute k2 = 4 * eps / (1 - eps)^2 | |
105 | and ep2 = 4 * n / (1 - n)^2 and expand in eps and n. | |
106 | ||
107 | The resulting sum of the series is of the form | |
108 | ||
109 | sum(C4[l] * cos((2*l+1)*sigma), l, 0, maxpow-1) ) | |
110 | ||
111 | The above expansion is performed in Computer Algebra System Maxima. | |
112 | The C++ code (that yields the function evaluate_coeffs_n below) is generated | |
113 | by the following Maxima script and is based on script: | |
114 | http://geographiclib.sourceforge.net/html/geod.mac | |
115 | ||
116 | // Maxima script begin | |
117 | taylordepth:5$ | |
118 | ataylor(expr,var,ord):=expand(ratdisrep(taylor(expr,var,0,ord)))$ | |
119 | jtaylor(expr,var1,var2,ord):=block([zz],expand(subst([zz=1], | |
120 | ratdisrep(taylor(subst([var1=zz*var1,var2=zz*var2],expr),zz,0,ord)))))$ | |
121 | ||
122 | compute(maxpow):=block([int,t,intexp,area, x,ep2,k2], | |
123 | maxpow:maxpow-1, | |
124 | t : sqrt(1+1/x) * asinh(sqrt(x)) + x, | |
125 | int:-(tf(ep2) - tf(k2*sin(sigma)^2)) / (ep2 - k2*sin(sigma)^2) | |
126 | * sin(sigma)/2, | |
127 | int:subst([tf(ep2)=subst([x=ep2],t), | |
128 | tf(k2*sin(sigma)^2)=subst([x=k2*sin(sigma)^2],t)], | |
129 | int), | |
130 | int:subst([abs(sin(sigma))=sin(sigma)],int), | |
131 | int:subst([k2=4*eps/(1-eps)^2,ep2=4*n/(1-n)^2],int), | |
132 | intexp:jtaylor(int,n,eps,maxpow), | |
133 | area:trigreduce(integrate(intexp,sigma)), | |
134 | area:expand(area-subst(sigma=%pi/2,area)), | |
135 | for i:0 thru maxpow do C4[i]:coeff(area,cos((2*i+1)*sigma)), | |
136 | if expand(area-sum(C4[i]*cos((2*i+1)*sigma),i,0,maxpow)) # 0 | |
137 | then error("left over terms in I4"), | |
138 | 'done)$ | |
139 | ||
140 | printcode(maxpow):= | |
141 | block([tab2:" ",tab3:" "], | |
142 | print(" switch (SeriesOrder) {"), | |
143 | for nn:1 thru maxpow do block([c], | |
144 | print(concat(tab2,"case ",string(nn-1),":")), | |
145 | c:0, | |
146 | for m:0 thru nn-1 do block( | |
147 | [q:jtaylor(subst([n=n],C4[m]),n,eps,nn-1), | |
148 | linel:1200], | |
149 | for j:m thru nn-1 do ( | |
150 | print(concat(tab3,"coeffs_n[",c,"] = ", | |
151 | string(horner(coeff(q,eps,j))),";")), | |
152 | c:c+1) | |
153 | ), | |
154 | print(concat(tab3,"break;"))), | |
155 | print(" }"), | |
156 | 'done)$ | |
157 | ||
158 | maxpow:6$ | |
159 | compute(maxpow)$ | |
160 | printcode(maxpow)$ | |
161 | // Maxima script end | |
162 | ||
163 | In the resulting code we should replace each number x by CT(x) | |
164 | e.g. using the following scirpt: | |
165 | sed -e 's/[0-9]\+/CT(&)/g; s/\[CT(/\[/g; s/)\]/\]/g; | |
166 | s/case\sCT(/case /g; s/):/:/g' | |
167 | */ | |
168 | ||
169 | static inline void evaluate_coeffs_n(CT n, CT coeffs_n[]) | |
170 | { | |
171 | ||
172 | switch (SeriesOrder) { | |
173 | case 0: | |
174 | coeffs_n[0] = CT(2)/CT(3); | |
175 | break; | |
176 | case 1: | |
177 | coeffs_n[0] = (CT(10)-CT(4)*n)/CT(15); | |
178 | coeffs_n[1] = -CT(1)/CT(5); | |
179 | coeffs_n[2] = CT(1)/CT(45); | |
180 | break; | |
181 | case 2: | |
182 | coeffs_n[0] = (n*(CT(8)*n-CT(28))+CT(70))/CT(105); | |
183 | coeffs_n[1] = (CT(16)*n-CT(7))/CT(35); | |
184 | coeffs_n[2] = -CT(2)/CT(105); | |
185 | coeffs_n[3] = (CT(7)-CT(16)*n)/CT(315); | |
186 | coeffs_n[4] = -CT(2)/CT(105); | |
187 | coeffs_n[5] = CT(4)/CT(525); | |
188 | break; | |
189 | case 3: | |
190 | coeffs_n[0] = (n*(n*(CT(4)*n+CT(24))-CT(84))+CT(210))/CT(315); | |
191 | coeffs_n[1] = ((CT(48)-CT(32)*n)*n-CT(21))/CT(105); | |
192 | coeffs_n[2] = (-CT(32)*n-CT(6))/CT(315); | |
193 | coeffs_n[3] = CT(11)/CT(315); | |
194 | coeffs_n[4] = (n*(CT(32)*n-CT(48))+CT(21))/CT(945); | |
195 | coeffs_n[5] = (CT(64)*n-CT(18))/CT(945); | |
196 | coeffs_n[6] = -CT(1)/CT(105); | |
197 | coeffs_n[7] = (CT(12)-CT(32)*n)/CT(1575); | |
198 | coeffs_n[8] = -CT(8)/CT(1575); | |
199 | coeffs_n[9] = CT(8)/CT(2205); | |
200 | break; | |
201 | case 4: | |
202 | coeffs_n[0] = (n*(n*(n*(CT(16)*n+CT(44))+CT(264))-CT(924))+CT(2310))/CT(3465); | |
203 | coeffs_n[1] = (n*(n*(CT(48)*n-CT(352))+CT(528))-CT(231))/CT(1155); | |
204 | coeffs_n[2] = (n*(CT(1088)*n-CT(352))-CT(66))/CT(3465); | |
205 | coeffs_n[3] = (CT(121)-CT(368)*n)/CT(3465); | |
206 | coeffs_n[4] = CT(4)/CT(1155); | |
207 | coeffs_n[5] = (n*((CT(352)-CT(48)*n)*n-CT(528))+CT(231))/CT(10395); | |
208 | coeffs_n[6] = ((CT(704)-CT(896)*n)*n-CT(198))/CT(10395); | |
209 | coeffs_n[7] = (CT(80)*n-CT(99))/CT(10395); | |
210 | coeffs_n[8] = CT(4)/CT(1155); | |
211 | coeffs_n[9] = (n*(CT(320)*n-CT(352))+CT(132))/CT(17325); | |
212 | coeffs_n[10] = (CT(384)*n-CT(88))/CT(17325); | |
213 | coeffs_n[11] = -CT(8)/CT(1925); | |
214 | coeffs_n[12] = (CT(88)-CT(256)*n)/CT(24255); | |
215 | coeffs_n[13] = -CT(16)/CT(8085); | |
216 | coeffs_n[14] = CT(64)/CT(31185); | |
217 | break; | |
218 | case 5: | |
219 | coeffs_n[0] = (n*(n*(n*(n*(CT(100)*n+CT(208))+CT(572))+CT(3432))-CT(12012))+CT(30030)) | |
220 | /CT(45045); | |
221 | coeffs_n[1] = (n*(n*(n*(CT(64)*n+CT(624))-CT(4576))+CT(6864))-CT(3003))/CT(15015); | |
222 | coeffs_n[2] = (n*((CT(14144)-CT(10656)*n)*n-CT(4576))-CT(858))/CT(45045); | |
223 | coeffs_n[3] = ((-CT(224)*n-CT(4784))*n+CT(1573))/CT(45045); | |
224 | coeffs_n[4] = (CT(1088)*n+CT(156))/CT(45045); | |
225 | coeffs_n[5] = CT(97)/CT(15015); | |
226 | coeffs_n[6] = (n*(n*((-CT(64)*n-CT(624))*n+CT(4576))-CT(6864))+CT(3003))/CT(135135); | |
227 | coeffs_n[7] = (n*(n*(CT(5952)*n-CT(11648))+CT(9152))-CT(2574))/CT(135135); | |
228 | coeffs_n[8] = (n*(CT(5792)*n+CT(1040))-CT(1287))/CT(135135); | |
229 | coeffs_n[9] = (CT(468)-CT(2944)*n)/CT(135135); | |
230 | coeffs_n[10] = CT(1)/CT(9009); | |
231 | coeffs_n[11] = (n*((CT(4160)-CT(1440)*n)*n-CT(4576))+CT(1716))/CT(225225); | |
232 | coeffs_n[12] = ((CT(4992)-CT(8448)*n)*n-CT(1144))/CT(225225); | |
233 | coeffs_n[13] = (CT(1856)*n-CT(936))/CT(225225); | |
234 | coeffs_n[14] = CT(8)/CT(10725); | |
235 | coeffs_n[15] = (n*(CT(3584)*n-CT(3328))+CT(1144))/CT(315315); | |
236 | coeffs_n[16] = (CT(1024)*n-CT(208))/CT(105105); | |
237 | coeffs_n[17] = -CT(136)/CT(63063); | |
238 | coeffs_n[18] = (CT(832)-CT(2560)*n)/CT(405405); | |
239 | coeffs_n[19] = -CT(128)/CT(135135); | |
240 | coeffs_n[20] = CT(128)/CT(99099); | |
241 | break; | |
242 | } | |
243 | } | |
244 | ||
245 | /* | |
246 | Expand in k2 and ep2. | |
247 | */ | |
248 | static inline void evaluate_coeffs_ep(CT ep, CT coeffs_n[]) | |
249 | { | |
250 | switch (SeriesOrder) { | |
251 | case 0: | |
252 | coeffs_n[0] = CT(2)/CT(3); | |
253 | break; | |
254 | case 1: | |
255 | coeffs_n[0] = (CT(10)-ep)/CT(15); | |
256 | coeffs_n[1] = -CT(1)/CT(20); | |
257 | coeffs_n[2] = CT(1)/CT(180); | |
258 | break; | |
259 | case 2: | |
260 | coeffs_n[0] = (ep*(CT(4)*ep-CT(7))+CT(70))/CT(105); | |
261 | coeffs_n[1] = (CT(4)*ep-CT(7))/CT(140); | |
262 | coeffs_n[2] = CT(1)/CT(42); | |
263 | coeffs_n[3] = (CT(7)-CT(4)*ep)/CT(1260); | |
264 | coeffs_n[4] = -CT(1)/CT(252); | |
265 | coeffs_n[5] = CT(1)/CT(2100); | |
266 | break; | |
267 | case 3: | |
268 | coeffs_n[0] = (ep*((CT(12)-CT(8)*ep)*ep-CT(21))+CT(210))/CT(315); | |
269 | coeffs_n[1] = ((CT(12)-CT(8)*ep)*ep-CT(21))/CT(420); | |
270 | coeffs_n[2] = (CT(3)-CT(2)*ep)/CT(126); | |
271 | coeffs_n[3] = -CT(1)/CT(72); | |
272 | coeffs_n[4] = (ep*(CT(8)*ep-CT(12))+CT(21))/CT(3780); | |
273 | coeffs_n[5] = (CT(2)*ep-CT(3))/CT(756); | |
274 | coeffs_n[6] = CT(1)/CT(360); | |
275 | coeffs_n[7] = (CT(3)-CT(2)*ep)/CT(6300); | |
276 | coeffs_n[8] = -CT(1)/CT(1800); | |
277 | coeffs_n[9] = CT(1)/CT(17640); | |
278 | break; | |
279 | case 4: | |
280 | coeffs_n[0] = (ep*(ep*(ep*(CT(64)*ep-CT(88))+CT(132))-CT(231))+CT(2310))/CT(3465); | |
281 | coeffs_n[1] = (ep*(ep*(CT(64)*ep-CT(88))+CT(132))-CT(231))/CT(4620); | |
282 | coeffs_n[2] = (ep*(CT(16)*ep-CT(22))+CT(33))/CT(1386); | |
283 | coeffs_n[3] = (CT(8)*ep-CT(11))/CT(792); | |
284 | coeffs_n[4] = CT(1)/CT(110); | |
285 | coeffs_n[5] = (ep*((CT(88)-CT(64)*ep)*ep-CT(132))+CT(231))/CT(41580); | |
286 | coeffs_n[6] = ((CT(22)-CT(16)*ep)*ep-CT(33))/CT(8316); | |
287 | coeffs_n[7] = (CT(11)-CT(8)*ep)/CT(3960); | |
288 | coeffs_n[8] = -CT(1)/CT(495); | |
289 | coeffs_n[9] = (ep*(CT(16)*ep-CT(22))+CT(33))/CT(69300); | |
290 | coeffs_n[10] = (CT(8)*ep-CT(11))/CT(19800); | |
291 | coeffs_n[11] = CT(1)/CT(1925); | |
292 | coeffs_n[12] = (CT(11)-CT(8)*ep)/CT(194040); | |
293 | coeffs_n[13] = -CT(1)/CT(10780); | |
294 | coeffs_n[14] = CT(1)/CT(124740); | |
295 | break; | |
296 | case 5: | |
297 | coeffs_n[0] = (ep*(ep*(ep*((CT(832)-CT(640)*ep)*ep-CT(1144))+CT(1716))-CT(3003))+CT(30030))/CT(45045); | |
298 | coeffs_n[1] = (ep*(ep*((CT(832)-CT(640)*ep)*ep-CT(1144))+CT(1716))-CT(3003))/CT(60060); | |
299 | coeffs_n[2] = (ep*((CT(208)-CT(160)*ep)*ep-CT(286))+CT(429))/CT(18018); | |
300 | coeffs_n[3] = ((CT(104)-CT(80)*ep)*ep-CT(143))/CT(10296); | |
301 | coeffs_n[4] = (CT(13)-CT(10)*ep)/CT(1430); | |
302 | coeffs_n[5] = -CT(1)/CT(156); | |
303 | coeffs_n[6] = (ep*(ep*(ep*(CT(640)*ep-CT(832))+CT(1144))-CT(1716))+CT(3003))/CT(540540); | |
304 | coeffs_n[7] = (ep*(ep*(CT(160)*ep-CT(208))+CT(286))-CT(429))/CT(108108); | |
305 | coeffs_n[8] = (ep*(CT(80)*ep-CT(104))+CT(143))/CT(51480); | |
306 | coeffs_n[9] = (CT(10)*ep-CT(13))/CT(6435); | |
307 | coeffs_n[10] = CT(5)/CT(3276); | |
308 | coeffs_n[11] = (ep*((CT(208)-CT(160)*ep)*ep-CT(286))+CT(429))/CT(900900); | |
309 | coeffs_n[12] = ((CT(104)-CT(80)*ep)*ep-CT(143))/CT(257400); | |
310 | coeffs_n[13] = (CT(13)-CT(10)*ep)/CT(25025); | |
311 | coeffs_n[14] = -CT(1)/CT(2184); | |
312 | coeffs_n[15] = (ep*(CT(80)*ep-CT(104))+CT(143))/CT(2522520); | |
313 | coeffs_n[16] = (CT(10)*ep-CT(13))/CT(140140); | |
314 | coeffs_n[17] = CT(5)/CT(45864); | |
315 | coeffs_n[18] = (CT(13)-CT(10)*ep)/CT(1621620); | |
316 | coeffs_n[19] = -CT(1)/CT(58968); | |
317 | coeffs_n[20] = CT(1)/CT(792792); | |
318 | break; | |
319 | } | |
320 | } | |
321 | ||
322 | /* | |
323 | Given the set of coefficients coeffs1[] evaluate on var2 and return | |
324 | the set of coefficients coeffs2[] | |
325 | */ | |
326 | static inline void evaluate_coeffs_var2(CT var2, | |
327 | CT coeffs1[], | |
328 | CT coeffs2[]){ | |
329 | std::size_t begin(0), end(0); | |
330 | for(std::size_t i = 0; i <= SeriesOrder; i++){ | |
331 | end = begin + SeriesOrder + 1 - i; | |
332 | coeffs2[i] = ((i==0) ? CT(1) : pow(var2,int(i))) | |
333 | * horner_evaluate(var2, coeffs1 + begin, coeffs1 + end); | |
334 | begin = end; | |
335 | } | |
336 | } | |
337 | ||
338 | /* | |
339 | Compute the spherical excess of a geodesic (or shperical) segment | |
340 | */ | |
341 | template < | |
342 | bool LongSegment, | |
343 | typename PointOfSegment | |
344 | > | |
345 | static inline CT spherical(PointOfSegment const& p1, | |
346 | PointOfSegment const& p2) | |
347 | { | |
348 | CT excess; | |
349 | ||
350 | if(LongSegment) // not for segments parallel to equator | |
351 | { | |
352 | CT cbet1 = cos(geometry::get_as_radian<1>(p1)); | |
353 | CT sbet1 = sin(geometry::get_as_radian<1>(p1)); | |
354 | CT cbet2 = cos(geometry::get_as_radian<1>(p2)); | |
355 | CT sbet2 = sin(geometry::get_as_radian<1>(p2)); | |
356 | ||
357 | CT omg12 = geometry::get_as_radian<0>(p1) | |
358 | - geometry::get_as_radian<0>(p2); | |
359 | CT comg12 = cos(omg12); | |
360 | CT somg12 = sin(omg12); | |
361 | ||
362 | CT alp1 = atan2(cbet1 * sbet2 | |
363 | - sbet1 * cbet2 * comg12, | |
364 | cbet2 * somg12); | |
365 | ||
366 | CT alp2 = atan2(cbet1 * sbet2 * comg12 | |
367 | - sbet1 * cbet2, | |
368 | cbet1 * somg12); | |
369 | ||
370 | excess = alp2 - alp1; | |
371 | ||
372 | } else { | |
373 | ||
374 | // Trapezoidal formula | |
375 | ||
376 | CT tan_lat1 = | |
377 | tan(geometry::get_as_radian<1>(p1) / 2.0); | |
378 | CT tan_lat2 = | |
379 | tan(geometry::get_as_radian<1>(p2) / 2.0); | |
380 | ||
381 | excess = CT(2.0) | |
382 | * atan(((tan_lat1 + tan_lat2) / (CT(1) + tan_lat1 * tan_lat2)) | |
383 | * tan((geometry::get_as_radian<0>(p2) | |
384 | - geometry::get_as_radian<0>(p1)) / 2)); | |
385 | } | |
386 | ||
387 | return excess; | |
388 | } | |
389 | ||
390 | struct return_type_ellipsoidal | |
391 | { | |
392 | return_type_ellipsoidal() | |
393 | : spherical_term(0), | |
394 | ellipsoidal_term(0) | |
395 | {} | |
396 | ||
397 | CT spherical_term; | |
398 | CT ellipsoidal_term; | |
399 | }; | |
400 | ||
401 | /* | |
402 | Compute the ellipsoidal correction of a geodesic (or shperical) segment | |
403 | */ | |
404 | template < | |
405 | template <typename, bool, bool, bool, bool, bool> class Inverse, | |
406 | //typename AzimuthStrategy, | |
407 | typename PointOfSegment, | |
408 | typename SpheroidConst | |
409 | > | |
410 | static inline return_type_ellipsoidal ellipsoidal(PointOfSegment const& p1, | |
411 | PointOfSegment const& p2, | |
412 | SpheroidConst spheroid_const) | |
413 | { | |
414 | return_type_ellipsoidal result; | |
415 | ||
416 | // Azimuth Approximation | |
417 | ||
418 | typedef Inverse<CT, false, true, true, false, false> inverse_type; | |
419 | typedef typename inverse_type::result_type inverse_result; | |
420 | ||
421 | inverse_result i_res = inverse_type::apply(get_as_radian<0>(p1), | |
422 | get_as_radian<1>(p1), | |
423 | get_as_radian<0>(p2), | |
424 | get_as_radian<1>(p2), | |
425 | spheroid_const.m_spheroid); | |
426 | ||
427 | CT alp1 = i_res.azimuth; | |
428 | CT alp2 = i_res.reverse_azimuth; | |
429 | ||
430 | // Constants | |
431 | ||
432 | CT const ep = spheroid_const.m_ep; | |
433 | CT const f = formula::flattening<CT>(spheroid_const.m_spheroid); | |
434 | CT const one_minus_f = CT(1) - f; | |
435 | std::size_t const series_order_plus_one = SeriesOrder + 1; | |
436 | std::size_t const series_order_plus_two = SeriesOrder + 2; | |
437 | ||
438 | // Basic trigonometric computations | |
439 | ||
440 | CT tan_bet1 = tan(get_as_radian<1>(p1)) * one_minus_f; | |
441 | CT tan_bet2 = tan(get_as_radian<1>(p2)) * one_minus_f; | |
442 | CT cos_bet1 = cos(atan(tan_bet1)); | |
443 | CT cos_bet2 = cos(atan(tan_bet2)); | |
444 | CT sin_bet1 = tan_bet1 * cos_bet1; | |
445 | CT sin_bet2 = tan_bet2 * cos_bet2; | |
446 | CT sin_alp1 = sin(alp1); | |
447 | CT cos_alp1 = cos(alp1); | |
448 | CT cos_alp2 = cos(alp2); | |
449 | CT sin_alp0 = sin_alp1 * cos_bet1; | |
450 | ||
451 | // Spherical term computation | |
452 | ||
453 | CT sin_omg1 = sin_alp0 * sin_bet1; | |
454 | CT cos_omg1 = cos_alp1 * cos_bet1; | |
455 | CT sin_omg2 = sin_alp0 * sin_bet2; | |
456 | CT cos_omg2 = cos_alp2 * cos_bet2; | |
457 | CT cos_omg12 = cos_omg1 * cos_omg2 + sin_omg1 * sin_omg2; | |
458 | CT excess; | |
459 | ||
460 | bool meridian = get<0>(p2) - get<0>(p1) == CT(0) | |
461 | || get<1>(p1) == CT(90) || get<1>(p1) == -CT(90) | |
462 | || get<1>(p2) == CT(90) || get<1>(p2) == -CT(90); | |
463 | ||
464 | if (!meridian && cos_omg12 > -CT(0.7) | |
465 | && sin_bet2 - sin_bet1 < CT(1.75)) // short segment | |
466 | { | |
467 | CT sin_omg12 = cos_omg1 * sin_omg2 - sin_omg1 * cos_omg2; | |
468 | normalize(sin_omg12, cos_omg12); | |
469 | ||
470 | CT cos_omg12p1 = CT(1) + cos_omg12; | |
471 | CT cos_bet1p1 = CT(1) + cos_bet1; | |
472 | CT cos_bet2p1 = CT(1) + cos_bet2; | |
473 | excess = CT(2) * atan2(sin_omg12 * (sin_bet1 * cos_bet2p1 + sin_bet2 * cos_bet1p1), | |
474 | cos_omg12p1 * (sin_bet1 * sin_bet2 + cos_bet1p1 * cos_bet2p1)); | |
475 | } | |
476 | else | |
477 | { | |
478 | /* | |
479 | CT sin_alp2 = sin(alp2); | |
480 | CT sin_alp12 = sin_alp2 * cos_alp1 - cos_alp2 * sin_alp1; | |
481 | CT cos_alp12 = cos_alp2 * cos_alp1 + sin_alp2 * sin_alp1; | |
482 | excess = atan2(sin_alp12, cos_alp12); | |
483 | */ | |
484 | excess = alp2 - alp1; | |
485 | } | |
486 | ||
487 | result.spherical_term = excess; | |
488 | ||
489 | // Ellipsoidal term computation (uses integral approximation) | |
490 | ||
491 | CT cos_alp0 = math::sqrt(CT(1) - math::sqr(sin_alp0)); | |
492 | CT cos_sig1 = cos_alp1 * cos_bet1; | |
493 | CT cos_sig2 = cos_alp2 * cos_bet2; | |
494 | CT sin_sig1 = sin_bet1; | |
495 | CT sin_sig2 = sin_bet2; | |
496 | ||
497 | normalize(sin_sig1, cos_sig1); | |
498 | normalize(sin_sig2, cos_sig2); | |
499 | ||
500 | CT coeffs[SeriesOrder + 1]; | |
501 | const std::size_t coeffs_var_size = (series_order_plus_two | |
502 | * series_order_plus_one) / 2; | |
503 | CT coeffs_var[coeffs_var_size]; | |
504 | ||
505 | if(ExpandEpsN){ // expand by eps and n | |
506 | ||
507 | CT k2 = math::sqr(ep * cos_alp0); | |
508 | CT sqrt_k2_plus_one = math::sqrt(CT(1) + k2); | |
509 | CT eps = (sqrt_k2_plus_one - CT(1)) / (sqrt_k2_plus_one + CT(1)); | |
510 | CT n = f / (CT(2) - f); | |
511 | ||
512 | // Generate and evaluate the polynomials on n | |
513 | // to get the series coefficients (that depend on eps) | |
514 | evaluate_coeffs_n(n, coeffs_var); | |
515 | ||
516 | // Generate and evaluate the polynomials on eps (i.e. var2 = eps) | |
517 | // to get the final series coefficients | |
518 | evaluate_coeffs_var2(eps, coeffs_var, coeffs); | |
519 | ||
520 | }else{ // expand by k2 and ep | |
521 | ||
522 | CT k2 = math::sqr(ep * cos_alp0); | |
523 | CT ep2 = math::sqr(ep); | |
524 | ||
525 | // Generate and evaluate the polynomials on ep2 | |
526 | evaluate_coeffs_ep(ep2, coeffs_var); | |
527 | ||
528 | // Generate and evaluate the polynomials on k2 (i.e. var2 = k2) | |
529 | evaluate_coeffs_var2(k2, coeffs_var, coeffs); | |
530 | ||
531 | } | |
532 | ||
533 | // Evaluate the trigonometric sum | |
534 | CT I12 = clenshaw_sum(cos_sig2, coeffs, coeffs + series_order_plus_one) | |
535 | - clenshaw_sum(cos_sig1, coeffs, coeffs + series_order_plus_one); | |
536 | ||
537 | // The part of the ellipsodal correction that depends on | |
538 | // point coordinates | |
539 | result.ellipsoidal_term = cos_alp0 * sin_alp0 * I12; | |
540 | ||
541 | return result; | |
542 | } | |
543 | ||
544 | // Check whenever a segment crosses the prime meridian | |
545 | // First normalize to [0,360) | |
546 | template <typename PointOfSegment> | |
547 | static inline bool crosses_prime_meridian(PointOfSegment const& p1, | |
548 | PointOfSegment const& p2) | |
549 | { | |
550 | CT const pi | |
551 | = geometry::math::pi<CT>(); | |
552 | CT const two_pi | |
553 | = geometry::math::two_pi<CT>(); | |
554 | ||
555 | CT p1_lon = get_as_radian<0>(p1) | |
556 | - ( floor( get_as_radian<0>(p1) / two_pi ) | |
557 | * two_pi ); | |
558 | CT p2_lon = get_as_radian<0>(p2) | |
559 | - ( floor( get_as_radian<0>(p2) / two_pi ) | |
560 | * two_pi ); | |
561 | ||
562 | CT max_lon = (std::max)(p1_lon, p2_lon); | |
563 | CT min_lon = (std::min)(p1_lon, p2_lon); | |
564 | ||
565 | return max_lon > pi && min_lon < pi && max_lon - min_lon > pi; | |
566 | } | |
567 | ||
568 | }; | |
569 | ||
570 | }}} // namespace boost::geometry::formula | |
571 | ||
572 | ||
573 | #endif // BOOST_GEOMETRY_FORMULAS_AREA_FORMULAS_HPP |