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
|
// Boost.Geometry
// Copyright (c) 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_GNOMONIC_SPHEROID_HPP
#define BOOST_GEOMETRY_FORMULAS_GNOMONIC_SPHEROID_HPP
#include <boost/geometry/core/radius.hpp>
#include <boost/geometry/util/condition.hpp>
#include <boost/geometry/util/math.hpp>
#include <boost/geometry/formulas/andoyer_inverse.hpp>
#include <boost/geometry/formulas/flattening.hpp>
#include <boost/geometry/formulas/thomas_inverse.hpp>
#include <boost/geometry/formulas/vincenty_direct.hpp>
#include <boost/geometry/formulas/vincenty_inverse.hpp>
namespace boost { namespace geometry { namespace formula
{
/*!
\brief Gnomonic projection on spheroid (ellipsoid of revolution).
\author See
- Charles F.F Karney, Algorithms for geodesics, 2011
https://arxiv.org/pdf/1109.4448.pdf
*/
template <
typename CT,
template <typename, bool, bool, bool, bool ,bool> class Inverse,
template <typename, bool, bool, bool, bool> class Direct
>
class gnomonic_spheroid
{
typedef Inverse<CT, false, true, true, true, true> inverse_type;
typedef typename inverse_type::result_type inverse_result;
typedef Direct<CT, false, false, true, true> direct_quantities_type;
typedef Direct<CT, true, false, false, false> direct_coordinates_type;
typedef typename direct_coordinates_type::result_type direct_result;
public:
template <typename Spheroid>
static inline bool forward(CT const& lon0, CT const& lat0,
CT const& lon, CT const& lat,
CT & x, CT & y,
Spheroid const& spheroid)
{
inverse_result i_res = inverse_type::apply(lon0, lat0, lon, lat, spheroid);
CT const& m = i_res.reduced_length;
CT const& M = i_res.geodesic_scale;
if (math::smaller_or_equals(M, CT(0)))
{
return false;
}
CT rho = m / M;
x = sin(i_res.azimuth) * rho;
y = cos(i_res.azimuth) * rho;
return true;
}
template <typename Spheroid>
static inline bool inverse(CT const& lon0, CT const& lat0,
CT const& x, CT const& y,
CT & lon, CT & lat,
Spheroid const& spheroid)
{
CT const a = get_radius<0>(spheroid);
CT const ds_threshold = a * std::numeric_limits<CT>::epsilon(); // TODO: 0 for non-fundamental type
CT const azimuth = atan2(x, y);
CT const rho = math::sqrt(math::sqr(x) + math::sqr(y)); // use hypot?
CT distance = a * atan(rho / a);
bool found = false;
for (int i = 0 ; i < 10 ; ++i)
{
direct_result d_res = direct_quantities_type::apply(lon0, lat0, distance, azimuth, spheroid);
CT const& m = d_res.reduced_length;
CT const& M = d_res.geodesic_scale;
if (math::smaller_or_equals(M, CT(0)))
{
// found = false;
return found;
}
CT const drho = m / M - rho; // rho = m / M
CT const ds = drho * math::sqr(M); // drho/ds = 1/M^2
distance -= ds;
// ds_threshold may be 0
if (math::abs(ds) <= ds_threshold)
{
found = true;
break;
}
}
if (found)
{
direct_result d_res = direct_coordinates_type::apply(lon0, lat0, distance, azimuth, spheroid);
lon = d_res.lon2;
lat = d_res.lat2;
}
return found;
}
};
}}} // namespace boost::geometry::formula
#endif // BOOST_GEOMETRY_FORMULAS_GNOMONIC_SPHEROID_HPP
|
