geometry/formulas/thomas_inverse.hpp

0001 // Boost.Geometry
0002
0003 // Copyright (c) 2023 Adam Wulkiewicz, Lodz, Poland.
0004
0005 // Copyright (c) 2015-2018 Oracle and/or its affiliates.
0006 // Contributed and/or modified by Vissarion Fysikopoulos, on behalf of Oracle
0007 // Contributed and/or modified by Adam Wulkiewicz, on behalf of Oracle
0008
0009 // Use, modification and distribution is subject to the Boost Software License,
0010 // Version 1.0. (See accompanying file LICENSE_1_0.txt or copy at
0011 // http://www.boost.org/LICENSE_1_0.txt)
0012
0013 #ifndef BOOST_GEOMETRY_FORMULAS_THOMAS_INVERSE_HPP
0014 #define BOOST_GEOMETRY_FORMULAS_THOMAS_INVERSE_HPP
0015
0016
0017 #include <boost/math/constants/constants.hpp>
0018
0019 #include <boost/geometry/core/radius.hpp>
0020
0021 #include <boost/geometry/util/constexpr.hpp>
0022 #include <boost/geometry/util/math.hpp>
0023
0024 #include <boost/geometry/formulas/differential_quantities.hpp>
0025 #include <boost/geometry/formulas/flattening.hpp>
0026 #include <boost/geometry/formulas/result_inverse.hpp>
0027
0028
0029 namespace boost { namespace geometry { namespace formula
0030 {
0031
0032 /*!
0033 \brief The solution of the inverse problem of geodesics on latlong coordinates,
0034        Forsyth-Andoyer-Lambert type approximation with second order terms.
0035 \author See
0036     - Technical Report: PAUL D. THOMAS, MATHEMATICAL MODELS FOR NAVIGATION SYSTEMS, 1965
0037       http://www.dtic.mil/docs/citations/AD0627893
0038     - Technical Report: PAUL D. THOMAS, SPHEROIDAL GEODESICS, REFERENCE SYSTEMS, AND LOCAL GEOMETRY, 1970
0039       http://www.dtic.mil/docs/citations/AD0703541
0040 */
0041 template <
0042     typename CT,
0043     bool EnableDistance,
0044     bool EnableAzimuth,
0045     bool EnableReverseAzimuth = false,
0046     bool EnableReducedLength = false,
0047     bool EnableGeodesicScale = false
0048 >
0049 class thomas_inverse
0050 {
0051     static const bool CalcQuantities = EnableReducedLength || EnableGeodesicScale;
0052     static const bool CalcAzimuths = EnableAzimuth || EnableReverseAzimuth || CalcQuantities;
0053     static const bool CalcFwdAzimuth = EnableAzimuth || CalcQuantities;
0054     static const bool CalcRevAzimuth = EnableReverseAzimuth || CalcQuantities;
0055
0056 public:
0057     typedef result_inverse<CT> result_type;
0058
0059     template <typename T1, typename T2, typename Spheroid>
0060     static inline result_type apply(T1 const& lon1,
0061                                     T1 const& lat1,
0062                                     T2 const& lon2,
0063                                     T2 const& lat2,
0064                                     Spheroid const& spheroid)
0065     {
0066         result_type result;
0067
0068         // coordinates in radians
0069
0070         if ( math::equals(lon1, lon2) && math::equals(lat1, lat2) )
0071         {
0072             return result;
0073         }
0074
0075         CT const c0 = 0;
0076         CT const c1 = 1;
0077         CT const c2 = 2;
0078         CT const c4 = 4;
0079
0080         CT const pi_half = math::pi<CT>() / c2;
0081         CT const f = formula::flattening<CT>(spheroid);
0082         CT const one_minus_f = c1 - f;
0083
0084 //        CT const tan_theta1 = one_minus_f * tan(lat1);
0085 //        CT const tan_theta2 = one_minus_f * tan(lat2);
0086 //        CT const theta1 = atan(tan_theta1);
0087 //        CT const theta2 = atan(tan_theta2);
0088
0089         CT const theta1 = math::equals(lat1, pi_half) ? lat1 :
0090                           math::equals(lat1, -pi_half) ? lat1 :
0091                           atan(one_minus_f * tan(lat1));
0092         CT const theta2 = math::equals(lat2, pi_half) ? lat2 :
0093                           math::equals(lat2, -pi_half) ? lat2 :
0094                           atan(one_minus_f * tan(lat2));
0095
0096         CT const theta_m = (theta1 + theta2) / c2;
0097         CT const d_theta_m = (theta2 - theta1) / c2;
0098         CT const d_lambda = lon2 - lon1;
0099         CT const d_lambda_m = d_lambda / c2;
0100
0101         CT const sin_theta_m = sin(theta_m);
0102         CT const cos_theta_m = cos(theta_m);
0103         CT const sin_d_theta_m = sin(d_theta_m);
0104         CT const cos_d_theta_m = cos(d_theta_m);
0105         CT const sin2_theta_m = math::sqr(sin_theta_m);
0106         CT const cos2_theta_m = math::sqr(cos_theta_m);
0107         CT const sin2_d_theta_m = math::sqr(sin_d_theta_m);
0108         CT const cos2_d_theta_m = math::sqr(cos_d_theta_m);
0109         CT const sin_d_lambda_m = sin(d_lambda_m);
0110         CT const sin2_d_lambda_m = math::sqr(sin_d_lambda_m);
0111
0112         CT const H = cos2_theta_m - sin2_d_theta_m;
0113         CT const L = sin2_d_theta_m + H * sin2_d_lambda_m;
0114         CT const cos_d = c1 - c2 * L;
0115         CT const d = acos(cos_d);
0116         CT const sin_d = sin(d);
0117
0118         CT const one_minus_L = c1 - L;
0119
0120         if ( math::equals(sin_d, c0)
0121           || math::equals(L, c0)
0122           || math::equals(one_minus_L, c0) )
0123         {
0124             return result;
0125         }
0126
0127         CT const U = c2 * sin2_theta_m * cos2_d_theta_m / one_minus_L;
0128         CT const V = c2 * sin2_d_theta_m * cos2_theta_m / L;
0129         CT const X = U + V;
0130         CT const Y = U - V;
0131         CT const T = d / sin_d;
0132         CT const D = c4 * math::sqr(T);
0133         CT const E = c2 * cos_d;
0134         CT const A = D * E;
0135         CT const B = c2 * D;
0136         CT const C = T - (A - E) / c2;
0137
0138         CT const f_sqr = math::sqr(f);
0139         CT const f_sqr_per_64 = f_sqr / CT(64);
0140
0141         if BOOST_GEOMETRY_CONSTEXPR (EnableDistance)
0142         {
0143             CT const n1 = X * (A + C*X);
0144             CT const n2 = Y * (B + E*Y);
0145             CT const n3 = D*X*Y;
0146
0147             CT const delta1d = f * (T*X-Y) / c4;
0148             CT const delta2d = f_sqr_per_64 * (n1 - n2 + n3);
0149
0150             CT const a = get_radius<0>(spheroid);
0151
0152             //result.distance = a * sin_d * (T - delta1d);
0153             result.distance = a * sin_d * (T - delta1d + delta2d);
0154         }
0155
0156         if BOOST_GEOMETRY_CONSTEXPR (CalcAzimuths)
0157         {
0158             // NOTE: if both cos_latX == 0 then below we'd have 0 * INF
0159             // it's a situation when the endpoints are on the poles +-90 deg
0160             // in this case the azimuth could either be 0 or +-pi
0161             // but above always 0 is returned
0162
0163             CT const F = c2*Y-E*(c4-X);
0164             CT const M = CT(32)*T-(CT(20)*T-A)*X-(B+c4)*Y;
0165             CT const G = f*T/c2 + f_sqr_per_64 * M;
0166
0167             // TODO:
0168             // If d_lambda is close to 90 or -90 deg then tan(d_lambda) is big
0169             // and F is small. The result is not accurate.
0170             // In the edge case the result may be 2 orders of magnitude less
0171             // accurate than Andoyer's.
0172             CT const tan_d_lambda = tan(d_lambda);
0173             CT const Q = -(F*G*tan_d_lambda) / c4;
0174             CT const d_lambda_m_p = (d_lambda + Q) / c2;
0175             CT const tan_d_lambda_m_p = tan(d_lambda_m_p);
0176
0177             CT const v = atan2(cos_d_theta_m, sin_theta_m * tan_d_lambda_m_p);
0178             CT const u = atan2(-sin_d_theta_m, cos_theta_m * tan_d_lambda_m_p);
0179
0180             CT const pi = math::pi<CT>();
0181
0182             if BOOST_GEOMETRY_CONSTEXPR (CalcFwdAzimuth)
0183             {
0184                 CT alpha1 = v + u;
0185                 if (alpha1 > pi)
0186                 {
0187                     alpha1 -= c2 * pi;
0188                 }
0189
0190                 result.azimuth = alpha1;
0191             }
0192
0193             if BOOST_GEOMETRY_CONSTEXPR (CalcRevAzimuth)
0194             {
0195                 CT alpha2 = pi - (v - u);
0196                 if (alpha2 > pi)
0197                 {
0198                     alpha2 -= c2 * pi;
0199                 }
0200
0201                 result.reverse_azimuth = alpha2;
0202             }
0203         }
0204
0205         if BOOST_GEOMETRY_CONSTEXPR (CalcQuantities)
0206         {
0207             typedef differential_quantities<CT, EnableReducedLength, EnableGeodesicScale, 2> quantities;
0208             quantities::apply(lon1, lat1, lon2, lat2,
0209                               result.azimuth, result.reverse_azimuth,
0210                               get_radius<2>(spheroid), f,
0211                               result.reduced_length, result.geodesic_scale);
0212         }
0213
0214         return result;
0215     }
0216 };
0217
0218 }}} // namespace boost::geometry::formula
0219
0220
0221 #endif // BOOST_GEOMETRY_FORMULAS_THOMAS_INVERSE_HPP