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