Back to home page

EIC code displayed by LXR
 
 

    

File indexing completed on 2025-01-18 09:35:24 

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

This page was automatically generated by the 2.3.7 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!