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 { resultboost::geometry::formula::meridian_inverse::result40 result() 41 : distance(0) 42 , meridian(false) 43 {} 44 45 CT distance; 46 bool meridian; 47 }; 48 49 template <typename T> meridian_not_crossing_pole(T lat1,T lat2,CT diff)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 meridian_crossing_pole(CT diff)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> meridian_not_crossing_pole_dist(T lat1,T lat2,Spheroid const & 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> meridian_crossing_pole_dist(T lat1,T lat2,Spheroid const & 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> apply(T lon1,T lat1,T lon2,T lat2,Spheroid const & 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> apply(T lat,Spheroid const & 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 177