diff options
Diffstat (limited to 'boost/geometry/algorithms/detail/andoyer_inverse.hpp')
| -rw-r--r-- | boost/geometry/algorithms/detail/andoyer_inverse.hpp | 156 |
1 files changed, 99 insertions, 57 deletions
diff --git a/boost/geometry/algorithms/detail/andoyer_inverse.hpp b/boost/geometry/algorithms/detail/andoyer_inverse.hpp index eb8485b8e1..66ad6446e9 100644 --- a/boost/geometry/algorithms/detail/andoyer_inverse.hpp +++ b/boost/geometry/algorithms/detail/andoyer_inverse.hpp @@ -1,6 +1,6 @@ // Boost.Geometry -// Copyright (c) 2015 Oracle and/or its affiliates. +// Copyright (c) 2015-2016 Oracle and/or its affiliates. // Contributed and/or modified by Adam Wulkiewicz, on behalf of Oracle @@ -49,23 +49,17 @@ struct andoyer_inverse T2 const& lat2, Spheroid const& spheroid) { + CT const c0 = CT(0); + CT const c1 = CT(1); + CT const pi = math::pi<CT>(); + result_type result; // coordinates in radians - if ( math::equals(lon1, lon2) - && math::equals(lat1, lat2) ) + if ( math::equals(lon1, lon2) && math::equals(lat1, lat2) ) { - result.set(CT(0), CT(0)); - return result; - } - - CT const pi_half = math::pi<CT>() / CT(2); - - if ( math::equals(math::abs(lat1), pi_half) - && math::equals(math::abs(lat2), pi_half) ) - { - result.set(CT(0), CT(0)); + result.set(c0, c0); return result; } @@ -80,22 +74,16 @@ struct andoyer_inverse // H,G,T = infinity if cos_d = 1 or cos_d = -1 // lat1 == +-90 && lat2 == +-90 // lat1 == lat2 && lon1 == lon2 - CT const cos_d = sin_lat1*sin_lat2 + cos_lat1*cos_lat2*cos_dlon; - CT const d = acos(cos_d); - CT const sin_d = sin(d); - - // just in case since above lat1 and lat2 is checked - // the check below is equal to cos_d == 1 || cos_d == -1 || d == 0 - if ( math::equals(sin_d, CT(0)) ) - { - result.set(CT(0), CT(0)); - return result; - } - - // if the function returned before this place - // and endpoints were on the poles +-90 deg - // in this case the azimuth could either be 0 or +-pi - + CT cos_d = sin_lat1*sin_lat2 + cos_lat1*cos_lat2*cos_dlon; + // on some platforms cos_d may be outside valid range + if (cos_d < -c1) + cos_d = -c1; + else if (cos_d > c1) + cos_d = c1; + + CT const d = acos(cos_d); // [0, pi] + CT const sin_d = sin(d); // [-1, 1] + CT const f = detail::flattening<CT>(spheroid); if ( BOOST_GEOMETRY_CONDITION(EnableDistance) ) @@ -103,11 +91,18 @@ struct andoyer_inverse CT const K = math::sqr(sin_lat1-sin_lat2); CT const L = math::sqr(sin_lat1+sin_lat2); CT const three_sin_d = CT(3) * sin_d; - // H or G = infinity if cos_d = 1 or cos_d = -1 - CT const H = (d+three_sin_d)/(CT(1)-cos_d); - CT const G = (d-three_sin_d)/(CT(1)+cos_d); - // for e.g. lat1=-90 && lat2=90 here we have G*L=INF*0 + CT const one_minus_cos_d = c1 - cos_d; + CT const one_plus_cos_d = c1 + cos_d; + // cos_d = 1 or cos_d = -1 means that the points are antipodal + + CT const H = math::equals(one_minus_cos_d, c0) ? + c0 : + (d + three_sin_d) / one_minus_cos_d; + CT const G = math::equals(one_plus_cos_d, c0) ? + c0 : + (d - three_sin_d) / one_plus_cos_d; + CT const dd = -(f/CT(4))*(H*K+G*L); CT const a = get_radius<0>(spheroid); @@ -116,41 +111,88 @@ struct andoyer_inverse } else { - result.distance = CT(0); + result.distance = c0; } if ( BOOST_GEOMETRY_CONDITION(EnableAzimuth) ) { - CT A = CT(0); - CT U = CT(0); - if ( ! math::equals(cos_lat2, CT(0)) ) + // sin_d = 0 <=> antipodal points + if (math::equals(sin_d, c0)) { - CT const tan_lat2 = sin_lat2/cos_lat2; - CT const M = cos_lat1*tan_lat2-sin_lat1*cos_dlon; - A = atan2(sin_dlon, M); - CT const sin_2A = sin(CT(2)*A); - U = (f/CT(2))*math::sqr(cos_lat1)*sin_2A; + // T = inf + // dA = inf + // azimuth = -inf + if (lat1 <= lat2) + result.azimuth = c0; + else + result.azimuth = pi; } - - CT V = CT(0); - if ( ! math::equals(cos_lat1, CT(0)) ) + else { - CT const tan_lat1 = sin_lat1/cos_lat1; - CT const N = cos_lat2*tan_lat1-sin_lat2*cos_dlon; - CT const B = atan2(sin_dlon, N); - CT const sin_2B = sin(CT(2)*B); - V = (f/CT(2))*math::sqr(cos_lat2)*sin_2B; + CT const c2 = CT(2); + + CT A = c0; + CT U = c0; + if ( ! math::equals(cos_lat2, c0) ) + { + CT const tan_lat2 = sin_lat2/cos_lat2; + CT const M = cos_lat1*tan_lat2-sin_lat1*cos_dlon; + A = atan2(sin_dlon, M); + CT const sin_2A = sin(c2*A); + U = (f/ c2)*math::sqr(cos_lat1)*sin_2A; + } + + CT V = c0; + if ( ! math::equals(cos_lat1, c0) ) + { + CT const tan_lat1 = sin_lat1/cos_lat1; + CT const N = cos_lat2*tan_lat1-sin_lat2*cos_dlon; + CT const B = atan2(sin_dlon, N); + CT const sin_2B = sin(c2*B); + V = (f/ c2)*math::sqr(cos_lat2)*sin_2B; + } + + CT const T = d / sin_d; + CT const dA = V*T-U; + + result.azimuth = A - dA; + + // even with sin_d == 0 checked above if the second point + // is somewhere in the antipodal area T may still be great + // therefore dA may be great and the resulting azimuth + // may be some more or less arbitrary angle + if (A >= c0) // A indicates Eastern hemisphere + { + if (dA >= c0) // A altered towards 0 + { + if ((result.azimuth) < c0) + result.azimuth = c0; + } + else // dA < 0, A altered towards pi + { + if (result.azimuth > pi) + result.azimuth = pi; + } + } + else // A indicates Western hemisphere + { + if (dA <= c0) // A altered towards 0 + { + if (result.azimuth > c0) + result.azimuth = c0; + } + else // dA > 0, A altered towards -pi + { + CT const minus_pi = -pi; + if ((result.azimuth) < minus_pi) + result.azimuth = minus_pi; + } + } } - - // infinity if sin_d = 0, so cos_d = 1 or cos_d = -1 - CT const T = d / sin_d; - CT const dA = V*T-U; - - result.azimuth = A - dA; } else { - result.azimuth = CT(0); + result.azimuth = c0; } return result; |