geometry/formulas/karney_direct.hpp

0001 // Boost.Geometry
0002
0003 // Copyright (c) 2018 Adeel Ahmad, Islamabad, Pakistan.
0004 // Copyright (c) 2023 Adam Wulkiewicz, Lodz, Poland.
0005
0006 // Contributed and/or modified by Adeel Ahmad,
0007 //   as part of Google Summer of Code 2018 program.
0008
0009 // This file was modified by Oracle on 2018-2022.
0010 // Modifications copyright (c) 2018-2022 Oracle and/or its affiliates.
0011 // Contributed and/or modified by Vissarion Fysikopoulos, on behalf of Oracle
0012 // Contributed and/or modified by Adam Wulkiewicz, on behalf of Oracle
0013
0014 // Use, modification and distribution is subject to the Boost Software License,
0015 // Version 1.0. (See accompanying file LICENSE_1_0.txt or copy at
0016 // http://www.boost.org/LICENSE_1_0.txt)
0017
0018 // This file is converted from GeographicLib, https://geographiclib.sourceforge.io
0019 // GeographicLib is originally written by Charles Karney.
0020
0021 // Author: Charles Karney (2008-2017)
0022
0023 // Last updated version of GeographicLib: 1.49
0024
0025 // Original copyright notice:
0026
0027 // Copyright (c) Charles Karney (2008-2017) <charles@karney.com> and licensed
0028 // under the MIT/X11 License. For more information, see
0029 // https://geographiclib.sourceforge.io
0030
0031 #ifndef BOOST_GEOMETRY_FORMULAS_KARNEY_DIRECT_HPP
0032 #define BOOST_GEOMETRY_FORMULAS_KARNEY_DIRECT_HPP
0033
0034
0035 #include <boost/math/constants/constants.hpp>
0036 #include <boost/math/special_functions/hypot.hpp>
0037
0038 #include <boost/geometry/formulas/flattening.hpp>
0039 #include <boost/geometry/formulas/result_direct.hpp>
0040
0041 #include <boost/geometry/util/condition.hpp>
0042 #include <boost/geometry/util/math.hpp>
0043 #include <boost/geometry/util/normalize_spheroidal_coordinates.hpp>
0044 #include <boost/geometry/util/series_expansion.hpp>
0045
0046
0047 namespace boost { namespace geometry { namespace formula
0048 {
0049
0050 namespace se = series_expansion;
0051
0052 /*!
0053 \brief The solution of the direct problem of geodesics on latlong coordinates,
0054        after Karney (2011).
0055 \author See
0056 - Charles F.F Karney, Algorithms for geodesics, 2011
0057 https://arxiv.org/pdf/1109.4448.pdf
0058 */
0059 template <
0060     typename CT,
0061     bool EnableCoordinates = true,
0062     bool EnableReverseAzimuth = false,
0063     bool EnableReducedLength = false,
0064     bool EnableGeodesicScale = false,
0065     size_t SeriesOrder = 8
0066 >
0067 class karney_direct
0068 {
0069     static const bool CalcQuantities = EnableReducedLength || EnableGeodesicScale;
0070     static const bool CalcCoordinates = EnableCoordinates || CalcQuantities;
0071     static const bool CalcRevAzimuth = EnableReverseAzimuth || CalcCoordinates || CalcQuantities;
0072
0073 public:
0074     typedef result_direct<CT> result_type;
0075
0076     template <typename T, typename Dist, typename Azi, typename Spheroid>
0077     static inline result_type apply(T const& lo1,
0078                                     T const& la1,
0079                                     Dist const& distance,
0080                                     Azi const& azimuth12,
0081                                     Spheroid const& spheroid)
0082     {
0083         result_type result;
0084
0085         CT lon1 = lo1 * math::r2d<CT>();
0086         CT const lat1 = la1 * math::r2d<CT>();
0087
0088         Azi azi12 = azimuth12 * math::r2d<CT>();
0089         math::normalize_azimuth<degree, Azi>(azi12);
0090
0091         CT const c0 = 0;
0092         CT const c1 = 1;
0093         CT const c2 = 2;
0094
0095         CT const b = CT(get_radius<2>(spheroid));
0096         CT const f = formula::flattening<CT>(spheroid);
0097         CT const one_minus_f = c1 - f;
0098         CT const two_minus_f = c2 - f;
0099
0100         CT const n = f / two_minus_f;
0101         CT const e2 = f * two_minus_f;
0102         CT const ep2 = e2 / math::sqr(one_minus_f);
0103
0104         CT sin_alpha1, cos_alpha1;
0105         math::sin_cos_degrees<CT>(azi12, sin_alpha1, cos_alpha1);
0106
0107         // Find the reduced latitude.
0108         CT sin_beta1, cos_beta1;
0109         math::sin_cos_degrees<CT>(lat1, sin_beta1, cos_beta1);
0110         sin_beta1 *= one_minus_f;
0111
0112         math::normalize_unit_vector<CT>(sin_beta1, cos_beta1);
0113
0114         cos_beta1 = (std::max)(c0, cos_beta1);
0115
0116         // Obtain alpha 0 by solving the spherical triangle.
0117         CT const sin_alpha0 = sin_alpha1 * cos_beta1;
0118         CT const cos_alpha0 = boost::math::hypot(cos_alpha1, sin_alpha1 * sin_beta1);
0119
0120         CT const k2 = math::sqr(cos_alpha0) * ep2;
0121
0122         CT const epsilon = k2 / (c2 * (c1 + math::sqrt(c1 + k2)) + k2);
0123
0124         // Find the coefficients for A1 by computing the
0125         // series expansion using Horner scehme.
0126         CT const expansion_A1 = se::evaluate_A1<SeriesOrder>(epsilon);
0127
0128         // Index zero element of coeffs_C1 is unused.
0129         se::coeffs_C1<SeriesOrder, CT> const coeffs_C1(epsilon);
0130
0131         // Tau is an integration variable.
0132         CT const tau12 = distance / (b * (c1 + expansion_A1));
0133
0134         CT const sin_tau12 = sin(tau12);
0135         CT const cos_tau12 = cos(tau12);
0136
0137         CT sin_sigma1 = sin_beta1;
0138         CT sin_omega1 = sin_alpha0 * sin_beta1;
0139
0140         CT cos_sigma1, cos_omega1;
0141         cos_sigma1 = cos_omega1 = sin_beta1 != c0 || cos_alpha1 != c0 ? cos_beta1 * cos_alpha1 : c1;
0142         math::normalize_unit_vector<CT>(sin_sigma1, cos_sigma1);
0143
0144         CT const B11 = se::sin_cos_series(sin_sigma1, cos_sigma1, coeffs_C1);
0145         CT const sin_B11 = sin(B11);
0146         CT const cos_B11 = cos(B11);
0147
0148         CT const sin_tau1 = sin_sigma1 * cos_B11 + cos_sigma1 * sin_B11;
0149         CT const cos_tau1 = cos_sigma1 * cos_B11 - sin_sigma1 * sin_B11;
0150
0151         // Index zero element of coeffs_C1p is unused.
0152         se::coeffs_C1p<SeriesOrder, CT> const coeffs_C1p(epsilon);
0153
0154         CT const B12 = - se::sin_cos_series(sin_tau1 * cos_tau12 + cos_tau1 * sin_tau12,
0155                                             cos_tau1 * cos_tau12 - sin_tau1 * sin_tau12,
0156                                             coeffs_C1p);
0157
0158         CT const sigma12 = tau12 - (B12 - B11);
0159         CT const sin_sigma12 = sin(sigma12);
0160         CT const cos_sigma12 = cos(sigma12);
0161
0162         CT const sin_sigma2 = sin_sigma1 * cos_sigma12 + cos_sigma1 * sin_sigma12;
0163         CT const cos_sigma2 = cos_sigma1 * cos_sigma12 - sin_sigma1 * sin_sigma12;
0164
0165         if (BOOST_GEOMETRY_CONDITION(CalcRevAzimuth))
0166         {
0167             CT const sin_alpha2 = sin_alpha0;
0168             CT const cos_alpha2 = cos_alpha0 * cos_sigma2;
0169
0170             result.reverse_azimuth = atan2(sin_alpha2, cos_alpha2);
0171         }
0172
0173         if (BOOST_GEOMETRY_CONDITION(CalcCoordinates))
0174         {
0175             // Find the latitude at the second point.
0176             CT const sin_beta2 = cos_alpha0 * sin_sigma2;
0177             CT const cos_beta2 = boost::math::hypot(sin_alpha0, cos_alpha0 * cos_sigma2);
0178
0179             result.lat2 = atan2(sin_beta2, one_minus_f * cos_beta2);
0180
0181             // Find the longitude at the second point.
0182             CT const sin_omega2 = sin_alpha0 * sin_sigma2;
0183             CT const cos_omega2 = cos_sigma2;
0184
0185             CT const omega12 = atan2(sin_omega2 * cos_omega1 - cos_omega2 * sin_omega1,
0186                                      cos_omega2 * cos_omega1 + sin_omega2 * sin_omega1);
0187
0188             se::coeffs_A3<SeriesOrder, CT> const coeffs_A3(n);
0189
0190             CT const A3 = math::horner_evaluate(epsilon, coeffs_A3.begin(), coeffs_A3.end());
0191             CT const A3c = -f * sin_alpha0 * A3;
0192
0193             se::coeffs_C3<SeriesOrder, CT> const coeffs_C3(n, epsilon);
0194
0195             CT const B31 = se::sin_cos_series(sin_sigma1, cos_sigma1, coeffs_C3);
0196
0197             CT const sin_cos_res = se::sin_cos_series(sin_sigma2, cos_sigma2, coeffs_C3);
0198             CT const lam12 = omega12 + A3c * (sigma12 + (sin_cos_res - B31));
0199
0200             // Convert to degrees to get the longitudinal difference.
0201             CT lon12 = lam12 * math::r2d<CT>();
0202
0203             // Add the longitude at first point to the longitudinal
0204             // difference and normalize the result.
0205             math::normalize_longitude<degree, CT>(lon1);
0206             math::normalize_longitude<degree, CT>(lon12);
0207
0208             result.lon2 = lon1 + lon12;
0209
0210             // For longitudes close to the antimeridian the result can be out
0211             // of range. Therefore normalize.
0212             // In other formulas this has to be done at the end because
0213             // otherwise differential quantities are calculated incorrectly.
0214             // But here it's ok since result.lon2 is not used after this point.
0215             math::normalize_longitude<degree, CT>(result.lon2);
0216
0217             result.lon2 *= math::d2r<CT>();
0218         }
0219
0220         if (BOOST_GEOMETRY_CONDITION(CalcQuantities))
0221         {
0222             // Evaluate the coefficients for C2.
0223             // Index zero element of coeffs_C2 is unused.
0224             se::coeffs_C2<SeriesOrder, CT> const coeffs_C2(epsilon);
0225
0226             CT const B21 = se::sin_cos_series(sin_sigma1, cos_sigma1, coeffs_C2);
0227             CT const B22 = se::sin_cos_series(sin_sigma2, cos_sigma2, coeffs_C2);
0228
0229             // Find the coefficients for A2 by computing the
0230             // series expansion using Horner scehme.
0231             CT const expansion_A2 = se::evaluate_A2<SeriesOrder>(epsilon);
0232
0233             CT const AB1 = (c1 + expansion_A1) * (B12 - B11);
0234             CT const AB2 = (c1 + expansion_A2) * (B22 - B21);
0235             CT const J12 = (expansion_A1 - expansion_A2) * sigma12 + (AB1 - AB2);
0236
0237             CT const dn1 = math::sqrt(c1 + ep2 * math::sqr(sin_beta1));
0238             CT const dn2 = math::sqrt(c1 + k2 * math::sqr(sin_sigma2));
0239
0240             // Find the reduced length.
0241             result.reduced_length = b * ((dn2 * (cos_sigma1 * sin_sigma2) -
0242                                           dn1 * (sin_sigma1 * cos_sigma2)) -
0243                                           cos_sigma1 * cos_sigma2 * J12);
0244
0245             // Find the geodesic scale.
0246             CT const t = k2 * (sin_sigma2 - sin_sigma1) * (sin_sigma2 + sin_sigma1) / (dn1 + dn2);
0247
0248             result.geodesic_scale = cos_sigma12 + (t * sin_sigma2 - cos_sigma2 * J12) *
0249                 sin_sigma1 / dn1;
0250         }
0251
0252         return result;
0253     }
0254 };
0255
0256 }}} // namespace boost::geometry::formula
0257
0258
0259 #endif // BOOST_GEOMETRY_FORMULAS_KARNEY_DIRECT_HPP