vertex_latitude.hpp 4.4 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147
  1. // Boost.Geometry
  2. // Copyright (c) 2016-2020 Oracle and/or its affiliates.
  3. // Contributed and/or modified by Vissarion Fysikopoulos, on behalf of Oracle
  4. // Contributed and/or modified by Adam Wulkiewicz, on behalf of Oracle
  5. // Use, modification and distribution is subject to the Boost Software License,
  6. // Version 1.0. (See accompanying file LICENSE_1_0.txt or copy at
  7. // http://www.boost.org/LICENSE_1_0.txt)
  8. #ifndef BOOST_GEOMETRY_FORMULAS_MAXIMUM_LATITUDE_HPP
  9. #define BOOST_GEOMETRY_FORMULAS_MAXIMUM_LATITUDE_HPP
  10. #include <boost/geometry/core/static_assert.hpp>
  11. #include <boost/geometry/formulas/flattening.hpp>
  12. #include <boost/geometry/formulas/spherical.hpp>
  13. namespace boost { namespace geometry { namespace formula
  14. {
  15. /*!
  16. \brief Algorithm to compute the vertex latitude of a geodesic segment. Vertex is
  17. a point on the geodesic that maximizes (or minimizes) the latitude.
  18. \author See
  19. [Wood96] Wood - Vertex Latitudes on Ellipsoid Geodesics, SIAM Rev., 38(4),
  20. 637–644, 1996
  21. */
  22. template <typename CT>
  23. class vertex_latitude_on_sphere
  24. {
  25. public:
  26. template<typename T1, typename T2>
  27. static inline CT apply(T1 const& lat1,
  28. T2 const& alp1)
  29. {
  30. return std::acos( math::abs(cos(lat1) * sin(alp1)) );
  31. }
  32. };
  33. template <typename CT>
  34. class vertex_latitude_on_spheroid
  35. {
  36. public:
  37. /*
  38. * formula based on paper
  39. * [Wood96] Wood - Vertex Latitudes on Ellipsoid Geodesics, SIAM Rev., 38(4),
  40. * 637–644, 1996
  41. template <typename T1, typename T2, typename Spheroid>
  42. static inline CT apply(T1 const& lat1,
  43. T2 const& alp1,
  44. Spheroid const& spheroid)
  45. {
  46. CT const f = formula::flattening<CT>(spheroid);
  47. CT const e2 = f * (CT(2) - f);
  48. CT const sin_alp1 = sin(alp1);
  49. CT const sin2_lat1 = math::sqr(sin(lat1));
  50. CT const cos2_lat1 = CT(1) - sin2_lat1;
  51. CT const e2_sin2 = CT(1) - e2 * sin2_lat1;
  52. CT const cos2_sin2 = cos2_lat1 * math::sqr(sin_alp1);
  53. CT const vertex_lat = std::asin( math::sqrt((e2_sin2 - cos2_sin2)
  54. / (e2_sin2 - e2 * cos2_sin2)));
  55. return vertex_lat;
  56. }
  57. */
  58. // simpler formula based on Clairaut relation for spheroids
  59. template <typename T1, typename T2, typename Spheroid>
  60. static inline CT apply(T1 const& lat1,
  61. T2 const& alp1,
  62. Spheroid const& spheroid)
  63. {
  64. CT const f = formula::flattening<CT>(spheroid);
  65. CT const one_minus_f = (CT(1) - f);
  66. //get the reduced latitude
  67. CT const bet1 = atan( one_minus_f * tan(lat1) );
  68. //apply Clairaut relation
  69. CT const betv = vertex_latitude_on_sphere<CT>::apply(bet1, alp1);
  70. //return the spheroid latitude
  71. return atan( tan(betv) / one_minus_f );
  72. }
  73. /*
  74. template <typename T>
  75. inline static void sign_adjustment(CT lat1, CT lat2, CT vertex_lat, T& vrt_result)
  76. {
  77. // signbit returns a non-zero value (true) if the sign is negative;
  78. // and zero (false) otherwise.
  79. bool sign = std::signbit(std::abs(lat1) > std::abs(lat2) ? lat1 : lat2);
  80. vrt_result.north = sign ? std::max(lat1, lat2) : vertex_lat;
  81. vrt_result.south = sign ? vertex_lat * CT(-1) : std::min(lat1, lat2);
  82. }
  83. template <typename T>
  84. inline static bool vertex_on_segment(CT alp1, CT alp2, CT lat1, CT lat2, T& vrt_result)
  85. {
  86. CT const half_pi = math::pi<CT>() / CT(2);
  87. // if the segment does not contain the vertex of the geodesic
  88. // then return the endpoint of max (min) latitude
  89. if ((alp1 < half_pi && alp2 < half_pi)
  90. || (alp1 > half_pi && alp2 > half_pi))
  91. {
  92. vrt_result.north = std::max(lat1, lat2);
  93. vrt_result.south = std::min(lat1, lat2);
  94. return false;
  95. }
  96. return true;
  97. }
  98. */
  99. };
  100. template <typename CT, typename CS_Tag>
  101. struct vertex_latitude
  102. {
  103. BOOST_GEOMETRY_STATIC_ASSERT_FALSE(
  104. "Not implemented for this coordinate system.",
  105. CT, CS_Tag);
  106. };
  107. template <typename CT>
  108. struct vertex_latitude<CT, spherical_equatorial_tag>
  109. : vertex_latitude_on_sphere<CT>
  110. {};
  111. template <typename CT>
  112. struct vertex_latitude<CT, geographic_tag>
  113. : vertex_latitude_on_spheroid<CT>
  114. {};
  115. }}} // namespace boost::geometry::formula
  116. #endif // BOOST_GEOMETRY_FORMULAS_MAXIMUM_LATITUDE_HPP