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_SPHERICAL_AREA_HPP
12 #define BOOST_GEOMETRY_STRATEGIES_SPHERICAL_AREA_HPP
13 
14 
15 #include <boost/geometry/formulas/area_formulas.hpp>
16 #include <boost/geometry/core/radius.hpp>
17 #include <boost/geometry/core/srs.hpp>
18 #include <boost/geometry/strategies/area.hpp>
19 
20 
21 namespace boost { namespace geometry
22 {
23 
24 namespace strategy { namespace area
25 {
26 
27 /*!
28 \brief Spherical area calculation
29 \ingroup strategies
30 \details Calculates area on the surface of a sphere using the trapezoidal rule
31 \tparam PointOfSegment \tparam_segment_point
32 \tparam CalculationType \tparam_calculation
33 
34 \qbk{
35 [heading See also]
36 [link geometry.reference.algorithms.area.area_2_with_strategy area (with strategy)]
37 }
38 */
39 template
40 <
41     typename PointOfSegment,
42     typename CalculationType = void
43 >
44 class spherical
45 {
46     // Enables special handling of long segments
47     static const bool LongSegment = false;
48 
49 typedef typename boost::mpl::if_c
50     <
51         boost::is_void<CalculationType>::type::value,
52         typename select_most_precise
53             <
54                 typename coordinate_type<PointOfSegment>::type,
55                 double
56             >::type,
57         CalculationType
58     >::type CT;
59 
60 protected :
61     struct excess_sum
62     {
63         CT m_sum;
64 
65         // Keep track if encircles some pole
66         size_t m_crosses_prime_meridian;
67 
excess_sumboost::geometry::strategy::area::spherical::excess_sum68         inline excess_sum()
69             : m_sum(0)
70             , m_crosses_prime_meridian(0)
71         {}
72         template <typename SphereType>
areaboost::geometry::strategy::area::spherical::excess_sum73         inline CT area(SphereType sphere) const
74         {
75             CT result;
76             CT radius = geometry::get_radius<0>(sphere);
77 
78             // Encircles pole
79             if(m_crosses_prime_meridian % 2 == 1)
80             {
81                 size_t times_crosses_prime_meridian
82                         = 1 + (m_crosses_prime_meridian / 2);
83 
84                 result = CT(2)
85                          * geometry::math::pi<CT>()
86                          * times_crosses_prime_meridian
87                          - geometry::math::abs(m_sum);
88 
89                 if(geometry::math::sign<CT>(m_sum) == 1)
90                 {
91                     result = - result;
92                 }
93 
94             } else {
95                 result =  m_sum;
96             }
97 
98             result *= radius * radius;
99 
100             return result;
101         }
102     };
103 
104 public :
105     typedef CT return_type;
106     typedef PointOfSegment segment_point_type;
107     typedef excess_sum state_type;
108     typedef geometry::srs::sphere<CT> sphere_type;
109 
110     // For backward compatibility reasons the radius is set to 1
spherical()111     inline spherical()
112         : m_sphere(1.0)
113     {}
114 
115     template <typename T>
spherical(geometry::srs::sphere<T> const & sphere)116     explicit inline spherical(geometry::srs::sphere<T> const& sphere)
117         : m_sphere(geometry::get_radius<0>(sphere))
118     {}
119 
spherical(CT const & radius)120     explicit inline spherical(CT const& radius)
121         : m_sphere(radius)
122     {}
123 
apply(PointOfSegment const & p1,PointOfSegment const & p2,excess_sum & state) const124     inline void apply(PointOfSegment const& p1,
125                       PointOfSegment const& p2,
126                       excess_sum& state) const
127     {
128         if (! geometry::math::equals(get<0>(p1), get<0>(p2)))
129         {
130 
131             state.m_sum += geometry::formula::area_formulas
132                              <CT>::template spherical<LongSegment>(p1, p2);
133 
134             // Keep track whenever a segment crosses the prime meridian
135             geometry::formula::area_formulas
136                          <CT>::crosses_prime_meridian(p1, p2, state);
137 
138         }
139     }
140 
result(excess_sum const & state) const141     inline return_type result(excess_sum const& state) const
142     {
143         return state.area(m_sphere);
144     }
145 
146 private :
147     /// srs Sphere
148     sphere_type m_sphere;
149 };
150 
151 #ifndef DOXYGEN_NO_STRATEGY_SPECIALIZATIONS
152 
153 namespace services
154 {
155 
156 
157 template <typename Point>
158 struct default_strategy<spherical_equatorial_tag, Point>
159 {
160     typedef strategy::area::spherical<Point> type;
161 };
162 
163 // Note: spherical polar coordinate system requires "get_as_radian_equatorial"
164 template <typename Point>
165 struct default_strategy<spherical_polar_tag, Point>
166 {
167     typedef strategy::area::spherical<Point> type;
168 };
169 
170 } // namespace services
171 
172 #endif // DOXYGEN_NO_STRATEGY_SPECIALIZATIONS
173 
174 
175 }} // namespace strategy::area
176 
177 
178 
179 
180 }} // namespace boost::geometry
181 
182 #endif // BOOST_GEOMETRY_STRATEGIES_SPHERICAL_AREA_HPP
183