ceph/src/boost/boost/geometry/formulas/gnomonic_spheroid.hpp

   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