ceph/src/boost/boost/geometry/strategies/geographic/area.hpp

   1 // Boost.Geometry (aka GGL, Generic Geometry Library)
   2
   3 // Copyright (c) 2016-2017 Oracle and/or its affiliates.
   4 // Contributed and/or modified by Vissarion Fisikopoulos, on behalf of Oracle
   5 // Contributed and/or modified by Adam Wulkiewicz, 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_STRATEGIES_GEOGRAPHIC_AREA_HPP
  12 #define BOOST_GEOMETRY_STRATEGIES_GEOGRAPHIC_AREA_HPP
  13
  14
  15 #include <boost/geometry/core/srs.hpp>
  16
  17 #include <boost/geometry/formulas/area_formulas.hpp>
  18 #include <boost/geometry/formulas/authalic_radius_sqr.hpp>
  19 #include <boost/geometry/formulas/eccentricity_sqr.hpp>
  20
  21 #include <boost/geometry/strategies/geographic/parameters.hpp>
  22
  23
  24 namespace boost { namespace geometry
  25 {
  26
  27 namespace strategy { namespace area
  28 {
  29
  30 /*!
  31 \brief Geographic area calculation
  32 \ingroup strategies
  33 \details Geographic area calculation by trapezoidal rule plus integral
  34          approximation that gives the ellipsoidal correction
  35 \tparam PointOfSegment \tparam_segment_point
  36 \tparam FormulaPolicy Formula used to calculate azimuths
  37 \tparam SeriesOrder The order of approximation of the geodesic integral
  38 \tparam Spheroid The spheroid model
  39 \tparam CalculationType \tparam_calculation
  40 \author See
  41 - Danielsen JS, The area under the geodesic. Surv Rev 30(232): 61–66, 1989
  42 - Charles F.F Karney, Algorithms for geodesics, 2011 https://arxiv.org/pdf/1109.4448.pdf
  43
  44 \qbk{
  45 [heading See also]
  46 [link geometry.reference.algorithms.area.area_2_with_strategy area (with strategy)]
  47 }
  48 */
  49 template
  50 <
  51     typename PointOfSegment,
  52     typename FormulaPolicy = strategy::andoyer,
  53     std::size_t SeriesOrder = strategy::default_order<FormulaPolicy>::value,
  54     typename Spheroid = srs::spheroid<double>,
  55     typename CalculationType = void
  56 >
  57 class geographic
  58 {
  59     // Switch between two kinds of approximation(series in eps and n v.s.series in k ^ 2 and e'^2)
  60     static const bool ExpandEpsN = true;
  61     // LongSegment Enables special handling of long segments
  62     static const bool LongSegment = false;
  63
  64     //Select default types in case they are not set
  65
  66     typedef typename boost::mpl::if_c
  67     <
  68         boost::is_void<CalculationType>::type::value,
  69         typename select_most_precise
  70             <
  71                 typename coordinate_type<PointOfSegment>::type,
  72                 double
  73             >::type,
  74         CalculationType
  75     >::type CT;
  76
  77 protected :
  78     struct spheroid_constants
  79     {
  80         Spheroid m_spheroid;
  81         CT const m_a2;  // squared equatorial radius
  82         CT const m_e2;  // squared eccentricity
  83         CT const m_ep2; // squared second eccentricity
  84         CT const m_ep;  // second eccentricity
  85         CT const m_c2;  // squared authalic radius
  86
  87         inline spheroid_constants(Spheroid const& spheroid)
  88             : m_spheroid(spheroid)
  89             , m_a2(math::sqr(get_radius<0>(spheroid)))
  90             , m_e2(formula::eccentricity_sqr<CT>(spheroid))
  91             , m_ep2(m_e2 / (CT(1.0) - m_e2))
  92             , m_ep(math::sqrt(m_ep2))
  93             , m_c2(formula_dispatch::authalic_radius_sqr
  94                     <
  95                         CT, Spheroid, srs_spheroid_tag
  96                     >::apply(m_a2, m_e2))
  97         {}
  98     };
  99
 100     struct area_sums
 101     {
 102         CT m_excess_sum;
 103         CT m_correction_sum;
 104
 105         // Keep track if encircles some pole
 106         std::size_t m_crosses_prime_meridian;
 107
 108         inline area_sums()
 109             : m_excess_sum(0)
 110             , m_correction_sum(0)
 111             , m_crosses_prime_meridian(0)
 112         {}
 113         inline CT area(spheroid_constants spheroid_const) const
 114         {
 115             CT result;
 116
 117             CT sum = spheroid_const.m_c2 * m_excess_sum
 118                    + spheroid_const.m_e2 * spheroid_const.m_a2 * m_correction_sum;
 119
 120             // If encircles some pole
 121             if (m_crosses_prime_meridian % 2 == 1)
 122             {
 123                 std::size_t times_crosses_prime_meridian
 124                         = 1 + (m_crosses_prime_meridian / 2);
 125
 126                 result = CT(2.0)
 127                          * geometry::math::pi<CT>()
 128                          * spheroid_const.m_c2
 129                          * CT(times_crosses_prime_meridian)
 130                          - geometry::math::abs(sum);
 131
 132                 if (geometry::math::sign<CT>(sum) == 1)
 133                 {
 134                     result = - result;
 135                 }
 136
 137             }
 138             else
 139             {
 140                 result = sum;
 141             }
 142
 143             return result;
 144         }
 145     };
 146
 147 public :
 148     typedef CT return_type;
 149     typedef PointOfSegment segment_point_type;
 150     typedef area_sums state_type;
 151
 152     explicit inline geographic(Spheroid const& spheroid = Spheroid())
 153         : m_spheroid_constants(spheroid)
 154     {}
 155
 156     inline void apply(PointOfSegment const& p1,
 157                       PointOfSegment const& p2,
 158                       area_sums& state) const
 159     {
 160
 161         if (! geometry::math::equals(get<0>(p1), get<0>(p2)))
 162         {
 163
 164             typedef geometry::formula::area_formulas
 165                 <
 166                     CT, SeriesOrder, ExpandEpsN
 167                 > area_formulas;
 168
 169             typename area_formulas::return_type_ellipsoidal result =
 170                      area_formulas::template ellipsoidal<FormulaPolicy::template inverse>
 171                                              (p1, p2, m_spheroid_constants);
 172
 173             state.m_excess_sum += result.spherical_term;
 174             state.m_correction_sum += result.ellipsoidal_term;
 175
 176             // Keep track whenever a segment crosses the prime meridian
 177             if (area_formulas::crosses_prime_meridian(p1, p2))
 178             {
 179                 state.m_crosses_prime_meridian++;
 180             }
 181         }
 182     }
 183
 184     inline return_type result(area_sums const& state) const
 185     {
 186         return state.area(m_spheroid_constants);
 187     }
 188
 189 private:
 190     spheroid_constants m_spheroid_constants;
 191
 192 };
 193
 194 #ifndef DOXYGEN_NO_STRATEGY_SPECIALIZATIONS
 195
 196 namespace services
 197 {
 198
 199
 200 template <typename Point>
 201 struct default_strategy<geographic_tag, Point>
 202 {
 203     typedef strategy::area::geographic<Point> type;
 204 };
 205
 206 #endif // DOXYGEN_NO_STRATEGY_SPECIALIZATIONS
 207
 208 }
 209
 210 }} // namespace strategy::area
 211
 212
 213
 214
 215 }} // namespace boost::geometry
 216
 217 #endif // BOOST_GEOMETRY_STRATEGIES_GEOGRAPHIC_AREA_HPP