Changeset 63839


Ignore:
Timestamp:
Jul 11, 2010, 9:40:32 AM (8 years ago)
Author:
Barend Gehrels
Message:

Updated definition of PI to support templated UDT
Updated Andoyer for high precision

Location:
sandbox/geometry/boost/geometry
Files:
4 edited

Legend:

Unmodified
Added
Removed
  • sandbox/geometry/boost/geometry/extensions/contrib/ttmath_stub.hpp

    r63707 r63839  
    33
    44#include <boost/math/constants/constants.hpp>
     5#include <boost/geometry/util/math.hpp>
    56
    67
     
    108109};
    109110
    110 namespace boost{ namespace math { namespace constants
     111namespace boost{ namespace geometry { namespace math
     112{
     113
     114namespace detail
    111115{
    112116    // Workaround for boost::math::constants::pi:
     
    114118    // 2) because it is implemented as a function, generic implementation not possible
    115119
     120    // Partial specialization for ttmath
    116121    template <ttmath::uint Exponent, ttmath::uint Mantissa>
    117     inline ttmath::Big<Exponent, Mantissa> ttmath_pi()
     122    struct define_pi<ttmath::Big<Exponent, Mantissa> >
    118123    {
    119         static ttmath::Big<Exponent, Mantissa> const the_pi("3.14159265358979323846264338327950288419716939937510582097494459230781640628620899862803482534211706798214808651328230664709384460955058223172535940812848111745028410270193852110555964462294895493038196");
    120         return the_pi;
    121     }
     124        static inline ttmath::Big<Exponent, Mantissa> apply()
     125        {
     126            static ttmath::Big<Exponent, Mantissa> const the_pi(
     127                "3.14159265358979323846264338327950288419716939937510582097494459230781640628620899862803482534211706798214808651328230664709384460955058223172535940812848111745028410270193852110555964462294895493038196");
     128            return the_pi;
     129        }
     130    };
    122131
     132    template <>
     133    struct define_pi<ttmath_big>
     134    {
     135        static inline ttmath_big apply()
     136        {
     137            return define_pi<ttmath::Big<1,4> >::apply();
     138        }
     139    };
    123140
    124     template <> inline ttmath::Big<1,4> pi<ttmath::Big<1,4> >(BOOST_MATH_EXPLICIT_TEMPLATE_TYPE_SPEC( (ttmath::Big<1,4>) ))
    125     {
    126         return ttmath_pi<1,4>();
    127     }
    128     template <> inline ttmath::Big<1,8> pi<ttmath::Big<1,8> >(BOOST_MATH_EXPLICIT_TEMPLATE_TYPE_SPEC( (ttmath::Big<1,8>) ))
    129     {
    130         return ttmath_pi<1,8>();
    131     }
     141}
     142}}} // boost::geometry::math
    132143
    133     // and so on...
    134     // or fix lexical_cast but probably has the same problem
    135 
    136 
    137     template <> inline ttmath_big pi<ttmath_big >(BOOST_MATH_EXPLICIT_TEMPLATE_TYPE_SPEC(ttmath_big))
    138     {
    139         return ttmath_pi<1,4>();
    140     }
    141 
    142 }}}
    143144
    144145
  • sandbox/geometry/boost/geometry/extensions/gis/geographic/strategies/andoyer.hpp

    r63707 r63839  
    1414#include <boost/geometry/core/radian_access.hpp>
    1515#include <boost/geometry/core/coordinate_type.hpp>
     16#include <boost/geometry/util/select_calculation_type.hpp>
     17#include <boost/geometry/util/promote_floating_point.hpp>
     18#include <boost/geometry/util/math.hpp>
    1619
    1720#include <boost/geometry/extensions/gis/geographic/detail/ellipsoid.hpp>
     
    2831    \brief Point-point distance approximation taking flattening into account
    2932    \ingroup distance
    30     \tparam P1 first point type
    31     \tparam P2 optional second point type
     33    \tparam Point1 first point type
     34    \tparam Point2 optional second point type
    3235    \author After Andoyer, 19xx, republished 1950, republished by Meeus, 1999
    3336    \note Although not so well-known, the approximation is very good: in all cases the results
     
    4245template
    4346<
    44     typename P1,
    45     typename P2 = P1
    46     // calculation_type
     47    typename Point1,
     48    typename Point2 = Point1,
     49    typename CalculationType = void
    4750>
    4851class andoyer
    4952{
    5053    public :
    51         typedef double return_type;
    52 
    53         andoyer()
     54    typedef typename promote_floating_point
     55        <
     56            typename select_calculation_type
     57                <
     58                    Point1,
     59                    Point2,
     60                    CalculationType
     61                >::type
     62        >::type calculation_type;
     63
     64        inline andoyer()
    5465            : m_ellipsoid()
    5566        {}
    5667
    57         explicit inline andoyer(double f)
     68        explicit inline andoyer(calculation_type f)
    5869            : m_ellipsoid(f)
    5970        {}
     
    6475
    6576
    66         inline return_type apply(P1 const& p1, P2 const& p2) const
     77        inline calculation_type apply(Point1 const& point1, Point2 const& point2) const
    6778        {
    68             return calc(get_as_radian<0>(p1), get_as_radian<1>(p1),
    69                             get_as_radian<0>(p2), get_as_radian<1>(p2));
     79            return calc(get_as_radian<0>(point1), get_as_radian<1>(point1),
     80                            get_as_radian<0>(point2), get_as_radian<1>(point2));
    7081        }
    7182
     
    7788
    7889    private :
    79         typedef typename coordinate_type<P1>::type T1;
    80         typedef typename coordinate_type<P2>::type T2;
     90        typedef typename coordinate_type<Point1>::type T1;
     91        typedef typename coordinate_type<Point2>::type T2;
    8192        geometry::detail::ellipsoid m_ellipsoid;
    8293
    83         inline return_type calc(T1 const& lon1, T1 const& lat1, T2 const& lon2, T2 const& lat2) const
     94        inline calculation_type calc(calculation_type const& lon1,
     95                    calculation_type const& lat1,
     96                    calculation_type const& lon2,
     97                    calculation_type const& lat2) const
    8498        {
    85             typedef double calculation_type;
    8699            calculation_type G = (lat1 - lat2) / 2.0;
    87100            calculation_type lambda = (lon1 - lon2) / 2.0;
     
    110123            }
    111124
     125            calculation_type const c1 = 1;
     126            calculation_type const c2 = 2;
     127            calculation_type const c3 = 3;
     128
    112129            calculation_type omega = atan(sqrt(S / C));
    113130            calculation_type r = sqrt(S * C) / omega; // not sure if this is r or greek nu
    114131
    115             calculation_type D = 2.0 * omega * m_ellipsoid.a();
    116             calculation_type H1 = (3 * r - 1.0) / (2.0 * C);
    117             calculation_type H2 = (3 * r + 1.0) / (2.0 * S);
    118 
    119             return return_type(D
    120                 * (1.0 + m_ellipsoid.f() * H1 * sinF2 * cosG2
    121                             - m_ellipsoid.f() * H2 * cosF2 * sinG2));
     132            calculation_type D = c2 * omega * m_ellipsoid.a();
     133            calculation_type H1 = (c3 * r - c1) / (c2 * C);
     134            calculation_type H2 = (c3 * r + c1) / (c2 * S);
     135
     136            calculation_type f = m_ellipsoid.f();
     137
     138            return D * (c1 + f * H1 * sinF2 * cosG2 - f * H2 * cosF2 * sinG2);
    122139        }
    123140};
     
    138155struct return_type<strategy::distance::andoyer<Point1, Point2> >
    139156{
    140     typedef typename strategy::distance::andoyer<Point1, Point2>::return_type type;
     157    typedef typename strategy::distance::andoyer<Point1, Point2>::calculation_type type;
    141158};
    142159
  • sandbox/geometry/boost/geometry/extensions/gis/geographic/strategies/vincenty.hpp

    r63707 r63839  
    6262                            CalculationType
    6363                        >::type,
    64                     double // it should be at least double otherwise Vincenty does not run
     64                    double // to avoid bad results in float
    6565                >::type
    66         >::type return_type;
     66        >::type calculation_type;
    6767
    6868    inline vincenty()
     
    7373    {}
    7474
    75     inline return_type apply(Point1 const& p1, Point2 const& p2) const
     75    inline calculation_type apply(Point1 const& p1, Point2 const& p2) const
    7676    {
    7777        return calculate(get_as_radian<0>(p1), get_as_radian<1>(p1),
     
    8686
    8787private :
    88     typedef return_type promoted_type;
    8988    geometry::detail::ellipsoid m_ellipsoid;
    9089
    91     inline return_type calculate(promoted_type const& lon1,
    92                 promoted_type const& lat1,
    93                 promoted_type const& lon2,
    94                 promoted_type const& lat2) const
    95     {
    96         namespace mc = boost::math::constants;
    97 
    98         promoted_type const c2 = 2;
    99         promoted_type const two_pi = c2 * mc::pi<promoted_type>();
     90    inline calculation_type calculate(calculation_type const& lon1,
     91                calculation_type const& lat1,
     92                calculation_type const& lon2,
     93                calculation_type const& lat2) const
     94    {
     95        calculation_type const c2 = 2;
     96        calculation_type const pi = geometry::math::pi<calculation_type>();
     97        calculation_type const two_pi = c2 * pi;
    10098
    10199        // lambda: difference in longitude on an auxiliary sphere
    102         promoted_type L = lon2 - lon1;
    103         promoted_type lambda = L;
    104 
    105         if (L < -mc::pi<promoted_type>()) L += two_pi;
    106         if (L > mc::pi<promoted_type>()) L -= two_pi;
     100        calculation_type L = lon2 - lon1;
     101        calculation_type lambda = L;
     102
     103        if (L < -pi) L += two_pi;
     104        if (L > pi) L -= two_pi;
    107105
    108106        if (math::equals(lat1, lat2) && math::equals(lon1, lon2))
    109107        {
    110             return return_type(0);
     108            return calculation_type(0);
    111109        }
    112110
    113111        // TODO: give ellipsoid a template-parameter
    114         promoted_type const ellipsoid_f = m_ellipsoid.f();
    115         promoted_type const ellipsoid_a = m_ellipsoid.a();
    116         promoted_type const ellipsoid_b = m_ellipsoid.b();
     112        calculation_type const ellipsoid_f = m_ellipsoid.f();
     113        calculation_type const ellipsoid_a = m_ellipsoid.a();
     114        calculation_type const ellipsoid_b = m_ellipsoid.b();
    117115
    118116        // U: reduced latitude, defined by tan U = (1-f) tan phi
    119         promoted_type const c1 = 1;
    120         promoted_type const one_min_f = c1 - ellipsoid_f;
    121 
    122         promoted_type const U1 = atan(one_min_f * tan(lat1)); // above (1)
    123         promoted_type const U2 = atan(one_min_f * tan(lat2)); // above (1)
    124 
    125         promoted_type const cos_U1 = cos(U1);
    126         promoted_type const cos_U2 = cos(U2);
    127         promoted_type const sin_U1 = sin(U1);
    128         promoted_type const sin_U2 = sin(U2);
     117        calculation_type const c1 = 1;
     118        calculation_type const one_min_f = c1 - ellipsoid_f;
     119
     120        calculation_type const U1 = atan(one_min_f * tan(lat1)); // above (1)
     121        calculation_type const U2 = atan(one_min_f * tan(lat2)); // above (1)
     122
     123        calculation_type const cos_U1 = cos(U1);
     124        calculation_type const cos_U2 = cos(U2);
     125        calculation_type const sin_U1 = sin(U1);
     126        calculation_type const sin_U2 = sin(U2);
    129127
    130128        // alpha: azimuth of the geodesic at the equator
    131         promoted_type cos2_alpha;
    132         promoted_type sin_alpha;
     129        calculation_type cos2_alpha;
     130        calculation_type sin_alpha;
    133131
    134132        // sigma: angular distance p1,p2 on the sphere
    135133        // sigma1: angular distance on the sphere from the equator to p1
    136134        // sigma_m: angular distance on the sphere from the equator to the midpoint of the line
    137         promoted_type sigma;
    138         promoted_type sin_sigma;
    139         promoted_type cos2_sigma_m;
    140 
    141         promoted_type previous_lambda;
    142 
    143         promoted_type const c3 = 3;
    144         promoted_type const c4 = 4;
    145         promoted_type const c6 = 6;
    146         promoted_type const c16 = 16;
    147 
    148         promoted_type const c_e_12 = 1e-12;
     135        calculation_type sigma;
     136        calculation_type sin_sigma;
     137        calculation_type cos2_sigma_m;
     138
     139        calculation_type previous_lambda;
     140
     141        calculation_type const c3 = 3;
     142        calculation_type const c4 = 4;
     143        calculation_type const c6 = 6;
     144        calculation_type const c16 = 16;
     145
     146        calculation_type const c_e_12 = 1e-12;
    149147
    150148        do
    151149        {
    152150            previous_lambda = lambda; // (13)
    153             promoted_type sin_lambda = sin(lambda);
    154             promoted_type cos_lambda = cos(lambda);
     151            calculation_type sin_lambda = sin(lambda);
     152            calculation_type cos_lambda = cos(lambda);
    155153            sin_sigma = sqrt(math::sqr(cos_U2 * sin_lambda) + math::sqr(cos_U1 * sin_U2 - sin_U1 * cos_U2 * cos_lambda)); // (14)
    156             promoted_type cos_sigma = sin_U1 * sin_U2 + cos_U1 * cos_U2 * cos_lambda; // (15)
     154            calculation_type cos_sigma = sin_U1 * sin_U2 + cos_U1 * cos_U2 * cos_lambda; // (15)
    157155            sin_alpha = cos_U1 * cos_U2 * sin_lambda / sin_sigma; // (17)
    158156            cos2_alpha = c1 - math::sqr(sin_alpha);
    159             cos2_sigma_m = cos2_alpha == 0 ? 0 : cos_sigma - c2 * sin_U1 * sin_U2 / cos2_alpha; // (18)
    160 
    161 
    162             promoted_type C = ellipsoid_f/c16 * cos2_alpha * (c4 + ellipsoid_f * (c4 - c3 * cos2_alpha)); // (10)
     157            cos2_sigma_m = math::equals(cos2_alpha, 0) ? 0 : cos_sigma - c2 * sin_U1 * sin_U2 / cos2_alpha; // (18)
     158
     159            calculation_type C = ellipsoid_f/c16 * cos2_alpha * (c4 + ellipsoid_f * (c4 - c3 * cos2_alpha)); // (10)
    163160            sigma = atan2(sin_sigma, cos_sigma); // (16)
    164161            lambda = L + (c1 - C) * ellipsoid_f * sin_alpha *
     
    166163
    167164        } while (geometry::math::abs(previous_lambda - lambda) > c_e_12
    168                 && geometry::math::abs(lambda) < mc::pi<promoted_type>());
    169 
    170         promoted_type sqr_u = cos2_alpha * (math::sqr(ellipsoid_a) - math::sqr(ellipsoid_b)) / math::sqr(ellipsoid_b); // above (1)
     165                && geometry::math::abs(lambda) < pi);
     166
     167        calculation_type sqr_u = cos2_alpha * (math::sqr(ellipsoid_a) - math::sqr(ellipsoid_b)) / math::sqr(ellipsoid_b); // above (1)
    171168
    172169        // Oops getting hard here
    173170        // (again, problem is that ttmath cannot divide by doubles, which is OK)
    174         promoted_type const c47 = 47;
    175         promoted_type const c74 = 74;
    176         promoted_type const c128 = 128;
    177         promoted_type const c256 = 256;
    178         promoted_type const c175 = 175;
    179         promoted_type const c320 = 320;
    180         promoted_type const c768 = 768;
    181         promoted_type const c1024 = 1024;
    182         promoted_type const c4096 = 4096;
    183         promoted_type const c16384 = 16384;
    184 
    185         promoted_type A = c1 + sqr_u/c16384 * (c4096 + sqr_u * (-c768 + sqr_u * (c320 - c175 * sqr_u))); // (3)
    186         promoted_type B = sqr_u/c1024 * (c256 + sqr_u * ( -c128 + sqr_u * (c74 - c47 * sqr_u))); // (4)
    187         promoted_type delta_sigma = B * sin_sigma * ( cos2_sigma_m + (B/c4) * (cos(sigma)* (-c1 + c2 * cos2_sigma_m)
     171        calculation_type const c47 = 47;
     172        calculation_type const c74 = 74;
     173        calculation_type const c128 = 128;
     174        calculation_type const c256 = 256;
     175        calculation_type const c175 = 175;
     176        calculation_type const c320 = 320;
     177        calculation_type const c768 = 768;
     178        calculation_type const c1024 = 1024;
     179        calculation_type const c4096 = 4096;
     180        calculation_type const c16384 = 16384;
     181
     182        calculation_type A = c1 + sqr_u/c16384 * (c4096 + sqr_u * (-c768 + sqr_u * (c320 - c175 * sqr_u))); // (3)
     183        calculation_type B = sqr_u/c1024 * (c256 + sqr_u * ( -c128 + sqr_u * (c74 - c47 * sqr_u))); // (4)
     184        calculation_type delta_sigma = B * sin_sigma * ( cos2_sigma_m + (B/c4) * (cos(sigma)* (-c1 + c2 * cos2_sigma_m)
    188185                - (B/c6) * cos2_sigma_m * (-c3 + c4 * math::sqr(sin_sigma)) * (-c3 + c4 * cos2_sigma_m))); // (6)
    189186
    190         promoted_type dist = ellipsoid_b * A * (sigma - delta_sigma); // (19)
    191 
    192         return return_type(dist);
     187        return ellipsoid_b * A * (sigma - delta_sigma); // (19)
    193188    }
    194189};
     
    208203struct return_type<strategy::distance::vincenty<Point1, Point2> >
    209204{
    210     typedef typename strategy::distance::vincenty<Point1, Point2>::return_type type;
     205    typedef typename strategy::distance::vincenty<Point1, Point2>::calculation_type type;
    211206};
    212207
     
    248243{
    249244    template <typename T>
    250     static inline typename vincenty<Point1, Point2>::return_type apply(vincenty<Point1, Point2> const& , T const& value)
     245    static inline typename return_type<vincenty<Point1, Point2> >::type apply(vincenty<Point1, Point2> const& , T const& value)
    251246    {
    252247        return value;
  • sandbox/geometry/boost/geometry/util/math.hpp

    r63593 r63839  
    2020{
    2121
     22namespace math
     23{
    2224
    2325#ifndef DOXYGEN_NO_DETAIL
    24 namespace detail { namespace math
     26namespace detail
    2527{
    2628
     
    4749
    4850
    49 }} // namespace detail::math
     51/*!
     52\brief Short construct to enable partial specialization for PI, currently not possible in Math.
     53*/
     54template <typename T>
     55struct define_pi
     56{
     57    static inline T apply()
     58    {
     59        // Default calls Boost.Math
     60        return boost::math::constants::pi<T>();
     61    }
     62};
     63
     64
     65} // namespace detail
    5066#endif
    5167
    52 namespace math
    53 {
     68
     69template <typename T>
     70inline T pi() { return detail::define_pi<T>::apply(); }
    5471
    5572
     
    7390{
    7491    typedef typename select_most_precise<T1, T2>::type select_type;
    75     return detail::math::equals
     92    return detail::equals
    7693        <
    7794            select_type,
     
    8198
    8299
    83 double const d2r = boost::math::constants::pi<double>() / 180.0;
     100double const d2r = geometry::math::pi<double>() / 180.0;
    84101double const r2d = 1.0 / d2r;
    85102
     
    122139}
    123140
     141
    124142} // namespace math
    125143
Note: See TracChangeset for help on using the changeset viewer.