quarter_meridian.hpp 2.9 KB

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