geometry/formulas/gnomonic_spheroid.hpp

0001 // Boost.Geometry
0002
0003 // Copyright (c) 2016 Oracle and/or its affiliates.
0004
0005 // Contributed and/or modified by Adam Wulkiewicz, on behalf of Oracle
0006
0007 // Use, modification and distribution is subject to the Boost Software License,
0008 // Version 1.0. (See accompanying file LICENSE_1_0.txt or copy at
0009 // http://www.boost.org/LICENSE_1_0.txt)
0010
0011 #ifndef BOOST_GEOMETRY_FORMULAS_GNOMONIC_SPHEROID_HPP
0012 #define BOOST_GEOMETRY_FORMULAS_GNOMONIC_SPHEROID_HPP
0013
0014
0015 #include <boost/geometry/core/radius.hpp>
0016
0017 #include <boost/geometry/util/condition.hpp>
0018 #include <boost/geometry/util/math.hpp>
0019
0020 #include <boost/geometry/formulas/andoyer_inverse.hpp>
0021 #include <boost/geometry/formulas/flattening.hpp>
0022 #include <boost/geometry/formulas/thomas_inverse.hpp>
0023 #include <boost/geometry/formulas/vincenty_direct.hpp>
0024 #include <boost/geometry/formulas/vincenty_inverse.hpp>
0025
0026
0027 namespace boost { namespace geometry { namespace formula
0028 {
0029
0030 /*!
0031 \brief Gnomonic projection on spheroid (ellipsoid of revolution).
0032 \author See
0033 - Charles F.F Karney, Algorithms for geodesics, 2011
0034 https://arxiv.org/pdf/1109.4448.pdf
0035 */
0036 template <
0037     typename CT,
0038     template <typename, bool, bool, bool, bool ,bool> class Inverse,
0039     template <typename, bool, bool, bool, bool> class Direct
0040 >
0041 class gnomonic_spheroid
0042 {
0043     typedef Inverse<CT, false, true, true, true, true> inverse_type;
0044     typedef typename inverse_type::result_type inverse_result;
0045
0046     typedef Direct<CT, false, false, true, true> direct_quantities_type;
0047     typedef Direct<CT, true, false, false, false> direct_coordinates_type;
0048     typedef typename direct_coordinates_type::result_type direct_result;
0049
0050 public:
0051     template <typename Spheroid>
0052     static inline bool forward(CT const& lon0, CT const& lat0,
0053                                CT const& lon, CT const& lat,
0054                                CT & x, CT & y,
0055                                Spheroid const& spheroid)
0056     {
0057         inverse_result i_res = inverse_type::apply(lon0, lat0, lon, lat, spheroid);
0058         CT const& m = i_res.reduced_length;
0059         CT const& M = i_res.geodesic_scale;
0060
0061         if (math::smaller_or_equals(M, CT(0)))
0062         {
0063             return false;
0064         }
0065
0066         CT rho = m / M;
0067         x = sin(i_res.azimuth) * rho;
0068         y = cos(i_res.azimuth) * rho;
0069
0070         return true;
0071     }
0072
0073     template <typename Spheroid>
0074     static inline bool inverse(CT const& lon0, CT const& lat0,
0075                                CT const& x, CT const& y,
0076                                CT & lon, CT & lat,
0077                                Spheroid const& spheroid)
0078     {
0079         CT const a = get_radius<0>(spheroid);
0080         CT const ds_threshold = a * std::numeric_limits<CT>::epsilon(); // TODO: 0 for non-fundamental type
0081
0082         CT const azimuth = atan2(x, y);
0083         CT const rho = math::sqrt(math::sqr(x) + math::sqr(y)); // use hypot?
0084         CT distance = a * atan(rho / a);
0085
0086         bool found = false;
0087         for (int i = 0 ; i < 10 ; ++i)
0088         {
0089             direct_result d_res = direct_quantities_type::apply(lon0, lat0, distance, azimuth, spheroid);
0090             CT const& m = d_res.reduced_length;
0091             CT const& M = d_res.geodesic_scale;
0092
0093             if (math::smaller_or_equals(M, CT(0)))
0094             {
0095                 // found = false;
0096                 return found;
0097             }
0098
0099             CT const drho = m / M - rho; // rho = m / M
0100             CT const ds = drho * math::sqr(M); // drho/ds = 1/M^2
0101             distance -= ds;
0102
0103             // ds_threshold may be 0
0104             if (math::abs(ds) <= ds_threshold)
0105             {
0106                 found = true;
0107                 break;
0108             }
0109         }
0110
0111         if (found)
0112         {
0113             direct_result d_res = direct_coordinates_type::apply(lon0, lat0, distance, azimuth, spheroid);
0114             lon = d_res.lon2;
0115             lat = d_res.lat2;
0116         }
0117
0118         return found;
0119     }
0120 };
0121
0122 }}} // namespace boost::geometry::formula
0123
0124
0125 #endif // BOOST_GEOMETRY_FORMULAS_GNOMONIC_SPHEROID_HPP