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