summaryrefslogtreecommitdiff log msg author committer range
 ```1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177 178 179 180 181 182 183 184 ``` ``````// Boost.Geometry // 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. // Contributed and/or modified by Adam Wulkiewicz, on behalf of Oracle // Use, modification and distribution is subject to the Boost Software License, // Version 1.0. (See accompanying file LICENSE_1_0.txt or copy at // http://www.boost.org/LICENSE_1_0.txt) #ifndef BOOST_GEOMETRY_FORMULAS_VINCENTY_DIRECT_HPP #define BOOST_GEOMETRY_FORMULAS_VINCENTY_DIRECT_HPP #include #include #include #include #include #include #include #include #ifndef BOOST_GEOMETRY_DETAIL_VINCENTY_MAX_STEPS #define BOOST_GEOMETRY_DETAIL_VINCENTY_MAX_STEPS 1000 #endif namespace boost { namespace geometry { namespace formula { /*! \brief The solution of the direct 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 \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 - http://futureboy.homeip.net/fsp/colorize.fsp?fileName=navigation.frink */ template < typename CT, bool EnableCoordinates = true, bool EnableReverseAzimuth = false, bool EnableReducedLength = false, bool EnableGeodesicScale = false > class vincenty_direct { static const bool CalcQuantities = EnableReducedLength || EnableGeodesicScale; static const bool CalcCoordinates = EnableCoordinates || CalcQuantities; static const bool CalcRevAzimuth = EnableReverseAzimuth || CalcQuantities; public: typedef result_direct result_type; template static inline result_type apply(T const& lo1, T const& la1, Dist const& distance, Azi const& azimuth12, Spheroid const& spheroid) { result_type result; CT const lon1 = lo1; CT const lat1 = la1; if ( math::equals(distance, Dist(0)) || distance < Dist(0) ) { result.lon2 = lon1; result.lat2 = lat1; return result; } CT const radius_a = CT(get_radius<0>(spheroid)); CT const radius_b = CT(get_radius<2>(spheroid)); CT const flattening = formula::flattening(spheroid); CT const sin_azimuth12 = sin(azimuth12); CT const cos_azimuth12 = cos(azimuth12); // U: reduced latitude, defined by tan U = (1-f) tan phi CT const one_min_f = CT(1) - flattening; CT const tan_U1 = one_min_f * tan(lat1); CT const sigma1 = atan2(tan_U1, cos_azimuth12); // (1) // may be calculated from tan using 1 sqrt() CT const U1 = atan(tan_U1); CT const sin_U1 = sin(U1); CT const cos_U1 = cos(U1); CT const sin_alpha = cos_U1 * sin_azimuth12; // (2) CT const sin_alpha_sqr = math::sqr(sin_alpha); CT const cos_alpha_sqr = CT(1) - sin_alpha_sqr; CT const b_sqr = radius_b * radius_b; CT const u_sqr = cos_alpha_sqr * (radius_a * radius_a - b_sqr) / b_sqr; CT const A = CT(1) + (u_sqr/CT(16384)) * (CT(4096) + u_sqr*(CT(-768) + u_sqr*(CT(320) - u_sqr*CT(175)))); // (3) CT const B = (u_sqr/CT(1024))*(CT(256) + u_sqr*(CT(-128) + u_sqr*(CT(74) - u_sqr*CT(47)))); // (4) CT s_div_bA = distance / (radius_b * A); CT sigma = s_div_bA; // (7) CT previous_sigma; CT sin_sigma; CT cos_sigma; CT cos_2sigma_m; CT cos_2sigma_m_sqr; int counter = 0; // robustness do { previous_sigma = sigma; CT const two_sigma_m = CT(2) * sigma1 + sigma; // (5) sin_sigma = sin(sigma); cos_sigma = cos(sigma); CT const sin_sigma_sqr = math::sqr(sin_sigma); cos_2sigma_m = cos(two_sigma_m); cos_2sigma_m_sqr = math::sqr(cos_2sigma_m); CT const delta_sigma = B * sin_sigma * (cos_2sigma_m + (B/CT(4)) * ( cos_sigma * (CT(-1) + CT(2)*cos_2sigma_m_sqr) - (B/CT(6) * cos_2sigma_m * (CT(-3)+CT(4)*sin_sigma_sqr) * (CT(-3)+CT(4)*cos_2sigma_m_sqr)) )); // (6) sigma = s_div_bA + delta_sigma; // (7) ++counter; // robustness } while ( geometry::math::abs(previous_sigma - sigma) > CT(1e-12) //&& geometry::math::abs(sigma) < pi && counter < BOOST_GEOMETRY_DETAIL_VINCENTY_MAX_STEPS ); // robustness if (BOOST_GEOMETRY_CONDITION(CalcCoordinates)) { result.lat2 = atan2( sin_U1 * cos_sigma + cos_U1 * sin_sigma * cos_azimuth12, one_min_f * math::sqrt(sin_alpha_sqr + math::sqr(sin_U1 * sin_sigma - cos_U1 * cos_sigma * cos_azimuth12))); // (8) CT const lambda = atan2( sin_sigma * sin_azimuth12, cos_U1 * cos_sigma - sin_U1 * sin_sigma * cos_azimuth12); // (9) CT const C = (flattening/CT(16)) * cos_alpha_sqr * ( CT(4) + flattening * ( CT(4) - CT(3) * cos_alpha_sqr ) ); // (10) CT const L = lambda - (CT(1) - C) * flattening * sin_alpha * ( sigma + C * sin_sigma * ( cos_2sigma_m + C * cos_sigma * ( CT(-1) + CT(2) * cos_2sigma_m_sqr ) ) ); // (11) result.lon2 = lon1 + L; } if (BOOST_GEOMETRY_CONDITION(CalcRevAzimuth)) { result.reverse_azimuth = atan2(sin_alpha, -sin_U1 * sin_sigma + cos_U1 * cos_sigma * cos_azimuth12); // (12) } if (BOOST_GEOMETRY_CONDITION(CalcQuantities)) { typedef differential_quantities quantities; quantities::apply(lon1, lat1, result.lon2, result.lat2, azimuth12, result.reverse_azimuth, radius_b, flattening, result.reduced_length, result.geodesic_scale); } return result; } }; }}} // namespace boost::geometry::formula #endif // BOOST_GEOMETRY_FORMULAS_VINCENTY_DIRECT_HPP ``````