3 // Copyright (c) 2015-2016 Oracle and/or its affiliates.
5 // Contributed and/or modified by Adam Wulkiewicz, on behalf of Oracle
7 // Use, modification and distribution is subject to the Boost Software License,
8 // Version 1.0. (See accompanying file LICENSE_1_0.txt or copy at
9 // http://www.boost.org/LICENSE_1_0.txt)
11 #ifndef BOOST_GEOMETRY_FORMULAS_THOMAS_INVERSE_HPP
12 #define BOOST_GEOMETRY_FORMULAS_THOMAS_INVERSE_HPP
15 #include <boost/math/constants/constants.hpp>
17 #include <boost/geometry/core/radius.hpp>
18 #include <boost/geometry/core/srs.hpp>
20 #include <boost/geometry/util/condition.hpp>
21 #include <boost/geometry/util/math.hpp>
23 #include <boost/geometry/formulas/differential_quantities.hpp>
24 #include <boost/geometry/formulas/flattening.hpp>
25 #include <boost/geometry/formulas/result_inverse.hpp>
28 namespace boost { namespace geometry { namespace formula
32 \brief The solution of the inverse problem of geodesics on latlong coordinates,
33 Forsyth-Andoyer-Lambert type approximation with second order terms.
35 - Technical Report: PAUL D. THOMAS, MATHEMATICAL MODELS FOR NAVIGATION SYSTEMS, 1965
36 http://www.dtic.mil/docs/citations/AD0627893
37 - Technical Report: PAUL D. THOMAS, SPHEROIDAL GEODESICS, REFERENCE SYSTEMS, AND LOCAL GEOMETRY, 1970
38 http://www.dtic.mil/docs/citations/AD0703541
44 bool EnableReverseAzimuth = false,
45 bool EnableReducedLength = false,
46 bool EnableGeodesicScale = false
50 static const bool CalcQuantities = EnableReducedLength || EnableGeodesicScale;
51 static const bool CalcAzimuths = EnableAzimuth || EnableReverseAzimuth || CalcQuantities;
52 static const bool CalcFwdAzimuth = EnableAzimuth || CalcQuantities;
53 static const bool CalcRevAzimuth = EnableReverseAzimuth || CalcQuantities;
56 typedef result_inverse<CT> result_type;
58 template <typename T1, typename T2, typename Spheroid>
59 static inline result_type apply(T1 const& lon1,
63 Spheroid const& spheroid)
67 // coordinates in radians
69 if ( math::equals(lon1, lon2) && math::equals(lat1, lat2) )
79 CT const pi_half = math::pi<CT>() / c2;
80 CT const f = formula::flattening<CT>(spheroid);
81 CT const one_minus_f = c1 - f;
83 // CT const tan_theta1 = one_minus_f * tan(lat1);
84 // CT const tan_theta2 = one_minus_f * tan(lat2);
85 // CT const theta1 = atan(tan_theta1);
86 // CT const theta2 = atan(tan_theta2);
88 CT const theta1 = math::equals(lat1, pi_half) ? lat1 :
89 math::equals(lat1, -pi_half) ? lat1 :
90 atan(one_minus_f * tan(lat1));
91 CT const theta2 = math::equals(lat2, pi_half) ? lat2 :
92 math::equals(lat2, -pi_half) ? lat2 :
93 atan(one_minus_f * tan(lat2));
95 CT const theta_m = (theta1 + theta2) / c2;
96 CT const d_theta_m = (theta2 - theta1) / c2;
97 CT const d_lambda = lon2 - lon1;
98 CT const d_lambda_m = d_lambda / c2;
100 CT const sin_theta_m = sin(theta_m);
101 CT const cos_theta_m = cos(theta_m);
102 CT const sin_d_theta_m = sin(d_theta_m);
103 CT const cos_d_theta_m = cos(d_theta_m);
104 CT const sin2_theta_m = math::sqr(sin_theta_m);
105 CT const cos2_theta_m = math::sqr(cos_theta_m);
106 CT const sin2_d_theta_m = math::sqr(sin_d_theta_m);
107 CT const cos2_d_theta_m = math::sqr(cos_d_theta_m);
108 CT const sin_d_lambda_m = sin(d_lambda_m);
109 CT const sin2_d_lambda_m = math::sqr(sin_d_lambda_m);
111 CT const H = cos2_theta_m - sin2_d_theta_m;
112 CT const L = sin2_d_theta_m + H * sin2_d_lambda_m;
113 CT const cos_d = c1 - c2 * L;
114 CT const d = acos(cos_d);
115 CT const sin_d = sin(d);
117 CT const one_minus_L = c1 - L;
119 if ( math::equals(sin_d, c0)
120 || math::equals(L, c0)
121 || math::equals(one_minus_L, c0) )
126 CT const U = c2 * sin2_theta_m * cos2_d_theta_m / one_minus_L;
127 CT const V = c2 * sin2_d_theta_m * cos2_theta_m / L;
130 CT const T = d / sin_d;
131 CT const D = c4 * math::sqr(T);
132 CT const E = c2 * cos_d;
135 CT const C = T - (A - E) / c2;
137 CT const f_sqr = math::sqr(f);
138 CT const f_sqr_per_64 = f_sqr / CT(64);
140 if ( BOOST_GEOMETRY_CONDITION(EnableDistance) )
142 CT const n1 = X * (A + C*X);
143 CT const n2 = Y * (B + E*Y);
146 CT const delta1d = f * (T*X-Y) / c4;
147 CT const delta2d = f_sqr_per_64 * (n1 - n2 + n3);
149 CT const a = get_radius<0>(spheroid);
151 //result.distance = a * sin_d * (T - delta1d);
152 result.distance = a * sin_d * (T - delta1d + delta2d);
155 if ( BOOST_GEOMETRY_CONDITION(CalcAzimuths) )
157 // NOTE: if both cos_latX == 0 then below we'd have 0 * INF
158 // it's a situation when the endpoints are on the poles +-90 deg
159 // in this case the azimuth could either be 0 or +-pi
160 // but above always 0 is returned
162 CT const F = c2*Y-E*(c4-X);
163 CT const M = CT(32)*T-(CT(20)*T-A)*X-(B+c4)*Y;
164 CT const G = f*T/c2 + f_sqr_per_64 * M;
167 // If d_lambda is close to 90 or -90 deg then tan(d_lambda) is big
168 // and F is small. The result is not accurate.
169 // In the edge case the result may be 2 orders of magnitude less
170 // accurate than Andoyer's.
171 CT const tan_d_lambda = tan(d_lambda);
172 CT const Q = -(F*G*tan_d_lambda) / c4;
173 CT const d_lambda_m_p = (d_lambda + Q) / c2;
174 CT const tan_d_lambda_m_p = tan(d_lambda_m_p);
176 CT const v = atan2(cos_d_theta_m, sin_theta_m * tan_d_lambda_m_p);
177 CT const u = atan2(-sin_d_theta_m, cos_theta_m * tan_d_lambda_m_p);
179 CT const pi = math::pi<CT>();
181 if (BOOST_GEOMETRY_CONDITION(EnableAzimuth))
189 result.azimuth = alpha1;
192 if (BOOST_GEOMETRY_CONDITION(EnableReverseAzimuth))
194 CT alpha2 = pi - (v - u);
200 result.reverse_azimuth = alpha2;
204 if (BOOST_GEOMETRY_CONDITION(CalcQuantities))
206 typedef differential_quantities<CT, EnableReducedLength, EnableGeodesicScale, 2> quantities;
207 quantities::apply(lon1, lat1, lon2, lat2,
208 result.azimuth, result.reverse_azimuth,
209 get_radius<2>(spheroid), f,
210 result.reduced_length, result.geodesic_scale);
217 }}} // namespace boost::geometry::formula
220 #endif // BOOST_GEOMETRY_FORMULAS_THOMAS_INVERSE_HPP