ceph/src/boost/boost/geometry/formulas/area_formulas.hpp

   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