1 // Boost.Geometry
2 
3 // Copyright (c) 2018 Oracle and/or its affiliates.
4 
5 // Contributed and/or modified by Vissarion Fysikopoulos, on behalf of Oracle
6 // Contributed and/or modified by Adam Wulkiewicz, on behalf of Oracle
7 
8 // Use, modification and distribution is subject to the Boost Software License,
9 // Version 1.0. (See accompanying file LICENSE_1_0.txt or copy at
10 // http://www.boost.org/LICENSE_1_0.txt)
11 
12 #ifndef BOOST_GEOMETRY_FORMULAS_QUARTER_MERIDIAN_HPP
13 #define BOOST_GEOMETRY_FORMULAS_QUARTER_MERIDIAN_HPP
14 
15 #include <boost/geometry/algorithms/not_implemented.hpp>
16 
17 #include <boost/geometry/core/radius.hpp>
18 #include <boost/geometry/core/tag.hpp>
19 #include <boost/geometry/core/tags.hpp>
20 
21 #include <boost/geometry/formulas/flattening.hpp>
22 
23 #include <boost/geometry/util/math.hpp>
24 
25 namespace boost { namespace geometry
26 {
27 
28 #ifndef DOXYGEN_NO_DISPATCH
29 namespace formula_dispatch
30 {
31 
32 template <typename ResultType, typename Geometry, typename Tag = typename tag<Geometry>::type>
33 struct quarter_meridian
34     : not_implemented<Tag>
35 {};
36 
37 template <typename ResultType, typename Geometry>
38 struct quarter_meridian<ResultType, Geometry, srs_spheroid_tag>
39 {
40     //https://en.wikipedia.org/wiki/Meridian_arc#Generalized_series
41     //http://www.wolframalpha.com/input/?i=(sum(((2*j-3)!!%2F(2*j)!!)%5E2*n%5E(2*j),j,0,8))
applyboost::geometry::formula_dispatch::quarter_meridian42     static inline ResultType apply(Geometry const& geometry)
43     {
44         //order 8 expansion
45         ResultType const C[] =
46         {
47             1073741824,
48             268435456,
49             16777216,
50             4194304,
51             1638400,
52             802816,
53             451584,
54             278784,
55             184041
56         };
57 
58         ResultType const c2 = 2;
59         ResultType const c4 = 4;
60         ResultType const f = formula::flattening<ResultType>(geometry);
61         ResultType const n = f / (c2 - f);
62         ResultType const ab4 = (get_radius<0>(geometry)
63                                 + get_radius<2>(geometry)) / c4;
64         return geometry::math::pi<ResultType>() * ab4 *
65                  horner_evaluate(n*n, C, C+8) / C[0];
66     }
67 
68 private :
69     //TODO: move the following to a more general space to be used by other
70     //      classes as well
71     /*
72         Evaluate the polynomial in x using Horner's method.
73     */
74     template <typename NT, typename IteratorType>
horner_evaluateboost::geometry::formula_dispatch::quarter_meridian75     static inline NT horner_evaluate(NT x,
76                                      IteratorType begin,
77                                      IteratorType end)
78     {
79         NT result(0);
80         if (begin == end)
81         {
82             return result;
83         }
84         IteratorType it = end;
85         do
86         {
87             result = result * x + *--it;
88         }
89         while (it != begin);
90         return result;
91     }
92 };
93 
94 } // namespace formula_dispatch
95 #endif // DOXYGEN_NO_DISPATCH
96 
97 #ifndef DOXYGEN_NO_DETAIL
98 namespace formula
99 {
100 
101 template <typename ResultType, typename Geometry>
102 ResultType quarter_meridian(Geometry const& geometry)
103 {
104     return formula_dispatch::quarter_meridian<ResultType, Geometry>::apply(geometry);
105 }
106 
107 } // namespace formula
108 #endif // DOXYGEN_NO_DETAIL
109 
110 }} // namespace boost::geometry
111 
112 #endif // BOOST_GEOMETRY_FORMULAS_QUARTER_MERIDIAN_HPP
113