]>
Commit | Line | Data |
---|---|---|
b32b8144 FG |
1 | // Boost.Geometry |
2 | ||
3 | // Copyright (c) 2016 Oracle and/or its affiliates. | |
4 | ||
5 | // Contributed and/or modified by Adam Wulkiewicz, 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_GNOMONIC_SPHEROID_HPP | |
12 | #define BOOST_GEOMETRY_FORMULAS_GNOMONIC_SPHEROID_HPP | |
13 | ||
14 | ||
15 | #include <boost/geometry/core/radius.hpp> | |
16 | ||
17 | #include <boost/geometry/util/condition.hpp> | |
18 | #include <boost/geometry/util/math.hpp> | |
19 | ||
20 | #include <boost/geometry/formulas/andoyer_inverse.hpp> | |
21 | #include <boost/geometry/formulas/flattening.hpp> | |
22 | #include <boost/geometry/formulas/thomas_inverse.hpp> | |
23 | #include <boost/geometry/formulas/vincenty_direct.hpp> | |
24 | #include <boost/geometry/formulas/vincenty_inverse.hpp> | |
25 | ||
26 | ||
27 | namespace boost { namespace geometry { namespace formula | |
28 | { | |
29 | ||
30 | /*! | |
31 | \brief Gnomonic projection on spheroid (ellipsoid of revolution). | |
32 | \author See | |
33 | - Charles F.F Karney, Algorithms for geodesics, 2011 | |
34 | https://arxiv.org/pdf/1109.4448.pdf | |
35 | */ | |
36 | template < | |
37 | typename CT, | |
38 | template <typename, bool, bool, bool, bool ,bool> class Inverse, | |
39 | template <typename, bool, bool, bool, bool> class Direct | |
40 | > | |
41 | class gnomonic_spheroid | |
42 | { | |
43 | typedef Inverse<CT, false, true, true, true, true> inverse_type; | |
44 | typedef typename inverse_type::result_type inverse_result; | |
45 | ||
46 | typedef Direct<CT, false, false, true, true> direct_quantities_type; | |
47 | typedef Direct<CT, true, false, false, false> direct_coordinates_type; | |
48 | typedef typename direct_coordinates_type::result_type direct_result; | |
49 | ||
50 | public: | |
51 | template <typename Spheroid> | |
52 | static inline bool forward(CT const& lon0, CT const& lat0, | |
53 | CT const& lon, CT const& lat, | |
54 | CT & x, CT & y, | |
55 | Spheroid const& spheroid) | |
56 | { | |
57 | inverse_result i_res = inverse_type::apply(lon0, lat0, lon, lat, spheroid); | |
58 | CT const& m = i_res.reduced_length; | |
59 | CT const& M = i_res.geodesic_scale; | |
60 | ||
61 | if (math::smaller_or_equals(M, CT(0))) | |
62 | { | |
63 | return false; | |
64 | } | |
65 | ||
66 | CT rho = m / M; | |
67 | x = sin(i_res.azimuth) * rho; | |
68 | y = cos(i_res.azimuth) * rho; | |
69 | ||
70 | return true; | |
71 | } | |
72 | ||
73 | template <typename Spheroid> | |
74 | static inline bool inverse(CT const& lon0, CT const& lat0, | |
75 | CT const& x, CT const& y, | |
76 | CT & lon, CT & lat, | |
77 | Spheroid const& spheroid) | |
78 | { | |
79 | CT const a = get_radius<0>(spheroid); | |
80 | CT const ds_threshold = a * std::numeric_limits<CT>::epsilon(); // TODO: 0 for non-fundamental type | |
81 | ||
82 | CT const azimuth = atan2(x, y); | |
83 | CT const rho = math::sqrt(math::sqr(x) + math::sqr(y)); // use hypot? | |
84 | CT distance = a * atan(rho / a); | |
85 | ||
86 | bool found = false; | |
87 | for (int i = 0 ; i < 10 ; ++i) | |
88 | { | |
89 | direct_result d_res = direct_quantities_type::apply(lon0, lat0, distance, azimuth, spheroid); | |
90 | CT const& m = d_res.reduced_length; | |
91 | CT const& M = d_res.geodesic_scale; | |
92 | ||
93 | if (math::smaller_or_equals(M, CT(0))) | |
94 | { | |
95 | // found = false; | |
96 | return found; | |
97 | } | |
98 | ||
99 | CT const drho = m / M - rho; // rho = m / M | |
100 | CT const ds = drho * math::sqr(M); // drho/ds = 1/M^2 | |
101 | distance -= ds; | |
102 | ||
103 | // ds_threshold may be 0 | |
104 | if (math::abs(ds) <= ds_threshold) | |
105 | { | |
106 | found = true; | |
107 | break; | |
108 | } | |
109 | } | |
110 | ||
111 | if (found) | |
112 | { | |
113 | direct_result d_res = direct_coordinates_type::apply(lon0, lat0, distance, azimuth, spheroid); | |
114 | lon = d_res.lon2; | |
115 | lat = d_res.lat2; | |
116 | } | |
117 | ||
118 | return found; | |
119 | } | |
120 | }; | |
121 | ||
122 | }}} // namespace boost::geometry::formula | |
123 | ||
124 | ||
125 | #endif // BOOST_GEOMETRY_FORMULAS_GNOMONIC_SPHEROID_HPP |