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