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