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