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_ALGORITHMS_DETAIL_ANDOYER_INVERSE_HPP
12 #define BOOST_GEOMETRY_ALGORITHMS_DETAIL_ANDOYER_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/algorithms/detail/flattening.hpp>
24 #include <boost/geometry/algorithms/detail/result_inverse.hpp>
27 namespace boost { namespace geometry { namespace detail
31 \brief The solution of the inverse problem of geodesics on latlong coordinates,
32 Forsyth-Andoyer-Lambert type approximation with first order terms.
34 - Technical Report: PAUL D. THOMAS, MATHEMATICAL MODELS FOR NAVIGATION SYSTEMS, 1965
35 http://www.dtic.mil/docs/citations/AD0627893
36 - Technical Report: PAUL D. THOMAS, SPHEROIDAL GEODESICS, REFERENCE SYSTEMS, AND LOCAL GEOMETRY, 1970
37 http://www.dtic.mil/docs/citations/AD703541
40 template <typename CT, bool EnableDistance, bool EnableAzimuth>
41 struct andoyer_inverse
43 typedef result_inverse<CT> result_type;
45 template <typename T1, typename T2, typename Spheroid>
46 static inline result_type apply(T1 const& lon1,
50 Spheroid const& spheroid)
54 CT const pi = math::pi<CT>();
58 // coordinates in radians
60 if ( math::equals(lon1, lon2) && math::equals(lat1, lat2) )
66 CT const dlon = lon2 - lon1;
67 CT const sin_dlon = sin(dlon);
68 CT const cos_dlon = cos(dlon);
69 CT const sin_lat1 = sin(lat1);
70 CT const cos_lat1 = cos(lat1);
71 CT const sin_lat2 = sin(lat2);
72 CT const cos_lat2 = cos(lat2);
74 // H,G,T = infinity if cos_d = 1 or cos_d = -1
75 // lat1 == +-90 && lat2 == +-90
76 // lat1 == lat2 && lon1 == lon2
77 CT cos_d = sin_lat1*sin_lat2 + cos_lat1*cos_lat2*cos_dlon;
78 // on some platforms cos_d may be outside valid range
84 CT const d = acos(cos_d); // [0, pi]
85 CT const sin_d = sin(d); // [-1, 1]
87 CT const f = detail::flattening<CT>(spheroid);
89 if ( BOOST_GEOMETRY_CONDITION(EnableDistance) )
91 CT const K = math::sqr(sin_lat1-sin_lat2);
92 CT const L = math::sqr(sin_lat1+sin_lat2);
93 CT const three_sin_d = CT(3) * sin_d;
95 CT const one_minus_cos_d = c1 - cos_d;
96 CT const one_plus_cos_d = c1 + cos_d;
97 // cos_d = 1 or cos_d = -1 means that the points are antipodal
99 CT const H = math::equals(one_minus_cos_d, c0) ?
101 (d + three_sin_d) / one_minus_cos_d;
102 CT const G = math::equals(one_plus_cos_d, c0) ?
104 (d - three_sin_d) / one_plus_cos_d;
106 CT const dd = -(f/CT(4))*(H*K+G*L);
108 CT const a = get_radius<0>(spheroid);
110 result.distance = a * (d + dd);
114 result.distance = c0;
117 if ( BOOST_GEOMETRY_CONDITION(EnableAzimuth) )
119 // sin_d = 0 <=> antipodal points
120 if (math::equals(sin_d, c0))
136 if ( ! math::equals(cos_lat2, c0) )
138 CT const tan_lat2 = sin_lat2/cos_lat2;
139 CT const M = cos_lat1*tan_lat2-sin_lat1*cos_dlon;
140 A = atan2(sin_dlon, M);
141 CT const sin_2A = sin(c2*A);
142 U = (f/ c2)*math::sqr(cos_lat1)*sin_2A;
146 if ( ! math::equals(cos_lat1, c0) )
148 CT const tan_lat1 = sin_lat1/cos_lat1;
149 CT const N = cos_lat2*tan_lat1-sin_lat2*cos_dlon;
150 CT const B = atan2(sin_dlon, N);
151 CT const sin_2B = sin(c2*B);
152 V = (f/ c2)*math::sqr(cos_lat2)*sin_2B;
155 CT const T = d / sin_d;
158 result.azimuth = A - dA;
160 // even with sin_d == 0 checked above if the second point
161 // is somewhere in the antipodal area T may still be great
162 // therefore dA may be great and the resulting azimuth
163 // may be some more or less arbitrary angle
164 if (A >= c0) // A indicates Eastern hemisphere
166 if (dA >= c0) // A altered towards 0
168 if ((result.azimuth) < c0)
171 else // dA < 0, A altered towards pi
173 if (result.azimuth > pi)
177 else // A indicates Western hemisphere
179 if (dA <= c0) // A altered towards 0
181 if (result.azimuth > c0)
184 else // dA > 0, A altered towards -pi
186 CT const minus_pi = -pi;
187 if ((result.azimuth) < minus_pi)
188 result.azimuth = minus_pi;
202 }}} // namespace boost::geometry::detail
205 #endif // BOOST_GEOMETRY_ALGORITHMS_DETAIL_ANDOYER_INVERSE_HPP