ceph/src/boost/boost/geometry/formulas/vertex_latitude.hpp

   1 // Boost.Geometry
   2
   3 // Copyright (c) 2016-2017 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_MAXIMUM_LATITUDE_HPP
  13 #define BOOST_GEOMETRY_FORMULAS_MAXIMUM_LATITUDE_HPP
  14
  15 #include <boost/geometry/core/srs.hpp>
  16 #include <boost/geometry/formulas/flattening.hpp>
  17 #include <boost/geometry/formulas/spherical.hpp>
  18 #include <boost/mpl/assert.hpp>
  19
  20 namespace boost { namespace geometry { namespace formula
  21 {
  22
  23 /*!
  24 \brief Algorithm to compute the vertex latitude of a geodesic segment. Vertex is
  25 a point on the geodesic that maximizes (or minimizes) the latitude.
  26 \author See
  27     [Wood96] Wood - Vertex Latitudes on Ellipsoid Geodesics, SIAM Rev., 38(4),
  28              637–644, 1996
  29 */
  30
  31 template <typename CT>
  32 class vertex_latitude_on_sphere
  33 {
  34
  35 public:
  36     template<typename T1, typename T2>
  37     static inline CT apply(T1 const& lat1,
  38                            T2 const& alp1)
  39     {
  40         return std::acos( math::abs(cos(lat1) * sin(alp1)) );
  41     }
  42 };
  43
  44 template <typename CT>
  45 class vertex_latitude_on_spheroid
  46 {
  47
  48 public:
  49 /*
  50  * formula based on paper
  51  *   [Wood96] Wood - Vertex Latitudes on Ellipsoid Geodesics, SIAM Rev., 38(4),
  52  *            637–644, 1996
  53     template <typename T1, typename T2, typename Spheroid>
  54     static inline CT apply(T1 const& lat1,
  55                            T2 const& alp1,
  56                            Spheroid const& spheroid)
  57     {
  58         CT const f = formula::flattening<CT>(spheroid);
  59
  60         CT const e2 = f * (CT(2) - f);
  61         CT const sin_alp1 = sin(alp1);
  62         CT const sin2_lat1 = math::sqr(sin(lat1));
  63         CT const cos2_lat1 = CT(1) - sin2_lat1;
  64
  65         CT const e2_sin2 = CT(1) - e2 * sin2_lat1;
  66         CT const cos2_sin2 = cos2_lat1 * math::sqr(sin_alp1);
  67         CT const vertex_lat = std::asin( math::sqrt((e2_sin2 - cos2_sin2)
  68                                                     / (e2_sin2 - e2 * cos2_sin2)));
  69         return vertex_lat;
  70     }
  71 */
  72
  73     // simpler formula based on Clairaut relation for spheroids
  74     template <typename T1, typename T2, typename Spheroid>
  75     static inline CT apply(T1 const& lat1,
  76                            T2 const& alp1,
  77                            Spheroid const& spheroid)
  78     {
  79         CT const f = formula::flattening<CT>(spheroid);
  80
  81         CT const one_minus_f = (CT(1) - f);
  82
  83         //get the reduced latitude
  84         CT const bet1 = atan( one_minus_f * tan(lat1) );
  85
  86         //apply Clairaut relation
  87         CT const betv =  vertex_latitude_on_sphere<CT>::apply(bet1, alp1);
  88
  89         //return the spheroid latitude
  90         return atan( tan(betv) / one_minus_f );
  91     }
  92
  93     /*
  94     template <typename T>
  95     inline static void sign_adjustment(CT lat1, CT lat2, CT vertex_lat, T& vrt_result)
  96     {
  97         // signbit returns a non-zero value (true) if the sign is negative;
  98         // and zero (false) otherwise.
  99         bool sign = std::signbit(std::abs(lat1) > std::abs(lat2) ? lat1 : lat2);
 100
 101         vrt_result.north = sign ? std::max(lat1, lat2) : vertex_lat;
 102         vrt_result.south = sign ? vertex_lat * CT(-1) : std::min(lat1, lat2);
 103     }
 104
 105     template <typename T>
 106     inline static bool vertex_on_segment(CT alp1, CT alp2, CT lat1, CT lat2, T& vrt_result)
 107     {
 108         CT const half_pi = math::pi<CT>() / CT(2);
 109
 110         // if the segment does not contain the vertex of the geodesic
 111         // then return the endpoint of max (min) latitude
 112         if ((alp1 < half_pi && alp2 < half_pi)
 113                 || (alp1 > half_pi && alp2 > half_pi))
 114         {
 115             vrt_result.north = std::max(lat1, lat2);
 116             vrt_result.south = std::min(lat1, lat2);
 117             return false;
 118         }
 119         return true;
 120     }
 121     */
 122 };
 123
 124
 125 template <typename CT, typename CS_Tag>
 126 struct vertex_latitude
 127 {
 128     BOOST_MPL_ASSERT_MSG
 129          (
 130              false, NOT_IMPLEMENTED_FOR_THIS_COORDINATE_SYSTEM, (types<CS_Tag>)
 131          );
 132
 133 };
 134
 135 template <typename CT>
 136 struct vertex_latitude<CT, spherical_equatorial_tag>
 137         : vertex_latitude_on_sphere<CT>
 138 {};
 139
 140 template <typename CT>
 141 struct vertex_latitude<CT, geographic_tag>
 142         : vertex_latitude_on_spheroid<CT>
 143 {};
 144
 145
 146 }}} // namespace boost::geometry::formula
 147
 148 #endif // BOOST_GEOMETRY_FORMULAS_MAXIMUM_LATITUDE_HPP