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