boost/geometry/formulas/karney_direct.hpp

   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.
   9 // Modifications copyright (c) 2018 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>
  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
 222         if (BOOST_GEOMETRY_CONDITION(CalcQuantities))
 223         {
 224             // Evaluate the coefficients for C2.
 225             // Index zero element of coeffs_C2 is unused.
 226             se::coeffs_C2<SeriesOrder, CT> const coeffs_C2(epsilon);
 227
 228             CT const B21 = se::sin_cos_series(sin_sigma1, cos_sigma1, coeffs_C2);
 229             CT const B22 = se::sin_cos_series(sin_sigma2, cos_sigma2, coeffs_C2);
 230
 231             // Find the coefficients for A2 by computing the
 232             // series expansion using Horner scehme.
 233             CT const expansion_A2 = se::evaluate_A2<SeriesOrder>(epsilon);
 234
 235             CT const AB1 = (c1 + expansion_A1) * (B12 - B11);
 236             CT const AB2 = (c1 + expansion_A2) * (B22 - B21);
 237             CT const J12 = (expansion_A1 - expansion_A2) * sigma12 + (AB1 - AB2);
 238
 239             CT const dn1 = math::sqrt(c1 + ep2 * math::sqr(sin_beta1));
 240             CT const dn2 = math::sqrt(c1 + k2 * math::sqr(sin_sigma2));
 241
 242             // Find the reduced length.
 243             result.reduced_length = b * ((dn2 * (cos_sigma1 * sin_sigma2) -
 244                                           dn1 * (sin_sigma1 * cos_sigma2)) -
 245                                           cos_sigma1 * cos_sigma2 * J12);
 246
 247             // Find the geodesic scale.
 248             CT const t = k2 * (sin_sigma2 - sin_sigma1) *
 249                               (sin_sigma2 + sin_sigma1) / (dn1 + dn2);
 250
 251             result.geodesic_scale = cos_sigma12 +
 252                                     (t * sin_sigma2 - cos_sigma2 * J12) *
 253                                     sin_sigma1 / dn1;
 254         }
 255
 256         return result;
 257     }
 258 };
 259
 260 }}} // namespace boost::geometry::formula
 261
 262
 263 #endif // BOOST_GEOMETRY_FORMULAS_KARNEY_DIRECT_HPP