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>
apply(T1 const & lat1,T2 const & alp1)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>
apply(T1 const & lat1,T2 const & alp1,Spheroid const & 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
149