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
|
// Boost.Geometry (aka GGL, Generic Geometry Library)
// Copyright (c) 2007-2012 Barend Gehrels, Amsterdam, the Netherlands.
// This file was modified by Oracle on 2014.
// Modifications copyright (c) 2014, 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_ALGORITHMS_DETAIL_AZIMUTH_HPP
#define BOOST_GEOMETRY_ALGORITHMS_DETAIL_AZIMUTH_HPP
#include <boost/geometry/core/cs.hpp>
#include <boost/geometry/core/access.hpp>
#include <boost/geometry/core/radian_access.hpp>
#include <boost/geometry/core/tags.hpp>
#include <boost/geometry/util/math.hpp>
#include <boost/geometry/algorithms/not_implemented.hpp>
#include <boost/geometry/algorithms/detail/vincenty_inverse.hpp>
namespace boost { namespace geometry
{
// An azimuth is an angle between a vector/segment from origin to a point of
// interest and a reference vector. Typically north-based azimuth is used.
// North direction is used as a reference, angle is measured clockwise
// (North - 0deg, East - 90deg). For consistency in 2d cartesian CS
// the reference vector is Y axis, angle is measured clockwise.
// http://en.wikipedia.org/wiki/Azimuth
#ifndef DOXYGEN_NO_DISPATCH
namespace detail_dispatch
{
template <typename ReturnType, typename Tag>
struct azimuth
: not_implemented<Tag>
{};
template <typename ReturnType>
struct azimuth<ReturnType, geographic_tag>
{
template <typename P1, typename P2, typename Spheroid>
static inline ReturnType apply(P1 const& p1, P2 const& p2, Spheroid const& spheroid)
{
return geometry::detail::vincenty_inverse<ReturnType>
( get_as_radian<0>(p1), get_as_radian<1>(p1),
get_as_radian<0>(p2), get_as_radian<1>(p2),
spheroid ).azimuth12();
}
template <typename P1, typename P2>
static inline ReturnType apply(P1 const& p1, P2 const& p2)
{
return apply(p1, p2, srs::spheroid<ReturnType>());
}
};
template <typename ReturnType>
struct azimuth<ReturnType, spherical_equatorial_tag>
{
template <typename P1, typename P2, typename Sphere>
static inline ReturnType apply(P1 const& p1, P2 const& p2, Sphere const& /*unused*/)
{
// http://williams.best.vwh.net/avform.htm#Crs
ReturnType dlon = get_as_radian<0>(p2) - get_as_radian<0>(p1);
ReturnType cos_p2lat = cos(get_as_radian<1>(p2));
// An optimization which should kick in often for Boxes
//if ( math::equals(dlon, ReturnType(0)) )
//if ( get<0>(p1) == get<0>(p2) )
//{
// return - sin(get_as_radian<1>(p1)) * cos_p2lat);
//}
// "An alternative formula, not requiring the pre-computation of d"
// In the formula below dlon is used as "d"
return atan2(sin(dlon) * cos_p2lat,
cos(get_as_radian<1>(p1)) * sin(get_as_radian<1>(p2))
- sin(get_as_radian<1>(p1)) * cos_p2lat * cos(dlon));
}
template <typename P1, typename P2>
static inline ReturnType apply(P1 const& p1, P2 const& p2)
{
return apply(p1, p2, 0); // dummy model
}
};
template <typename ReturnType>
struct azimuth<ReturnType, spherical_polar_tag>
: azimuth<ReturnType, spherical_equatorial_tag>
{};
template <typename ReturnType>
struct azimuth<ReturnType, cartesian_tag>
{
template <typename P1, typename P2, typename Plane>
static inline ReturnType apply(P1 const& p1, P2 const& p2, Plane const& /*unused*/)
{
ReturnType x = get<0>(p2) - get<0>(p1);
ReturnType y = get<1>(p2) - get<1>(p1);
// NOTE: azimuth 0 is at Y axis, increasing right
// as in spherical/geographic where 0 is at North axis
return atan2(x, y);
}
template <typename P1, typename P2>
static inline ReturnType apply(P1 const& p1, P2 const& p2)
{
return apply(p1, p2, 0); // dummy model
}
};
} // detail_dispatch
#endif // DOXYGEN_NO_DISPATCH
#ifndef DOXYGEN_NO_DETAIL
namespace detail
{
/// Calculate azimuth between two points.
/// The result is in radians.
template <typename ReturnType, typename Point1, typename Point2>
inline ReturnType azimuth(Point1 const& p1, Point2 const& p2)
{
return detail_dispatch::azimuth
<
ReturnType,
typename geometry::cs_tag<Point1>::type
>::apply(p1, p2);
}
/// Calculate azimuth between two points.
/// The result is in radians.
template <typename ReturnType, typename Point1, typename Point2, typename Model>
inline ReturnType azimuth(Point1 const& p1, Point2 const& p2, Model const& model)
{
return detail_dispatch::azimuth
<
ReturnType,
typename geometry::cs_tag<Point1>::type
>::apply(p1, p2, model);
}
} // namespace detail
#endif // DOXYGEN_NO_DETAIL
}} // namespace boost::geometry
#endif // BOOST_GEOMETRY_ALGORITHMS_DETAIL_AZIMUTH_HPP
|
