1 // Boost.Geometry 2 3 // Copyright (c) 2018 Adeel Ahmad, Islamabad, Pakistan. 4 5 // Contributed and/or modified by Adeel Ahmad, 6 // as part of Google Summer of Code 2018 program. 7 8 // This file was modified by Oracle on 2018-2020. 9 // Modifications copyright (c) 2018-2020 Oracle and/or its affiliates. 10 // Contributed and/or modified by Adam Wulkiewicz, on behalf of Oracle 11 12 // Use, modification and distribution is subject to the Boost Software License, 13 // Version 1.0. (See accompanying file LICENSE_1_0.txt or copy at 14 // http://www.boost.org/LICENSE_1_0.txt) 15 16 // This file is converted from GeographicLib, https://geographiclib.sourceforge.io 17 // GeographicLib is originally written by Charles Karney. 18 19 // Author: Charles Karney (2008-2017) 20 21 // Last updated version of GeographicLib: 1.49 22 23 // Original copyright notice: 24 25 // Copyright (c) Charles Karney (2008-2017) <charles@karney.com> and licensed 26 // under the MIT/X11 License. For more information, see 27 // https://geographiclib.sourceforge.io 28 29 #ifndef BOOST_GEOMETRY_FORMULAS_KARNEY_DIRECT_HPP 30 #define BOOST_GEOMETRY_FORMULAS_KARNEY_DIRECT_HPP 31 32 33 #include <boost/array.hpp> 34 35 #include <boost/math/constants/constants.hpp> 36 #include <boost/math/special_functions/hypot.hpp> 37 38 #include <boost/geometry/formulas/flattening.hpp> 39 #include <boost/geometry/formulas/result_direct.hpp> 40 41 #include <boost/geometry/util/condition.hpp> 42 #include <boost/geometry/util/math.hpp> 43 #include <boost/geometry/util/normalize_spheroidal_coordinates.hpp> 44 #include <boost/geometry/util/series_expansion.hpp> 45 46 47 namespace boost { namespace geometry { namespace formula 48 { 49 50 namespace se = series_expansion; 51 52 /*! 53 \brief The solution of the direct problem of geodesics on latlong coordinates, 54 after Karney (2011). 55 \author See 56 - Charles F.F Karney, Algorithms for geodesics, 2011 57 https://arxiv.org/pdf/1109.4448.pdf 58 */ 59 template < 60 typename CT, 61 bool EnableCoordinates = true, 62 bool EnableReverseAzimuth = false, 63 bool EnableReducedLength = false, 64 bool EnableGeodesicScale = false, 65 size_t SeriesOrder = 8 66 > 67 class karney_direct 68 { 69 static const bool CalcQuantities = EnableReducedLength || EnableGeodesicScale; 70 static const bool CalcCoordinates = EnableCoordinates || CalcQuantities; 71 static const bool CalcRevAzimuth = EnableReverseAzimuth || CalcCoordinates || CalcQuantities; 72 73 public: 74 typedef result_direct<CT> result_type; 75 76 template <typename T, typename Dist, typename Azi, typename Spheroid> apply(T const & lo1,T const & la1,Dist const & distance,Azi const & azimuth12,Spheroid const & spheroid)77 static inline result_type apply(T const& lo1, 78 T const& la1, 79 Dist const& distance, 80 Azi const& azimuth12, 81 Spheroid const& spheroid) 82 { 83 result_type result; 84 85 CT lon1 = lo1; 86 CT const lat1 = la1; 87 88 Azi azi12 = azimuth12; 89 math::normalize_azimuth<degree, Azi>(azi12); 90 91 CT const c0 = 0; 92 CT const c1 = 1; 93 CT const c2 = 2; 94 95 CT const b = CT(get_radius<2>(spheroid)); 96 CT const f = formula::flattening<CT>(spheroid); 97 CT const one_minus_f = c1 - f; 98 CT const two_minus_f = c2 - f; 99 100 CT const n = f / two_minus_f; 101 CT const e2 = f * two_minus_f; 102 CT const ep2 = e2 / math::sqr(one_minus_f); 103 104 CT sin_alpha1, cos_alpha1; 105 math::sin_cos_degrees<CT>(azi12, sin_alpha1, cos_alpha1); 106 107 // Find the reduced latitude. 108 CT sin_beta1, cos_beta1; 109 math::sin_cos_degrees<CT>(lat1, sin_beta1, cos_beta1); 110 sin_beta1 *= one_minus_f; 111 112 math::normalize_unit_vector<CT>(sin_beta1, cos_beta1); 113 114 cos_beta1 = (std::max)(c0, cos_beta1); 115 116 // Obtain alpha 0 by solving the spherical triangle. 117 CT const sin_alpha0 = sin_alpha1 * cos_beta1; 118 CT const cos_alpha0 = boost::math::hypot(cos_alpha1, sin_alpha1 * sin_beta1); 119 120 CT const k2 = math::sqr(cos_alpha0) * ep2; 121 122 CT const epsilon = k2 / (c2 * (c1 + math::sqrt(c1 + k2)) + k2); 123 124 // Find the coefficients for A1 by computing the 125 // series expansion using Horner scehme. 126 CT const expansion_A1 = se::evaluate_A1<SeriesOrder>(epsilon); 127 128 // Index zero element of coeffs_C1 is unused. 129 se::coeffs_C1<SeriesOrder, CT> const coeffs_C1(epsilon); 130 131 // Tau is an integration variable. 132 CT const tau12 = distance / (b * (c1 + expansion_A1)); 133 134 CT const sin_tau12 = sin(tau12); 135 CT const cos_tau12 = cos(tau12); 136 137 CT sin_sigma1 = sin_beta1; 138 CT sin_omega1 = sin_alpha0 * sin_beta1; 139 140 CT cos_sigma1, cos_omega1; 141 cos_sigma1 = cos_omega1 = sin_beta1 != c0 || cos_alpha1 != c0 ? cos_beta1 * cos_alpha1 : c1; 142 math::normalize_unit_vector<CT>(sin_sigma1, cos_sigma1); 143 144 CT const B11 = se::sin_cos_series(sin_sigma1, cos_sigma1, coeffs_C1); 145 CT const sin_B11 = sin(B11); 146 CT const cos_B11 = cos(B11); 147 148 CT const sin_tau1 = sin_sigma1 * cos_B11 + cos_sigma1 * sin_B11; 149 CT const cos_tau1 = cos_sigma1 * cos_B11 - sin_sigma1 * sin_B11; 150 151 // Index zero element of coeffs_C1p is unused. 152 se::coeffs_C1p<SeriesOrder, CT> const coeffs_C1p(epsilon); 153 154 CT const B12 = - se::sin_cos_series 155 (sin_tau1 * cos_tau12 + cos_tau1 * sin_tau12, 156 cos_tau1 * cos_tau12 - sin_tau1 * sin_tau12, 157 coeffs_C1p); 158 159 CT const sigma12 = tau12 - (B12 - B11); 160 CT const sin_sigma12 = sin(sigma12); 161 CT const cos_sigma12 = cos(sigma12); 162 163 CT const sin_sigma2 = sin_sigma1 * cos_sigma12 + cos_sigma1 * sin_sigma12; 164 CT const cos_sigma2 = cos_sigma1 * cos_sigma12 - sin_sigma1 * sin_sigma12; 165 166 if (BOOST_GEOMETRY_CONDITION(CalcRevAzimuth)) 167 { 168 CT const sin_alpha2 = sin_alpha0; 169 CT const cos_alpha2 = cos_alpha0 * cos_sigma2; 170 171 result.reverse_azimuth = atan2(sin_alpha2, cos_alpha2); 172 173 // Convert the angle to radians. 174 result.reverse_azimuth /= math::d2r<CT>(); 175 } 176 177 if (BOOST_GEOMETRY_CONDITION(CalcCoordinates)) 178 { 179 // Find the latitude at the second point. 180 CT const sin_beta2 = cos_alpha0 * sin_sigma2; 181 CT const cos_beta2 = boost::math::hypot(sin_alpha0, cos_alpha0 * cos_sigma2); 182 183 result.lat2 = atan2(sin_beta2, one_minus_f * cos_beta2); 184 185 // Convert the coordinate to radians. 186 result.lat2 /= math::d2r<CT>(); 187 188 // Find the longitude at the second point. 189 CT const sin_omega2 = sin_alpha0 * sin_sigma2; 190 CT const cos_omega2 = cos_sigma2; 191 192 CT const omega12 = atan2(sin_omega2 * cos_omega1 - cos_omega2 * sin_omega1, 193 cos_omega2 * cos_omega1 + sin_omega2 * sin_omega1); 194 195 se::coeffs_A3<SeriesOrder, CT> const coeffs_A3(n); 196 197 CT const A3 = math::horner_evaluate(epsilon, coeffs_A3.begin(), coeffs_A3.end()); 198 CT const A3c = -f * sin_alpha0 * A3; 199 200 se::coeffs_C3<SeriesOrder, CT> const coeffs_C3(n, epsilon); 201 202 CT const B31 = se::sin_cos_series(sin_sigma1, cos_sigma1, coeffs_C3); 203 204 CT const lam12 = omega12 + A3c * 205 (sigma12 + (se::sin_cos_series 206 (sin_sigma2, 207 cos_sigma2, 208 coeffs_C3) - B31)); 209 210 // Convert to radians to get the 211 // longitudinal difference. 212 CT lon12 = lam12 / math::d2r<CT>(); 213 214 // Add the longitude at first point to the longitudinal 215 // difference and normalize the result. 216 math::normalize_longitude<degree, CT>(lon1); 217 math::normalize_longitude<degree, CT>(lon12); 218 219 result.lon2 = lon1 + lon12; 220 221 // For longitudes close to the antimeridian the result can be out 222 // of range. Therefore normalize. 223 // In other formulas this has to be done at the end because 224 // otherwise differential quantities are calculated incorrectly. 225 // But here it's ok since result.lon2 is not used after this point. 226 math::normalize_longitude<degree, CT>(result.lon2); 227 } 228 229 if (BOOST_GEOMETRY_CONDITION(CalcQuantities)) 230 { 231 // Evaluate the coefficients for C2. 232 // Index zero element of coeffs_C2 is unused. 233 se::coeffs_C2<SeriesOrder, CT> const coeffs_C2(epsilon); 234 235 CT const B21 = se::sin_cos_series(sin_sigma1, cos_sigma1, coeffs_C2); 236 CT const B22 = se::sin_cos_series(sin_sigma2, cos_sigma2, coeffs_C2); 237 238 // Find the coefficients for A2 by computing the 239 // series expansion using Horner scehme. 240 CT const expansion_A2 = se::evaluate_A2<SeriesOrder>(epsilon); 241 242 CT const AB1 = (c1 + expansion_A1) * (B12 - B11); 243 CT const AB2 = (c1 + expansion_A2) * (B22 - B21); 244 CT const J12 = (expansion_A1 - expansion_A2) * sigma12 + (AB1 - AB2); 245 246 CT const dn1 = math::sqrt(c1 + ep2 * math::sqr(sin_beta1)); 247 CT const dn2 = math::sqrt(c1 + k2 * math::sqr(sin_sigma2)); 248 249 // Find the reduced length. 250 result.reduced_length = b * ((dn2 * (cos_sigma1 * sin_sigma2) - 251 dn1 * (sin_sigma1 * cos_sigma2)) - 252 cos_sigma1 * cos_sigma2 * J12); 253 254 // Find the geodesic scale. 255 CT const t = k2 * (sin_sigma2 - sin_sigma1) * 256 (sin_sigma2 + sin_sigma1) / (dn1 + dn2); 257 258 result.geodesic_scale = cos_sigma12 + 259 (t * sin_sigma2 - cos_sigma2 * J12) * 260 sin_sigma1 / dn1; 261 } 262 263 return result; 264 } 265 }; 266 267 }}} // namespace boost::geometry::formula 268 269 270 #endif // BOOST_GEOMETRY_FORMULAS_KARNEY_DIRECT_HPP 271