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

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