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

   1 // Boost.Geometry
   2
   3 // Copyright (c) 2017-2018 Oracle and/or its affiliates.
   4
   5 // Contributed and/or modified by Vissarion Fysikopoulos, 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_MERIDIAN_INVERSE_HPP
  12 #define BOOST_GEOMETRY_FORMULAS_MERIDIAN_INVERSE_HPP
  13
  14 #include <boost/math/constants/constants.hpp>
  15
  16 #include <boost/geometry/core/radius.hpp>
  17
  18 #include <boost/geometry/util/condition.hpp>
  19 #include <boost/geometry/util/math.hpp>
  20 #include <boost/geometry/util/normalize_spheroidal_coordinates.hpp>
  21
  22 #include <boost/geometry/formulas/flattening.hpp>
  23 #include <boost/geometry/formulas/meridian_segment.hpp>
  24
  25 namespace boost { namespace geometry { namespace formula
  26 {
  27
  28 /*!
  29 \brief Compute the arc length of an ellipse.
  30 */
  31
  32 template <typename CT, unsigned int Order = 1>
  33 class meridian_inverse
  34 {
  35
  36 public :
  37
  38     struct result
  39     {
  40         result()
  41             : distance(0)
  42             , meridian(false)
  43         {}
  44
  45         CT distance;
  46         bool meridian;
  47     };
  48
  49     template <typename T>
  50     static bool meridian_not_crossing_pole(T lat1, T lat2, CT diff)
  51     {
  52         CT half_pi = math::pi<CT>()/CT(2);
  53         return math::equals(diff, CT(0)) ||
  54                     (math::equals(lat2, half_pi) && math::equals(lat1, -half_pi));
  55     }
  56
  57     static bool meridian_crossing_pole(CT diff)
  58     {
  59         return math::equals(math::abs(diff), math::pi<CT>());
  60     }
  61
  62
  63     template <typename T, typename Spheroid>
  64     static CT meridian_not_crossing_pole_dist(T lat1, T lat2, Spheroid const& spheroid)
  65     {
  66         return math::abs(apply(lat2, spheroid) - apply(lat1, spheroid));
  67     }
  68
  69     template <typename T, typename Spheroid>
  70     static CT meridian_crossing_pole_dist(T lat1, T lat2, Spheroid const& spheroid)
  71     {
  72         CT c0 = 0;
  73         CT half_pi = math::pi<CT>()/CT(2);
  74         CT lat_sign = 1;
  75         if (lat1+lat2 < c0)
  76         {
  77             lat_sign = CT(-1);
  78         }
  79         return math::abs(lat_sign * CT(2) * apply(half_pi, spheroid)
  80                          - apply(lat1, spheroid) - apply(lat2, spheroid));
  81     }
  82
  83     template <typename T, typename Spheroid>
  84     static result apply(T lon1, T lat1, T lon2, T lat2, Spheroid const& spheroid)
  85     {
  86         result res;
  87
  88         CT diff = geometry::math::longitude_distance_signed<geometry::radian>(lon1, lon2);
  89
  90         if (lat1 > lat2)
  91         {
  92             std::swap(lat1, lat2);
  93         }
  94
  95         if ( meridian_not_crossing_pole(lat1, lat2, diff) )
  96         {
  97             res.distance = meridian_not_crossing_pole_dist(lat1, lat2, spheroid);
  98             res.meridian = true;
  99         }
 100         else if ( meridian_crossing_pole(diff) )
 101         {
 102             res.distance = meridian_crossing_pole_dist(lat1, lat2, spheroid);
 103             res.meridian = true;
 104         }
 105         return res;
 106     }
 107
 108     // Distance computation on meridians using series approximations
 109     // to elliptic integrals. Formula to compute distance from lattitude 0 to lat
 110     // https://en.wikipedia.org/wiki/Meridian_arc
 111     // latitudes are assumed to be in radians and in [-pi/2,pi/2]
 112     template <typename T, typename Spheroid>
 113     static CT apply(T lat, Spheroid const& spheroid)
 114     {
 115         CT const a = get_radius<0>(spheroid);
 116         CT const f = formula::flattening<CT>(spheroid);
 117         CT n = f / (CT(2) - f);
 118         CT M = a/(1+n);
 119         CT C0 = 1;
 120
 121         if (Order == 0)
 122         {
 123            return M * C0 * lat;
 124         }
 125
 126         CT C2 = -1.5 * n;
 127
 128         if (Order == 1)
 129         {
 130             return M * (C0 * lat + C2 * sin(2*lat));
 131         }
 132
 133         CT n2 = n * n;
 134         C0 += .25 * n2;
 135         CT C4 = 0.9375 * n2;
 136
 137         if (Order == 2)
 138         {
 139             return M * (C0 * lat + C2 * sin(2*lat) + C4 * sin(4*lat));
 140         }
 141
 142         CT n3 = n2 * n;
 143         C2 += 0.1875 * n3;
 144         CT C6 = -0.729166667 * n3;
 145
 146         if (Order == 3)
 147         {
 148             return M * (C0 * lat + C2 * sin(2*lat) + C4 * sin(4*lat)
 149                       + C6 * sin(6*lat));
 150         }
 151
 152         CT n4 = n2 * n2;
 153         C4 -= 0.234375 * n4;
 154         CT C8 = 0.615234375 * n4;
 155
 156         if (Order == 4)
 157         {
 158             return M * (C0 * lat + C2 * sin(2*lat) + C4 * sin(4*lat)
 159                       + C6 * sin(6*lat) + C8 * sin(8*lat));
 160         }
 161
 162         CT n5 = n4 * n;
 163         C6 += 0.227864583 * n5;
 164         CT C10 = -0.54140625 * n5;
 165
 166         // Order 5 or higher
 167         return M * (C0 * lat + C2 * sin(2*lat) + C4 * sin(4*lat)
 168                   + C6 * sin(6*lat) + C8 * sin(8*lat) + C10 * sin(10*lat));
 169
 170     }
 171 };
 172
 173 }}} // namespace boost::geometry::formula
 174
 175
 176 #endif // BOOST_GEOMETRY_FORMULAS_MERIDIAN_INVERSE_HPP