diff options
Diffstat (limited to 'boost/geometry/formulas/vincenty_inverse.hpp')
| -rw-r--r-- | boost/geometry/formulas/vincenty_inverse.hpp | 38 |
1 files changed, 21 insertions, 17 deletions
diff --git a/boost/geometry/formulas/vincenty_inverse.hpp b/boost/geometry/formulas/vincenty_inverse.hpp index bbda00036b..032e16e291 100644 --- a/boost/geometry/formulas/vincenty_inverse.hpp +++ b/boost/geometry/formulas/vincenty_inverse.hpp @@ -2,8 +2,8 @@ // Copyright (c) 2007-2012 Barend Gehrels, Amsterdam, the Netherlands. -// This file was modified by Oracle on 2014, 2016. -// Modifications copyright (c) 2014-2016 Oracle and/or its affiliates. +// This file was modified by Oracle on 2014, 2016, 2017. +// Modifications copyright (c) 2014-2017 Oracle and/or its affiliates. // Contributed and/or modified by Adam Wulkiewicz, on behalf of Oracle @@ -23,9 +23,8 @@ #include <boost/geometry/util/condition.hpp> #include <boost/geometry/util/math.hpp> -#include <boost/geometry/algorithms/detail/flattening.hpp> - #include <boost/geometry/formulas/differential_quantities.hpp> +#include <boost/geometry/formulas/flattening.hpp> #include <boost/geometry/formulas/result_inverse.hpp> @@ -41,7 +40,7 @@ namespace boost { namespace geometry { namespace formula \brief The solution of the inverse problem of geodesics on latlong coordinates, after Vincenty, 1975 \author See - http://www.ngs.noaa.gov/PUBS_LIB/inverse.pdf - - http://www.icsm.gov.au/gda/gdav2.3.pdf + - http://www.icsm.gov.au/gda/gda-v_2.4.pdf \author Adapted from various implementations to get it close to the original document - http://www.movable-type.co.uk/scripts/LatLongVincenty.html - http://exogen.case.edu/projects/geopy/source/geopy.distance.html @@ -99,10 +98,10 @@ public: CT const radius_a = CT(get_radius<0>(spheroid)); CT const radius_b = CT(get_radius<2>(spheroid)); - CT const flattening = geometry::detail::flattening<CT>(spheroid); + CT const f = formula::flattening<CT>(spheroid); // U: reduced latitude, defined by tan U = (1-f) tan phi - CT const one_min_f = c1 - flattening; + CT const one_min_f = c1 - f; CT const tan_U1 = one_min_f * tan(lat1); // above (1) CT const tan_U2 = one_min_f * tan(lat2); // above (1) @@ -113,8 +112,9 @@ public: CT const cos_U1 = c1 / temp_den_U1; CT const cos_U2 = c1 / temp_den_U2; // sin = tan / sqrt(1 + tan^2) - CT const sin_U1 = tan_U1 / temp_den_U1; - CT const sin_U2 = tan_U2 / temp_den_U2; + // sin = tan * cos + CT const sin_U1 = tan_U1 * cos_U1; + CT const sin_U2 = tan_U2 * cos_U2; // calculate sin U and cos U directly //CT const U1 = atan(tan_U1); @@ -130,7 +130,8 @@ public: CT sin_sigma; CT sin_alpha; CT cos2_alpha; - CT cos2_sigma_m; + CT cos_2sigma_m; + CT cos2_2sigma_m; CT sigma; int counter = 0; // robustness @@ -144,12 +145,13 @@ public: CT cos_sigma = sin_U1 * sin_U2 + cos_U1 * cos_U2 * cos_lambda; // (15) sin_alpha = cos_U1 * cos_U2 * sin_lambda / sin_sigma; // (17) cos2_alpha = c1 - math::sqr(sin_alpha); - cos2_sigma_m = math::equals(cos2_alpha, 0) ? 0 : cos_sigma - c2 * sin_U1 * sin_U2 / cos2_alpha; // (18) + cos_2sigma_m = math::equals(cos2_alpha, 0) ? 0 : cos_sigma - c2 * sin_U1 * sin_U2 / cos2_alpha; // (18) + cos2_2sigma_m = math::sqr(cos_2sigma_m); - CT C = flattening/c16 * cos2_alpha * (c4 + flattening * (c4 - c3 * cos2_alpha)); // (10) + CT C = f/c16 * cos2_alpha * (c4 + f * (c4 - c3 * cos2_alpha)); // (10) sigma = atan2(sin_sigma, cos_sigma); // (16) - lambda = L + (c1 - C) * flattening * sin_alpha * - (sigma + C * sin_sigma * ( cos2_sigma_m + C * cos_sigma * (-c1 + c2 * math::sqr(cos2_sigma_m)))); // (11) + lambda = L + (c1 - C) * f * sin_alpha * + (sigma + C * sin_sigma * (cos_2sigma_m + C * cos_sigma * (-c1 + c2 * cos2_2sigma_m))); // (11) ++counter; // robustness @@ -182,8 +184,10 @@ public: CT A = c1 + sqr_u/c16384 * (c4096 + sqr_u * (-c768 + sqr_u * (c320 - c175 * sqr_u))); // (3) CT B = sqr_u/c1024 * (c256 + sqr_u * ( -c128 + sqr_u * (c74 - c47 * sqr_u))); // (4) - CT delta_sigma = B * sin_sigma * ( cos2_sigma_m + (B/c4) * (cos(sigma)* (-c1 + c2 * cos2_sigma_m) - - (B/c6) * cos2_sigma_m * (-c3 + c4 * math::sqr(sin_sigma)) * (-c3 + c4 * cos2_sigma_m))); // (6) + CT const cos_sigma = cos(sigma); + CT const sin2_sigma = math::sqr(sin_sigma); + CT delta_sigma = B * sin_sigma * (cos_2sigma_m + (B/c4) * (cos_sigma* (-c1 + c2 * cos2_2sigma_m) + - (B/c6) * cos_2sigma_m * (-c3 + c4 * sin2_sigma) * (-c3 + c4 * cos2_2sigma_m))); // (6) result.distance = radius_b * A * (sigma - delta_sigma); // (19) } @@ -206,7 +210,7 @@ public: typedef differential_quantities<CT, EnableReducedLength, EnableGeodesicScale, 2> quantities; quantities::apply(lon1, lat1, lon2, lat2, result.azimuth, result.reverse_azimuth, - radius_b, flattening, + radius_b, f, result.reduced_length, result.geodesic_scale); } |