ceph/src/boost/boost/geometry/formulas/vincenty_direct.hpp

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