]>
Commit | Line | Data |
---|---|---|
1 | // Boost.Geometry | |
2 | ||
3 | // Copyright (c) 2017 Oracle and/or its affiliates. | |
4 | ||
5 | // Contributed and/or modified by Vissarion Fysikopoulos, on behalf of Oracle | |
6 | ||
7 | // Use, modification and distribution is subject to the Boost Software License, | |
8 | // Version 1.0. (See accompanying file LICENSE_1_0.txt or copy at | |
9 | // http://www.boost.org/LICENSE_1_0.txt) | |
10 | ||
11 | #ifndef BOOST_GEOMETRY_FORMULAS_ELLIPTIC_ARC_LENGTH_HPP | |
12 | #define BOOST_GEOMETRY_FORMULAS_ELLIPTIC_ARC_LENGTH_HPP | |
13 | ||
14 | #include <boost/math/constants/constants.hpp> | |
15 | ||
16 | #include <boost/geometry/core/radius.hpp> | |
17 | ||
18 | #include <boost/geometry/util/condition.hpp> | |
19 | #include <boost/geometry/util/math.hpp> | |
20 | #include <boost/geometry/util/normalize_spheroidal_coordinates.hpp> | |
21 | ||
22 | #include <boost/geometry/formulas/flattening.hpp> | |
23 | ||
24 | namespace boost { namespace geometry { namespace formula | |
25 | { | |
26 | ||
27 | /*! | |
28 | \brief Compute the arc length of an ellipse. | |
29 | */ | |
30 | ||
31 | template <typename CT, unsigned int Order = 1> | |
32 | class elliptic_arc_length | |
33 | { | |
34 | ||
35 | public : | |
36 | ||
37 | struct result | |
38 | { | |
39 | result() | |
40 | : distance(0) | |
41 | , meridian(false) | |
42 | {} | |
43 | ||
44 | CT distance; | |
45 | bool meridian; | |
46 | }; | |
47 | ||
48 | template <typename T> | |
49 | static bool meridian_not_crossing_pole(T lat1, T lat2, CT diff) | |
50 | { | |
51 | CT half_pi = math::pi<CT>()/CT(2); | |
52 | return math::equals(diff, CT(0)) || | |
53 | (math::equals(lat2, half_pi) && math::equals(lat1, -half_pi)); | |
54 | } | |
55 | ||
56 | static bool meridian_crossing_pole(CT diff) | |
57 | { | |
58 | return math::equals(math::abs(diff), math::pi<CT>()); | |
59 | } | |
60 | ||
61 | ||
62 | template <typename T, typename Spheroid> | |
63 | static CT meridian_not_crossing_pole_dist(T lat1, T lat2, Spheroid const& spheroid) | |
64 | { | |
65 | return math::abs(apply(lat2, spheroid) - apply(lat1, spheroid)); | |
66 | } | |
67 | ||
68 | template <typename T, typename Spheroid> | |
69 | static CT meridian_crossing_pole_dist(T lat1, T lat2, Spheroid const& spheroid) | |
70 | { | |
71 | CT c0 = 0; | |
72 | CT half_pi = math::pi<CT>()/CT(2); | |
73 | CT lat_sign = 1; | |
74 | if (lat1+lat2 < c0) | |
75 | { | |
76 | lat_sign = CT(-1); | |
77 | } | |
78 | return math::abs(lat_sign * CT(2) * apply(half_pi, spheroid) | |
79 | - apply(lat1, spheroid) - apply(lat2, spheroid)); | |
80 | } | |
81 | ||
82 | template <typename T, typename Spheroid> | |
83 | static result apply(T lon1, T lat1, T lon2, T lat2, Spheroid const& spheroid) | |
84 | { | |
85 | result res; | |
86 | ||
87 | CT diff = geometry::math::longitude_distance_signed<geometry::radian>(lon1, lon2); | |
88 | ||
89 | if (lat1 > lat2) | |
90 | { | |
91 | std::swap(lat1, lat2); | |
92 | } | |
93 | ||
94 | if ( meridian_not_crossing_pole(lat1, lat2, diff) ) | |
95 | { | |
96 | res.distance = meridian_not_crossing_pole_dist(lat1, lat2, spheroid); | |
97 | res.meridian = true; | |
98 | } | |
99 | else if ( meridian_crossing_pole(diff) ) | |
100 | { | |
101 | res.distance = meridian_crossing_pole_dist(lat1, lat2, spheroid); | |
102 | res.meridian = true; | |
103 | } | |
104 | return res; | |
105 | } | |
106 | ||
107 | // Distance computation on meridians using series approximations | |
108 | // to elliptic integrals. Formula to compute distance from lattitude 0 to lat | |
109 | // https://en.wikipedia.org/wiki/Meridian_arc | |
110 | // latitudes are assumed to be in radians and in [-pi/2,pi/2] | |
111 | template <typename T, typename Spheroid> | |
112 | static CT apply(T lat, Spheroid const& spheroid) | |
113 | { | |
114 | CT const a = get_radius<0>(spheroid); | |
115 | CT const f = formula::flattening<CT>(spheroid); | |
116 | CT n = f / (CT(2) - f); | |
117 | CT M = a/(1+n); | |
118 | CT C0 = 1; | |
119 | ||
120 | if (Order == 0) | |
121 | { | |
122 | return M * C0 * lat; | |
123 | } | |
124 | ||
125 | CT C2 = -1.5 * n; | |
126 | ||
127 | if (Order == 1) | |
128 | { | |
129 | return M * (C0 * lat + C2 * sin(2*lat)); | |
130 | } | |
131 | ||
132 | CT n2 = n * n; | |
133 | C0 += .25 * n2; | |
134 | CT C4 = 0.9375 * n2; | |
135 | ||
136 | if (Order == 2) | |
137 | { | |
138 | return M * (C0 * lat + C2 * sin(2*lat) + C4 * sin(4*lat)); | |
139 | } | |
140 | ||
141 | CT n3 = n2 * n; | |
142 | C2 += 0.1875 * n3; | |
143 | CT C6 = -0.729166667 * n3; | |
144 | ||
145 | if (Order == 3) | |
146 | { | |
147 | return M * (C0 * lat + C2 * sin(2*lat) + C4 * sin(4*lat) | |
148 | + C6 * sin(6*lat)); | |
149 | } | |
150 | ||
151 | CT n4 = n2 * n2; | |
152 | C4 -= 0.234375 * n4; | |
153 | CT C8 = 0.615234375 * n4; | |
154 | ||
155 | if (Order == 4) | |
156 | { | |
157 | return M * (C0 * lat + C2 * sin(2*lat) + C4 * sin(4*lat) | |
158 | + C6 * sin(6*lat) + C8 * sin(8*lat)); | |
159 | } | |
160 | ||
161 | CT n5 = n4 * n; | |
162 | C6 += 0.227864583 * n5; | |
163 | CT C10 = -0.54140625 * n5; | |
164 | ||
165 | // Order 5 or higher | |
166 | return M * (C0 * lat + C2 * sin(2*lat) + C4 * sin(4*lat) | |
167 | + C6 * sin(6*lat) + C8 * sin(8*lat) + C10 * sin(10*lat)); | |
168 | ||
169 | } | |
170 | ||
171 | // Iterative method to elliptic arc length based on | |
172 | // http://www.codeguru.com/cpp/cpp/algorithms/article.php/c5115/ | |
173 | // Geographic-Distance-and-Azimuth-Calculations.htm | |
174 | // latitudes are assumed to be in radians and in [-pi/2,pi/2] | |
175 | template <typename T1, typename T2, typename Spheroid> | |
176 | CT interative_method(T1 lat1, | |
177 | T2 lat2, | |
178 | Spheroid const& spheroid) | |
179 | { | |
180 | CT result = 0; | |
181 | CT const zero = 0; | |
182 | CT const one = 1; | |
183 | CT const c1 = 2; | |
184 | CT const c2 = 0.5; | |
185 | CT const c3 = 4000; | |
186 | ||
187 | CT const a = get_radius<0>(spheroid); | |
188 | CT const f = formula::flattening<CT>(spheroid); | |
189 | ||
190 | // how many steps to use | |
191 | ||
192 | CT lat1_deg = lat1 * geometry::math::r2d<CT>(); | |
193 | CT lat2_deg = lat2 * geometry::math::r2d<CT>(); | |
194 | ||
195 | int steps = c1 + (c2 + (lat2_deg > lat1_deg) ? CT(lat2_deg - lat1_deg) | |
196 | : CT(lat1_deg - lat2_deg)); | |
197 | steps = (steps > c3) ? c3 : steps; | |
198 | ||
199 | //std::cout << "Steps=" << steps << std::endl; | |
200 | ||
201 | CT snLat1 = sin(lat1); | |
202 | CT snLat2 = sin(lat2); | |
203 | CT twoF = 2 * f - f * f; | |
204 | ||
205 | // limits of integration | |
206 | CT x1 = a * cos(lat1) / | |
207 | sqrt(1 - twoF * snLat1 * snLat1); | |
208 | CT x2 = a * cos(lat2) / | |
209 | sqrt(1 - twoF * snLat2 * snLat2); | |
210 | ||
211 | CT dx = (x2 - x1) / (steps - one); | |
212 | CT x, y1, y2, dy, dydx; | |
213 | CT adx = (dx < zero) ? -dx : dx; // absolute value of dx | |
214 | ||
215 | CT a2 = a * a; | |
216 | CT oneF = 1 - f; | |
217 | ||
218 | // now loop through each step adding up all the little | |
219 | // hypotenuses | |
220 | for (int i = 0; i < (steps - 1); i++){ | |
221 | x = x1 + dx * i; | |
222 | dydx = ((a * oneF * sqrt((one - ((x+dx)*(x+dx))/a2))) - | |
223 | (a * oneF * sqrt((one - (x*x)/a2)))) / dx; | |
224 | result += adx * sqrt(one + dydx*dydx); | |
225 | } | |
226 | ||
227 | return result; | |
228 | } | |
229 | }; | |
230 | ||
231 | }}} // namespace boost::geometry::formula | |
232 | ||
233 | ||
234 | #endif // BOOST_GEOMETRY_FORMULAS_ELLIPTIC_ARC_LENGTH_HPP |