geometry/formulas/quarter_meridian.hpp

0001 // Boost.Geometry
0002
0003 // Copyright (c) 2018 Oracle and/or its affiliates.
0004
0005 // Contributed and/or modified by Vissarion Fysikopoulos, on behalf of Oracle
0006 // Contributed and/or modified by Adam Wulkiewicz, on behalf of Oracle
0007
0008 // Use, modification and distribution is subject to the Boost Software License,
0009 // Version 1.0. (See accompanying file LICENSE_1_0.txt or copy at
0010 // http://www.boost.org/LICENSE_1_0.txt)
0011
0012 #ifndef BOOST_GEOMETRY_FORMULAS_QUARTER_MERIDIAN_HPP
0013 #define BOOST_GEOMETRY_FORMULAS_QUARTER_MERIDIAN_HPP
0014
0015 #include <boost/geometry/algorithms/not_implemented.hpp>
0016
0017 #include <boost/geometry/core/radius.hpp>
0018 #include <boost/geometry/core/tag.hpp>
0019 #include <boost/geometry/core/tags.hpp>
0020
0021 #include <boost/geometry/formulas/flattening.hpp>
0022
0023 #include <boost/geometry/util/math.hpp>
0024
0025 namespace boost { namespace geometry
0026 {
0027
0028 #ifndef DOXYGEN_NO_DISPATCH
0029 namespace formula_dispatch
0030 {
0031
0032 template <typename ResultType, typename Geometry, typename Tag = typename tag<Geometry>::type>
0033 struct quarter_meridian
0034     : not_implemented<Tag>
0035 {};
0036
0037 template <typename ResultType, typename Geometry>
0038 struct quarter_meridian<ResultType, Geometry, srs_spheroid_tag>
0039 {
0040     //https://en.wikipedia.org/wiki/Meridian_arc#Generalized_series
0041     //http://www.wolframalpha.com/input/?i=(sum(((2*j-3)!!%2F(2*j)!!)%5E2*n%5E(2*j),j,0,8))
0042     static inline ResultType apply(Geometry const& geometry)
0043     {
0044         //order 8 expansion
0045         ResultType const C[] =
0046         {
0047             1073741824,
0048             268435456,
0049             16777216,
0050             4194304,
0051             1638400,
0052             802816,
0053             451584,
0054             278784,
0055             184041
0056         };
0057
0058         ResultType const c2 = 2;
0059         ResultType const c4 = 4;
0060         ResultType const f = formula::flattening<ResultType>(geometry);
0061         ResultType const n = f / (c2 - f);
0062         ResultType const ab4 = (get_radius<0>(geometry)
0063                                 + get_radius<2>(geometry)) / c4;
0064         return geometry::math::pi<ResultType>() * ab4 *
0065                  horner_evaluate(n*n, C, C+8) / C[0];
0066     }
0067
0068 private :
0069     //TODO: move the following to a more general space to be used by other
0070     //      classes as well
0071     /*
0072         Evaluate the polynomial in x using Horner's method.
0073     */
0074     template <typename NT, typename IteratorType>
0075     static inline NT horner_evaluate(NT x,
0076                                      IteratorType begin,
0077                                      IteratorType end)
0078     {
0079         NT result(0);
0080         if (begin == end)
0081         {
0082             return result;
0083         }
0084         IteratorType it = end;
0085         do
0086         {
0087             result = result * x + *--it;
0088         }
0089         while (it != begin);
0090         return result;
0091     }
0092 };
0093
0094 } // namespace formula_dispatch
0095 #endif // DOXYGEN_NO_DISPATCH
0096
0097 #ifndef DOXYGEN_NO_DETAIL
0098 namespace formula
0099 {
0100
0101 template <typename ResultType, typename Geometry>
0102 ResultType quarter_meridian(Geometry const& geometry)
0103 {
0104     return formula_dispatch::quarter_meridian<ResultType, Geometry>::apply(geometry);
0105 }
0106
0107 } // namespace formula
0108 #endif // DOXYGEN_NO_DETAIL
0109
0110 }} // namespace boost::geometry
0111
0112 #endif // BOOST_GEOMETRY_FORMULAS_QUARTER_MERIDIAN_HPP