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>
forward(CT const & lon0,CT const & lat0,CT const & lon,CT const & lat,CT & x,CT & y,Spheroid const & 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>
inverse(CT const & lon0,CT const & lat0,CT const & x,CT const & y,CT & lon,CT & lat,Spheroid const & 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
126