geometry/formulas/area_formulas.hpp

0001 // Boost.Geometry
0002
0003 // Copyright (c) 2023-2024 Adam Wulkiewicz, Lodz, Poland.
0004
0005 // Copyright (c) 2015-2022 Oracle and/or its affiliates.
0006 // Contributed and/or modified by Vissarion Fysikopoulos, on behalf of Oracle
0007 // Contributed and/or modified by Adam Wulkiewicz, on behalf of Oracle
0008
0009 // Use, modification and distribution is subject to the Boost Software License,
0010 // Version 1.0. (See accompanying file LICENSE_1_0.txt or copy at
0011 // http://www.boost.org/LICENSE_1_0.txt)
0012
0013 #ifndef BOOST_GEOMETRY_FORMULAS_AREA_FORMULAS_HPP
0014 #define BOOST_GEOMETRY_FORMULAS_AREA_FORMULAS_HPP
0015
0016 #include <boost/geometry/core/radian_access.hpp>
0017 #include <boost/geometry/formulas/flattening.hpp>
0018 #include <boost/geometry/formulas/mean_radius.hpp>
0019 #include <boost/geometry/formulas/karney_inverse.hpp>
0020 #include <boost/geometry/util/constexpr.hpp>
0021 #include <boost/geometry/util/math.hpp>
0022 #include <boost/math/special_functions/hypot.hpp>
0023
0024 namespace boost { namespace geometry { namespace formula
0025 {
0026
0027 /*!
0028 \brief Formulas for computing spherical and ellipsoidal polygon area.
0029  The current class computes the area of the trapezoid defined by a segment
0030  the two meridians passing by the endpoints and the equator.
0031 \author See
0032 - Danielsen JS, The area under the geodesic. Surv Rev 30(232):
0033 61–66, 1989
0034 - Charles F.F Karney, Algorithms for geodesics, 2011
0035 https://arxiv.org/pdf/1109.4448.pdf
0036 */
0037
0038 template
0039 <
0040     typename CT,
0041     std::size_t SeriesOrder = 2,
0042     bool ExpandEpsN = true
0043 >
0044 class area_formulas
0045 {
0046
0047 public:
0048
0049     //TODO: move the following to a more general space to be used by other
0050     //      classes as well
0051     /*
0052         Evaluate the polynomial in x using Horner's method.
0053     */
0054     template <typename NT, typename IteratorType>
0055     static inline NT horner_evaluate(NT const& x,
0056                                      IteratorType begin,
0057                                      IteratorType end)
0058     {
0059         NT result(0);
0060         IteratorType it = end;
0061         do
0062         {
0063             result = result * x + *--it;
0064         }
0065         while (it != begin);
0066         return result;
0067     }
0068
0069     /*
0070         Clenshaw algorithm for summing trigonometric series
0071         https://en.wikipedia.org/wiki/Clenshaw_algorithm
0072     */
0073     template <typename NT, typename IteratorType>
0074     static inline NT clenshaw_sum(NT const& cosx,
0075                                   IteratorType begin,
0076                                   IteratorType end)
0077     {
0078         IteratorType it = end;
0079         bool odd = true;
0080         CT b_k, b_k1(0), b_k2(0);
0081         do
0082         {
0083             CT c_k = odd ? *--it : NT(0);
0084             b_k = c_k + NT(2) * cosx * b_k1 - b_k2;
0085             b_k2 = b_k1;
0086             b_k1 = b_k;
0087             odd = !odd;
0088         }
0089         while (it != begin);
0090
0091         return *begin + b_k1 * cosx - b_k2;
0092     }
0093
0094     template<typename T>
0095     static inline void normalize(T& x, T& y)
0096     {
0097         T h = boost::math::hypot(x, y);
0098         x /= h;
0099         y /= h;
0100     }
0101
0102     /*
0103      Generate and evaluate the series expansion of the following integral
0104
0105         I4 = -integrate( (t(ep2) - t(k2*sin(sigma1)^2)) / (ep2 - k2*sin(sigma1)^2)
0106            * sin(sigma1)/2, sigma1, pi/2, sigma )
0107      where
0108
0109         t(x) = sqrt(1+1/x)*asinh(sqrt(x)) + x
0110
0111      valid for ep2 and k2 small.  We substitute k2 = 4 * eps / (1 - eps)^2
0112      and ep2 = 4 * n / (1 - n)^2 and expand in eps and n.
0113
0114      The resulting sum of the series is of the form
0115
0116         sum(C4[l] * cos((2*l+1)*sigma), l, 0, maxpow-1) )
0117
0118      The above expansion is performed in Computer Algebra System Maxima.
0119      The C++ code (that yields the function evaluate_coeffs_n below) is generated
0120      by the following Maxima script and is based on script:
0121      http://geographiclib.sourceforge.net/html/geod.mac
0122
0123         // Maxima script begin
0124         taylordepth:5$
0125         ataylor(expr,var,ord):=expand(ratdisrep(taylor(expr,var,0,ord)))$
0126         jtaylor(expr,var1,var2,ord):=block([zz],expand(subst([zz=1],
0127         ratdisrep(taylor(subst([var1=zz*var1,var2=zz*var2],expr),zz,0,ord)))))$
0128
0129         compute(maxpow):=block([int,t,intexp,area, x,ep2,k2],
0130         maxpow:maxpow-1,
0131         t : sqrt(1+1/x) * asinh(sqrt(x)) + x,
0132         int:-(tf(ep2) - tf(k2*sin(sigma)^2)) / (ep2 - k2*sin(sigma)^2)
0133         * sin(sigma)/2,
0134         int:subst([tf(ep2)=subst([x=ep2],t),
0135         tf(k2*sin(sigma)^2)=subst([x=k2*sin(sigma)^2],t)],
0136         int),
0137         int:subst([abs(sin(sigma))=sin(sigma)],int),
0138         int:subst([k2=4*eps/(1-eps)^2,ep2=4*n/(1-n)^2],int),
0139         intexp:jtaylor(int,n,eps,maxpow),
0140         area:trigreduce(integrate(intexp,sigma)),
0141         area:expand(area-subst(sigma=%pi/2,area)),
0142         for i:0 thru maxpow do C4[i]:coeff(area,cos((2*i+1)*sigma)),
0143         if expand(area-sum(C4[i]*cos((2*i+1)*sigma),i,0,maxpow)) # 0
0144         then error("left over terms in I4"),
0145         'done)$
0146
0147         printcode(maxpow):=
0148         block([tab2:"    ",tab3:"      "],
0149         print(" switch (SeriesOrder) {"),
0150         for nn:1 thru maxpow do block([c],
0151         print(concat(tab2,"case ",string(nn-1),":")),
0152         c:0,
0153         for m:0 thru nn-1 do block(
0154           [q:jtaylor(subst([n=n],C4[m]),n,eps,nn-1),
0155           linel:1200],
0156           for j:m thru nn-1 do (
0157             print(concat(tab3,"coeffs_n[",c,"] = ",
0158                 string(horner(coeff(q,eps,j))),";")),
0159             c:c+1)
0160         ),
0161         print(concat(tab3,"break;"))),
0162         print("    }"),
0163         'done)$
0164
0165         maxpow:6$
0166         compute(maxpow)$
0167         printcode(maxpow)$
0168         // Maxima script end
0169
0170      In the resulting code we should replace each number x by CT(x)
0171      e.g. using the following scirpt:
0172        sed -e 's/[0-9]\+/CT(&)/g; s/\[CT(/\[/g; s/)\]/\]/g;
0173                s/case\sCT(/case /g; s/):/:/g'
0174     */
0175
0176     static inline void evaluate_coeffs_n(CT const& n, CT coeffs_n[])
0177     {
0178         switch (SeriesOrder) {
0179         case 0:
0180             coeffs_n[0] = CT(2)/CT(3);
0181             break;
0182         case 1:
0183             coeffs_n[0] = (CT(10)-CT(4)*n)/CT(15);
0184             coeffs_n[1] = -CT(1)/CT(5);
0185             coeffs_n[2] = CT(1)/CT(45);
0186             break;
0187         case 2:
0188             coeffs_n[0] = (n*(CT(8)*n-CT(28))+CT(70))/CT(105);
0189             coeffs_n[1] = (CT(16)*n-CT(7))/CT(35);
0190             coeffs_n[2] = -CT(2)/CT(105);
0191             coeffs_n[3] = (CT(7)-CT(16)*n)/CT(315);
0192             coeffs_n[4] = -CT(2)/CT(105);
0193             coeffs_n[5] = CT(4)/CT(525);
0194             break;
0195         case 3:
0196             coeffs_n[0] = (n*(n*(CT(4)*n+CT(24))-CT(84))+CT(210))/CT(315);
0197             coeffs_n[1] = ((CT(48)-CT(32)*n)*n-CT(21))/CT(105);
0198             coeffs_n[2] = (-CT(32)*n-CT(6))/CT(315);
0199             coeffs_n[3] = CT(11)/CT(315);
0200             coeffs_n[4] = (n*(CT(32)*n-CT(48))+CT(21))/CT(945);
0201             coeffs_n[5] = (CT(64)*n-CT(18))/CT(945);
0202             coeffs_n[6] = -CT(1)/CT(105);
0203             coeffs_n[7] = (CT(12)-CT(32)*n)/CT(1575);
0204             coeffs_n[8] = -CT(8)/CT(1575);
0205             coeffs_n[9] = CT(8)/CT(2205);
0206             break;
0207         case 4:
0208             coeffs_n[0] = (n*(n*(n*(CT(16)*n+CT(44))+CT(264))-CT(924))+CT(2310))/CT(3465);
0209             coeffs_n[1] = (n*(n*(CT(48)*n-CT(352))+CT(528))-CT(231))/CT(1155);
0210             coeffs_n[2] = (n*(CT(1088)*n-CT(352))-CT(66))/CT(3465);
0211             coeffs_n[3] = (CT(121)-CT(368)*n)/CT(3465);
0212             coeffs_n[4] = CT(4)/CT(1155);
0213             coeffs_n[5] = (n*((CT(352)-CT(48)*n)*n-CT(528))+CT(231))/CT(10395);
0214             coeffs_n[6] = ((CT(704)-CT(896)*n)*n-CT(198))/CT(10395);
0215             coeffs_n[7] = (CT(80)*n-CT(99))/CT(10395);
0216             coeffs_n[8] = CT(4)/CT(1155);
0217             coeffs_n[9] = (n*(CT(320)*n-CT(352))+CT(132))/CT(17325);
0218             coeffs_n[10] = (CT(384)*n-CT(88))/CT(17325);
0219             coeffs_n[11] = -CT(8)/CT(1925);
0220             coeffs_n[12] = (CT(88)-CT(256)*n)/CT(24255);
0221             coeffs_n[13] = -CT(16)/CT(8085);
0222             coeffs_n[14] = CT(64)/CT(31185);
0223             break;
0224         case 5:
0225             coeffs_n[0] = (n*(n*(n*(n*(CT(100)*n+CT(208))+CT(572))+CT(3432))-CT(12012))+CT(30030))
0226                           /CT(45045);
0227             coeffs_n[1] = (n*(n*(n*(CT(64)*n+CT(624))-CT(4576))+CT(6864))-CT(3003))/CT(15015);
0228             coeffs_n[2] = (n*((CT(14144)-CT(10656)*n)*n-CT(4576))-CT(858))/CT(45045);
0229             coeffs_n[3] = ((-CT(224)*n-CT(4784))*n+CT(1573))/CT(45045);
0230             coeffs_n[4] = (CT(1088)*n+CT(156))/CT(45045);
0231             coeffs_n[5] = CT(97)/CT(15015);
0232             coeffs_n[6] = (n*(n*((-CT(64)*n-CT(624))*n+CT(4576))-CT(6864))+CT(3003))/CT(135135);
0233             coeffs_n[7] = (n*(n*(CT(5952)*n-CT(11648))+CT(9152))-CT(2574))/CT(135135);
0234             coeffs_n[8] = (n*(CT(5792)*n+CT(1040))-CT(1287))/CT(135135);
0235             coeffs_n[9] = (CT(468)-CT(2944)*n)/CT(135135);
0236             coeffs_n[10] = CT(1)/CT(9009);
0237             coeffs_n[11] = (n*((CT(4160)-CT(1440)*n)*n-CT(4576))+CT(1716))/CT(225225);
0238             coeffs_n[12] = ((CT(4992)-CT(8448)*n)*n-CT(1144))/CT(225225);
0239             coeffs_n[13] = (CT(1856)*n-CT(936))/CT(225225);
0240             coeffs_n[14] = CT(8)/CT(10725);
0241             coeffs_n[15] = (n*(CT(3584)*n-CT(3328))+CT(1144))/CT(315315);
0242             coeffs_n[16] = (CT(1024)*n-CT(208))/CT(105105);
0243             coeffs_n[17] = -CT(136)/CT(63063);
0244             coeffs_n[18] = (CT(832)-CT(2560)*n)/CT(405405);
0245             coeffs_n[19] = -CT(128)/CT(135135);
0246             coeffs_n[20] = CT(128)/CT(99099);
0247             break;
0248         }
0249     }
0250
0251     /*
0252        Expand in k2 and ep2.
0253     */
0254     static inline void evaluate_coeffs_ep(CT const& ep, CT coeffs_n[])
0255     {
0256         switch (SeriesOrder) {
0257         case 0:
0258             coeffs_n[0] = CT(2)/CT(3);
0259             break;
0260         case 1:
0261             coeffs_n[0] = (CT(10)-ep)/CT(15);
0262             coeffs_n[1] = -CT(1)/CT(20);
0263             coeffs_n[2] = CT(1)/CT(180);
0264             break;
0265         case 2:
0266             coeffs_n[0] = (ep*(CT(4)*ep-CT(7))+CT(70))/CT(105);
0267             coeffs_n[1] = (CT(4)*ep-CT(7))/CT(140);
0268             coeffs_n[2] = CT(1)/CT(42);
0269             coeffs_n[3] = (CT(7)-CT(4)*ep)/CT(1260);
0270             coeffs_n[4] = -CT(1)/CT(252);
0271             coeffs_n[5] = CT(1)/CT(2100);
0272             break;
0273         case 3:
0274             coeffs_n[0] = (ep*((CT(12)-CT(8)*ep)*ep-CT(21))+CT(210))/CT(315);
0275             coeffs_n[1] = ((CT(12)-CT(8)*ep)*ep-CT(21))/CT(420);
0276             coeffs_n[2] = (CT(3)-CT(2)*ep)/CT(126);
0277             coeffs_n[3] = -CT(1)/CT(72);
0278             coeffs_n[4] = (ep*(CT(8)*ep-CT(12))+CT(21))/CT(3780);
0279             coeffs_n[5] = (CT(2)*ep-CT(3))/CT(756);
0280             coeffs_n[6] = CT(1)/CT(360);
0281             coeffs_n[7] = (CT(3)-CT(2)*ep)/CT(6300);
0282             coeffs_n[8] = -CT(1)/CT(1800);
0283             coeffs_n[9] = CT(1)/CT(17640);
0284             break;
0285         case 4:
0286             coeffs_n[0] = (ep*(ep*(ep*(CT(64)*ep-CT(88))+CT(132))-CT(231))+CT(2310))/CT(3465);
0287             coeffs_n[1] = (ep*(ep*(CT(64)*ep-CT(88))+CT(132))-CT(231))/CT(4620);
0288             coeffs_n[2] = (ep*(CT(16)*ep-CT(22))+CT(33))/CT(1386);
0289             coeffs_n[3] = (CT(8)*ep-CT(11))/CT(792);
0290             coeffs_n[4] = CT(1)/CT(110);
0291             coeffs_n[5] = (ep*((CT(88)-CT(64)*ep)*ep-CT(132))+CT(231))/CT(41580);
0292             coeffs_n[6] = ((CT(22)-CT(16)*ep)*ep-CT(33))/CT(8316);
0293             coeffs_n[7] = (CT(11)-CT(8)*ep)/CT(3960);
0294             coeffs_n[8] = -CT(1)/CT(495);
0295             coeffs_n[9] = (ep*(CT(16)*ep-CT(22))+CT(33))/CT(69300);
0296             coeffs_n[10] = (CT(8)*ep-CT(11))/CT(19800);
0297             coeffs_n[11] = CT(1)/CT(1925);
0298             coeffs_n[12] = (CT(11)-CT(8)*ep)/CT(194040);
0299             coeffs_n[13] = -CT(1)/CT(10780);
0300             coeffs_n[14] = CT(1)/CT(124740);
0301             break;
0302         case 5:
0303             coeffs_n[0] = (ep*(ep*(ep*((CT(832)-CT(640)*ep)*ep-CT(1144))+CT(1716))-CT(3003))+CT(30030))/CT(45045);
0304             coeffs_n[1] = (ep*(ep*((CT(832)-CT(640)*ep)*ep-CT(1144))+CT(1716))-CT(3003))/CT(60060);
0305             coeffs_n[2] = (ep*((CT(208)-CT(160)*ep)*ep-CT(286))+CT(429))/CT(18018);
0306             coeffs_n[3] = ((CT(104)-CT(80)*ep)*ep-CT(143))/CT(10296);
0307             coeffs_n[4] = (CT(13)-CT(10)*ep)/CT(1430);
0308             coeffs_n[5] = -CT(1)/CT(156);
0309             coeffs_n[6] = (ep*(ep*(ep*(CT(640)*ep-CT(832))+CT(1144))-CT(1716))+CT(3003))/CT(540540);
0310             coeffs_n[7] = (ep*(ep*(CT(160)*ep-CT(208))+CT(286))-CT(429))/CT(108108);
0311             coeffs_n[8] = (ep*(CT(80)*ep-CT(104))+CT(143))/CT(51480);
0312             coeffs_n[9] = (CT(10)*ep-CT(13))/CT(6435);
0313             coeffs_n[10] = CT(5)/CT(3276);
0314             coeffs_n[11] = (ep*((CT(208)-CT(160)*ep)*ep-CT(286))+CT(429))/CT(900900);
0315             coeffs_n[12] = ((CT(104)-CT(80)*ep)*ep-CT(143))/CT(257400);
0316             coeffs_n[13] = (CT(13)-CT(10)*ep)/CT(25025);
0317             coeffs_n[14] = -CT(1)/CT(2184);
0318             coeffs_n[15] = (ep*(CT(80)*ep-CT(104))+CT(143))/CT(2522520);
0319             coeffs_n[16] = (CT(10)*ep-CT(13))/CT(140140);
0320             coeffs_n[17] = CT(5)/CT(45864);
0321             coeffs_n[18] = (CT(13)-CT(10)*ep)/CT(1621620);
0322             coeffs_n[19] = -CT(1)/CT(58968);
0323             coeffs_n[20] = CT(1)/CT(792792);
0324             break;
0325         }
0326     }
0327
0328     /*
0329         Given the set of coefficients coeffs1[] evaluate on var2 and return
0330         the set of coefficients coeffs2[]
0331     */
0332     template <typename CoeffsType>
0333     static inline void evaluate_coeffs_var2(CT const& var2,
0334                                             CoeffsType const coeffs1[],
0335                                             CT coeffs2[])
0336     {
0337         std::size_t begin(0), end(0);
0338         for(std::size_t i = 0; i <= SeriesOrder; i++)
0339         {
0340             end = begin + SeriesOrder + 1 - i;
0341             coeffs2[i] = ((i==0) ? CT(1) : math::pow(var2, int(i)))
0342                         * horner_evaluate(var2, coeffs1 + begin, coeffs1 + end);
0343             begin = end;
0344         }
0345     }
0346
0347     static inline CT trapezoidal_formula(CT lat1r, CT lat2r, CT lon21r)
0348     {
0349         CT const c1 = CT(1);
0350         CT const c2 = CT(2);
0351         CT const tan_lat1 = tan(lat1r / c2);
0352         CT const tan_lat2 = tan(lat2r / c2);
0353
0354         return c2 * atan(((tan_lat1 + tan_lat2) / (c1 + tan_lat1 * tan_lat2))* tan(lon21r / c2));
0355     }
0356
0357     /*
0358         Compute the spherical excess of a geodesic (or shperical) segment
0359     */
0360     template
0361     <
0362         bool LongSegment,
0363         typename PointOfSegment
0364     >
0365     static inline CT spherical(PointOfSegment const& p1,
0366                                PointOfSegment const& p2)
0367     {
0368         CT const pi = math::pi<CT>();
0369
0370         CT excess;
0371
0372         CT const lon1r = get_as_radian<0>(p1);
0373         CT const lat1r = get_as_radian<1>(p1);
0374         CT const lon2r = get_as_radian<0>(p2);
0375         CT const lat2r = get_as_radian<1>(p2);
0376
0377         CT lon12r = lon2r - lon1r;
0378         math::normalize_longitude<radian, CT>(lon12r);
0379
0380         if (lon12r == pi || lon12r == -pi)
0381         {
0382             return pi;
0383         }
0384
0385         if BOOST_GEOMETRY_CONSTEXPR (LongSegment)
0386         {
0387             if (lat1r != lat2r) // not for segments parallel to equator
0388             {
0389                 CT const cbet1 = cos(lat1r);
0390                 CT const sbet1 = sin(lat1r);
0391                 CT const cbet2 = cos(lat2r);
0392                 CT const sbet2 = sin(lat2r);
0393
0394                 CT const omg12 = lon2r - lon1r;
0395                 CT const comg12 = cos(omg12);
0396                 CT const somg12 = sin(omg12);
0397
0398                 CT const cbet1_sbet2 = cbet1 * sbet2;
0399                 CT const sbet1_cbet2 = sbet1 * cbet2;
0400                 CT const alp1 = atan2(cbet1_sbet2 - sbet1_cbet2 * comg12, cbet2 * somg12);
0401                 CT const alp2 = atan2(cbet1_sbet2 * comg12 - sbet1_cbet2, cbet1 * somg12);
0402
0403                 excess = alp2 - alp1;
0404             }
0405             else
0406             {
0407                 excess = trapezoidal_formula(lat1r, lat2r, lon12r);
0408             }
0409         }
0410         else
0411         {
0412             excess = trapezoidal_formula(lat1r, lat2r, lon12r);
0413         }
0414
0415         return excess;
0416     }
0417
0418     struct return_type_ellipsoidal
0419     {
0420         return_type_ellipsoidal()
0421             :   spherical_term(0),
0422                 ellipsoidal_term(0)
0423         {}
0424
0425         CT spherical_term;
0426         CT ellipsoidal_term;
0427     };
0428
0429     /*
0430         Compute the ellipsoidal correction of a geodesic (or shperical) segment
0431     */
0432     template
0433     <
0434         template <typename, bool, bool, bool, bool, bool> class Inverse,
0435         typename PointOfSegment,
0436         typename SpheroidConst
0437     >
0438     static inline auto ellipsoidal(PointOfSegment const& p1,
0439                                    PointOfSegment const& p2,
0440                                    SpheroidConst const& spheroid_const)
0441     {
0442         return_type_ellipsoidal result;
0443
0444         CT const lon1r = get_as_radian<0>(p1);
0445         CT const lat1r = get_as_radian<1>(p1);
0446         CT const lon2r = get_as_radian<0>(p2);
0447         CT const lat2r = get_as_radian<1>(p2);
0448
0449         // Azimuth Approximation
0450
0451         using inverse_type = Inverse<CT, true, true, true, false, false>;
0452         auto i_res = inverse_type::apply(lon1r, lat1r, lon2r, lat2r, spheroid_const.m_spheroid);
0453
0454         CT const alp1 = i_res.azimuth;
0455         CT const alp2 = i_res.reverse_azimuth;
0456
0457         // Constants
0458
0459         CT const c0 = CT(0);
0460         CT const c1 = CT(1);
0461         CT const c2 = CT(2);
0462         CT const pi = math::pi<CT>();
0463         CT const half_pi = pi / c2;
0464         CT const ep = spheroid_const.m_ep;
0465         CT const one_minus_f = c1 - spheroid_const.m_f;
0466
0467         // Basic trigonometric computations
0468         // the compiler could optimize here using sincos function
0469         // TODO: optimization: those quantities are already computed in inverse formula
0470         // at least in some inverse formulas, so do not compute them again here
0471         /*
0472         CT sin_bet1 = sin(lat1r);
0473         CT cos_bet1 = cos(lat1r);
0474         CT sin_bet2 = sin(lat2r);
0475         CT cos_bet2 = cos(lat2r);
0476
0477         sin_bet1 *= one_minus_f;
0478         sin_bet2 *= one_minus_f;
0479         normalize(sin_bet1, cos_bet1);
0480         normalize(sin_bet2, cos_bet2);
0481         */
0482
0483         CT const tan_bet1 = tan(lat1r) * one_minus_f;
0484         CT const tan_bet2 = tan(lat2r) * one_minus_f;
0485         CT const cos_bet1 = cos(atan(tan_bet1));
0486         CT const cos_bet2 = cos(atan(tan_bet2));
0487         CT const sin_bet1 = tan_bet1 * cos_bet1;
0488         CT const sin_bet2 = tan_bet2 * cos_bet2;
0489
0490         CT const sin_alp1 = sin(alp1);
0491         CT const cos_alp1 = cos(alp1);
0492         CT const cos_alp2 = cos(alp2);
0493         CT const sin_alp0 = sin_alp1 * cos_bet1;
0494
0495         // Spherical term computation
0496
0497         CT excess;
0498
0499         CT lon12r = lon2r - lon1r;
0500         math::normalize_longitude<radian, CT>(lon12r);
0501
0502         // Comparing with "==" works with all test cases here, but could potential create numerical issues
0503         if (lon12r == pi || lon12r == -pi)
0504         {
0505             result.spherical_term = pi;
0506         }
0507         else
0508         {
0509             bool const meridian = lon12r == c0
0510                 || lat1r == half_pi || lat1r == -half_pi
0511                 || lat2r == half_pi || lat2r == -half_pi;
0512
0513             if (!meridian && (i_res.distance)
0514                 < mean_radius<CT>(spheroid_const.m_spheroid) / CT(638))  // short segment
0515             {
0516                 excess = trapezoidal_formula(lat1r, lat2r, lon12r);
0517             }
0518             else
0519             {
0520                 /* in some cases this formula gives more accurate results
0521                 CT sin_omg12 =  cos_omg1 * sin_omg2 - sin_omg1 * cos_omg2;
0522                 normalize(sin_omg12, cos_omg12);
0523
0524                 CT cos_omg12p1 = CT(1) + cos_omg12;
0525                 CT cos_bet1p1 = CT(1) + cos_bet1;
0526                 CT cos_bet2p1 = CT(1) + cos_bet2;
0527                 excess = CT(2) * atan2(sin_omg12 * (sin_bet1 * cos_bet2p1 + sin_bet2 * cos_bet1p1),
0528                     cos_omg12p1 * (sin_bet1 * sin_bet2 + cos_bet1p1 * cos_bet2p1));
0529                 */
0530
0531                 excess = alp2 - alp1;
0532             }
0533
0534             result.spherical_term = excess;
0535         }
0536
0537         // Ellipsoidal term computation (uses integral approximation)
0538
0539         CT const cos_alp0 = math::sqrt(c1 - math::sqr(sin_alp0));
0540         //CT const cos_alp0 = hypot(cos_alp1, sin_alp1 * sin_bet1);
0541         CT cos_sig1 = cos_alp1 * cos_bet1;
0542         CT cos_sig2 = cos_alp2 * cos_bet2;
0543         CT sin_sig1 = sin_bet1;
0544         CT sin_sig2 = sin_bet2;
0545
0546         normalize(sin_sig1, cos_sig1);
0547         normalize(sin_sig2, cos_sig2);
0548
0549         CT coeffs[SeriesOrder + 1];
0550
0551         if (ExpandEpsN) // expand by eps and n
0552         {
0553             CT const k2 = math::sqr(ep * cos_alp0);
0554             CT const sqrt_k2_plus_one = math::sqrt(c1 + k2);
0555             CT const eps = (sqrt_k2_plus_one - c1) / (sqrt_k2_plus_one + c1);
0556
0557             // Generate and evaluate the polynomials on eps (i.e. var2 = eps)
0558             // to get the final series coefficients
0559             evaluate_coeffs_var2(eps, spheroid_const.m_coeffs_var, coeffs);
0560         }
0561         else
0562         { // expand by k2 and ep
0563
0564             CT const k2 = math::sqr(ep * cos_alp0);
0565             CT const ep2 = math::sqr(ep);
0566
0567             CT coeffs_var[((SeriesOrder+2)*(SeriesOrder+1))/2];
0568
0569             // Generate and evaluate the polynomials on ep2
0570             evaluate_coeffs_ep(ep2, coeffs_var);
0571
0572             // Generate and evaluate the polynomials on k2 (i.e. var2 = k2)
0573             evaluate_coeffs_var2(k2, coeffs_var, coeffs);
0574         }
0575
0576         // Evaluate the trigonometric sum
0577         constexpr auto series_order_plus_one = SeriesOrder + 1;
0578         CT const I12 = clenshaw_sum(cos_sig2, coeffs, coeffs + series_order_plus_one)
0579             - clenshaw_sum(cos_sig1, coeffs, coeffs + series_order_plus_one);
0580
0581         // The part of the ellipsodal correction that depends on
0582         // point coordinates
0583         result.ellipsoidal_term = cos_alp0 * sin_alp0 * I12;
0584
0585         return result;
0586     }
0587
0588     // Check whenever a segment crosses the prime meridian
0589     // First normalize to [0,360)
0590     template <typename PointOfSegment>
0591     static inline bool crosses_prime_meridian(PointOfSegment const& p1,
0592                                               PointOfSegment const& p2)
0593     {
0594         CT const pi = geometry::math::pi<CT>();
0595         CT const two_pi = geometry::math::two_pi<CT>();
0596
0597         CT const lon1r = get_as_radian<0>(p1);
0598         CT const lon2r = get_as_radian<0>(p2);
0599
0600         CT lon12 = lon2r - lon1r;
0601         math::normalize_longitude<radian, CT>(lon12);
0602
0603         // Comparing with "==" works with all test cases here, but could potential create numerical issues
0604         if (lon12 == pi || lon12 == -pi)
0605         {
0606             return true;
0607         }
0608
0609         CT const p1_lon = lon1r - ( std::floor( lon1r / two_pi ) * two_pi );
0610         CT const p2_lon = lon2r - ( std::floor( lon2r / two_pi ) * two_pi );
0611
0612         CT const max_lon = (std::max)(p1_lon, p2_lon);
0613         CT const min_lon = (std::min)(p1_lon, p2_lon);
0614
0615         return max_lon > pi && min_lon < pi && max_lon - min_lon > pi;
0616     }
0617
0618 };
0619
0620 }}} // namespace boost::geometry::formula
0621
0622
0623 #endif // BOOST_GEOMETRY_FORMULAS_AREA_FORMULAS_HPP