geometry/formulas/thomas_inverse.hpp

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