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