GeographicLib 2.5
Rhumb.cpp
Go to the documentation of this file.
1/**
2 * \file Rhumb.cpp
3 * \brief Implementation for GeographicLib::Rhumb and GeographicLib::RhumbLine
4 * classes
5 *
6 * Copyright (c) Charles Karney (2014-2023) <karney@alum.mit.edu> and licensed
7 * under the MIT/X11 License. For more information, see
8 * https://geographiclib.sourceforge.io/
9 **********************************************************************/
10
12#include <GeographicLib/DST.hpp>
13
14namespace GeographicLib {
15
16 using namespace std;
17
18 Rhumb::Rhumb(real a, real f, bool exact)
19 : _aux(a, f)
20 , _exact(exact)
21 , _a(a)
22 , _f(f)
23 , _n(_f / (2 - _f))
24 , _rm(_aux.RectifyingRadius(_exact))
25 , _c2(_aux.AuthalicRadiusSquared(_exact) * Math::degree())
26 , _lL(_exact ? 8 : Lmax_) // 8 is starting size for DFT fit
27 , _pP(_lL)
28 {
29 AreaCoeffs();
30 }
31
33 static const Rhumb
35 return wgs84;
36 }
37
38 void Rhumb::AreaCoeffs() {
39 // Set up coefficients for area calculation
40 if (_exact) {
41 // Compute coefficients by Fourier transform of integrand
42 static const real eps = numeric_limits<real>::epsilon()/2;
43 qIntegrand f(_aux);
44 int L = 4;
45 vector<real> c(L);
46 DST fft(L); fft.transform(f, c.data()); L *= 2;
47 // For |n| <= 0.99, actual max for doubles is 2163. This scales as
48 // Math::digits() and for long doubles (GEOGRAPHICLIB_PRECISION = 3,
49 // digits = 64), this becomes 2163 * 64 / 53 = 2612. Round this up to
50 // 2^12 = 4096 and scale this by Math::digits()/64 if digits() > 64.
51 //
52 // 64 = digits for long double, 6 = 12 - log2(64)
53 int Lmax = 1<<(int(ceil(log2(max(Math::digits(), 64)))) + 6);
54 for (_lL = 0; L <= Lmax && _lL == 0; L *=2) {
55 fft.reset(L/2); c.resize(L); fft.refine(f, c.data());
56 _pP.resize(L);
57 for (int l = 0, k = -1; l < L; ++l) {
58 // Compute Fourier coefficients of integral
59 _pP[l] = (c[l] + (l+1 < L ? c[l+1] : 0)) / (-4 * (l+1));
60 if (fabs(_pP[l]) <= eps) {
61 if (k < 0) k = l; // mark as first small value
62 } else
63 k = -1; // run interrupted
64 if (k >= 0 && l - k + 1 >= (l + 1 + 7) / 8) {
65 // run of small values of at least l/8?
66 _lL = l + 1; _pP.resize(_lL); break;
67 }
68 }
69 // loop exits if _lL > 0
70 }
71 if (_lL == 0) // Hasn't converged -- just use the values we have
72 _lL = int(_pP.size());
73 } else {
74 // Use series expansions in n for Fourier coeffients of the integral
75 // See "Series expansions for computing rhumb areas"
76 // https://doi.org/10.5281/zenodo.7685484
77#if GEOGRAPHICLIB_RHUMBAREA_ORDER == 4
78 static const real coeffs[] = {
79 // Coefficients in matrix Q
80 596/real(2025), -398/real(945), 22/real(45), -1/real(3),
81 1543/real(4725), -118/real(315), 1/real(5),
82 152/real(945), -17/real(315),
83 5/real(252),
84 };
85#elif GEOGRAPHICLIB_RHUMBAREA_ORDER == 5
86 static const real coeffs[] = {
87 // Coefficients in matrix Q
88 -102614/real(467775), 596/real(2025), -398/real(945), 22/real(45),
89 -1/real(3),
90 -24562/real(155925), 1543/real(4725), -118/real(315), 1/real(5),
91 -38068/real(155925), 152/real(945), -17/real(315),
92 -752/real(10395), 5/real(252),
93 -101/real(17325),
94 };
95#elif GEOGRAPHICLIB_RHUMBAREA_ORDER == 6
96 static const real coeffs[] = {
97 // Coefficients in matrix Q
98 138734126/real(638512875), -102614/real(467775), 596/real(2025),
99 -398/real(945), 22/real(45), -1/real(3),
100 17749373/real(425675250), -24562/real(155925), 1543/real(4725),
101 -118/real(315), 1/real(5),
102 1882432/real(8513505), -38068/real(155925), 152/real(945),
103 -17/real(315),
104 268864/real(2027025), -752/real(10395), 5/real(252),
105 62464/real(2027025), -101/real(17325),
106 11537/real(4054050),
107 };
108#elif GEOGRAPHICLIB_RHUMBAREA_ORDER == 7
109 static const real coeffs[] = {
110 // Coefficients in matrix Q
111 -565017322/real(1915538625), 138734126/real(638512875),
112 -102614/real(467775), 596/real(2025), -398/real(945), 22/real(45),
113 -1/real(3),
114 -1969276/real(58046625), 17749373/real(425675250), -24562/real(155925),
115 1543/real(4725), -118/real(315), 1/real(5),
116 -58573784/real(638512875), 1882432/real(8513505), -38068/real(155925),
117 152/real(945), -17/real(315),
118 -6975184/real(42567525), 268864/real(2027025), -752/real(10395),
119 5/real(252),
120 -112832/real(1447875), 62464/real(2027025), -101/real(17325),
121 -4096/real(289575), 11537/real(4054050),
122 -311/real(525525),
123 };
124#elif GEOGRAPHICLIB_RHUMBAREA_ORDER == 8
125 static const real coeffs[] = {
126 // Coefficients in matrix Q
127 188270561816LL/real(488462349375LL), -565017322/real(1915538625),
128 138734126/real(638512875), -102614/real(467775), 596/real(2025),
129 -398/real(945), 22/real(45), -1/real(3),
130 2332829602LL/real(23260111875LL), -1969276/real(58046625),
131 17749373/real(425675250), -24562/real(155925), 1543/real(4725),
132 -118/real(315), 1/real(5),
133 -41570288/real(930404475), -58573784/real(638512875),
134 1882432/real(8513505), -38068/real(155925), 152/real(945),
135 -17/real(315),
136 1538774036/real(10854718875LL), -6975184/real(42567525),
137 268864/real(2027025), -752/real(10395), 5/real(252),
138 436821248/real(3618239625LL), -112832/real(1447875),
139 62464/real(2027025), -101/real(17325),
140 3059776/real(80405325), -4096/real(289575), 11537/real(4054050),
141 4193792/real(723647925), -311/real(525525),
142 1097653/real(1929727800),
143 };
144#else
145#error "Bad value for GEOGRAPHICLIB_RHUMBAREA_ORDER"
146#endif
147
148 static_assert(sizeof(coeffs) / sizeof(real) ==
149 (Lmax_ * (Lmax_ + 1))/2,
150 "Coefficient array size mismatch for Rhumb");
151 real d = 1;
152 int o = 0;
153 for (int l = 0; l < Lmax_; ++l) {
154 int m = Lmax_ - l - 1;
155 d *= _n;
156 _pP[l] = d * Math::polyval(m, coeffs + o, _n);
157 o += m + 1;
158 }
159 }
160 // Post condition: o == sizeof(coeffs) / sizeof(real)
161 }
162
163 Rhumb::qIntegrand::qIntegrand(const AuxLatitude& aux)
164 : _aux(aux) {}
165
166 Math::real Rhumb::qIntegrand::operator()(real beta) const {
167 // pbeta(beta) = integrate(q(beta), beta)
168 // q(beta) = (1-f) * (sin(xi) - sin(chi)) / cos(phi)
169 // = (1-f) * (cos(chi) - cos(xi)) / cos(phi) *
170 // (cos(xi) + cos(chi)) / (sin(xi) + sin(chi))
171 // Fit q(beta)/cos(beta) with Fourier transform
172 // q(beta)/cos(beta) = sum(c[k] * sin((2*k+1)*beta), k, 0, K-1)
173 // then the integral is
174 // pbeta = sum(d[k] * cos((2*k+2)*beta), k, 0, K-1)
175 // where
176 // d[k] = -1/(4*(k+1)) * (c[k] + c[k+1]) for k in 0..K-2
177 // d[K-1] = -1/(4*K) * c[K-1]
178 AuxAngle betaa(AuxAngle::radians(beta)),
179 phia(_aux.Convert(AuxLatitude::BETA, AuxLatitude::PHI,
180 betaa, true).normalized()),
181 chia(_aux.Convert(AuxLatitude::PHI , AuxLatitude::CHI,
182 phia , true).normalized()),
183 xia (_aux.Convert(AuxLatitude::PHI , AuxLatitude::XI ,
184 phia , true).normalized());
185 real schi = chia.y(), cchi = chia.x(), sxi = xia.y(), cxi = xia.x(),
186 cphi = phia.x(), cbeta = betaa.x();
187 return (1 - _aux.Flattening()) *
188 ( fabs(schi) < fabs(cchi) ? sxi - schi :
189 (cchi - cxi) * (cxi + cchi) / (sxi + schi) ) / (cphi * cbeta);
190 // Value for beta = pi/2. This isn't needed because beta is given in
191 // radians and cos(pi/2) is never exactly 0. See formulas in the auxlat
192 // paper for tan(chi)/tan(phi) and tan(xi)/tan(phi) in the limit phi ->
193 // pi/2.
194 //
195 // n = 1/2;
196 // e = 2*sqrt(n) / (1 + n);
197 // e1 = 2*sqrt(n) / (1 - n);
198 // at = f == 0 ? 1 : (f < 0 ? atan(|e|)/|e| : asinh(e1)/e);
199 // q1 = at + 1 / (1 - e^2);
200 // s1 = sinh(e^2 * at);
201 // h1 = f < 0 ? hypot(1, s1) - s1 : 1/(hypot(1, s1) + s1);
202 // v = (1 - e^2) / (2 * h1^2) - 1 / ((1 - e^2) * q1);
203 }
204
205 void Rhumb::GenInverse(real lat1, real lon1, real lat2, real lon2,
206 unsigned outmask,
207 real& s12, real& azi12, real& S12) const {
208 using std::isinf; // Needed for Centos 7, ubuntu 14
209 AuxAngle phi1(AuxAngle::degrees(lat1)), phi2(AuxAngle::degrees(lat2)),
210 chi1(_aux.Convert(_aux.PHI, _aux.CHI, phi1, _exact)),
211 chi2(_aux.Convert(_aux.PHI, _aux.CHI, phi2, _exact));
212 real
213 lon12 = Math::AngDiff(lon1, lon2),
214 lam12 = lon12 * Math::degree<real>(),
215 psi1 = chi1.lam(),
216 psi2 = chi2.lam(),
217 psi12 = psi2 - psi1;
218 if (outmask & AZIMUTH)
219 azi12 = Math::atan2d(lam12, psi12);
220 if (outmask & DISTANCE) {
221 if (isinf(psi1) || isinf(psi2)) {
223 phi2, _exact).radians() -
225 phi1, _exact).radians()) * _rm;
226 } else {
227 real h = hypot(lam12, psi12);
228 // dmu/dpsi = dmu/dchi / dpsi/dchi
229 real dmudpsi = _exact ?
230 _aux.DRectifying(phi1, phi2) / _aux.DIsometric(phi1, phi2) :
231 _aux.DConvert(AuxLatitude::CHI, AuxLatitude::MU, chi1, chi2)
232 / DAuxLatitude::Dlam(chi1.tan(), chi2.tan());
233 s12 = h * dmudpsi * _rm;
234 }
235 }
236 if (outmask & AREA)
237 S12 = _c2 * lon12 * MeanSinXi(chi1, chi2);
238 }
239
240 RhumbLine Rhumb::Line(real lat1, real lon1, real azi12) const
241 { return RhumbLine(*this, lat1, lon1, azi12); }
242
243 void Rhumb::GenDirect(real lat1, real lon1, real azi12, real s12,
244 unsigned outmask,
245 real& lat2, real& lon2, real& S12) const
246 { Line(lat1, lon1, azi12).GenPosition(s12, outmask, lat2, lon2, S12); }
247
248 Math::real Rhumb::MeanSinXi(const AuxAngle& chix, const AuxAngle& chiy)
249 const {
251 phix (_aux.Convert(_aux.CHI, _aux.PHI , chix, _exact)),
252 phiy (_aux.Convert(_aux.CHI, _aux.PHI , chiy, _exact)),
253 betax(_aux.Convert(_aux.PHI, _aux.BETA, phix, _exact).normalized()),
254 betay(_aux.Convert(_aux.PHI, _aux.BETA, phiy, _exact).normalized());
255 real DpbetaDbeta =
257 betay.radians() - betax.radians(),
258 betax.y(), betax.x(), betay.y(), betay.x(),
259 _pP.data(), _lL),
260 tx = chix.tan(), ty = chiy.tan(),
261 DbetaDpsi = _exact ?
262 _aux.DParametric(phix, phiy) / _aux.DIsometric(phix, phiy) :
263 _aux.DConvert(AuxLatitude::CHI, AuxLatitude::BETA, chix, chiy) /
264 DAuxLatitude::Dlam(tx, ty);
265 return DAuxLatitude::Dp0Dpsi(tx, ty) + DpbetaDbeta * DbetaDpsi;
266 }
267
268 RhumbLine::RhumbLine(const Rhumb& rh, real lat1, real lon1, real azi12)
269 : _rh(rh)
270 , _lat1(Math::LatFix(lat1))
271 , _lon1(lon1)
272 , _azi12(Math::AngNormalize(azi12))
273 {
274 Math::sincosd(_azi12, _salp, _calp);
275 _phi1 = AuxAngle::degrees(lat1);
276 _mu1 = _rh._aux.Convert(AuxLatitude::PHI, AuxLatitude::MU,
277 _phi1, _rh._exact).degrees();
278 _chi1 = _rh._aux.Convert(AuxLatitude::PHI, AuxLatitude::CHI,
279 _phi1, _rh._exact);
280 _psi1 = _chi1.lam();
281 }
282
283 void RhumbLine::GenPosition(real s12, unsigned outmask,
284 real& lat2, real& lon2, real& S12) const {
285 real
286 r12 = s12 / (_rh._rm * Math::degree()), // scaled distance in degrees
287 mu12 = r12 * _calp,
288 mu2 = _mu1 + mu12;
289 real lat2x, lon2x;
290 if (fabs(mu2) <= Math::qd) {
291 AuxAngle mu2a(AuxAngle::degrees(mu2)),
293 mu2a, _rh._exact)),
295 phi2, _rh._exact));
296 lat2x = phi2.degrees();
297 real dmudpsi = _rh._exact ?
298 _rh._aux.DRectifying(_phi1, phi2) / _rh._aux.DIsometric(_phi1, phi2) :
299 _rh._aux.DConvert(AuxLatitude::CHI, AuxLatitude::MU, _chi1, chi2)
300 / DAuxLatitude::Dlam(_chi1.tan(), chi2.tan());
301 lon2x = r12 * _salp / dmudpsi;
302 if (outmask & AREA)
303 S12 = _rh._c2 * lon2x * _rh.MeanSinXi(_chi1, chi2);
304 lon2x = outmask & LONG_UNROLL ? _lon1 + lon2x :
306 } else {
307 // Reduce to the interval [-180, 180)
308 mu2 = Math::AngNormalize(mu2);
309 // Deal with points on the anti-meridian
310 if (fabs(mu2) > Math::qd) mu2 = Math::AngNormalize(Math::hd - mu2);
311 lat2x = _rh._aux.Convert(AuxLatitude::MU, AuxLatitude::PHI,
312 AuxAngle::degrees(mu2), _rh._exact).degrees();
313 lon2x = Math::NaN();
314 if (outmask & AREA)
315 S12 = Math::NaN();
316 }
317 if (outmask & LATITUDE) lat2 = lat2x;
318 if (outmask & LONGITUDE) lon2 = lon2x;
319 }
320
321} // namespace GeographicLib
Header for GeographicLib::DST class.
GeographicLib::Math::real real
Definition GeodSolve.cpp:28
Header for GeographicLib::Rhumb and GeographicLib::RhumbLine classes.
An accurate representation of angles.
Definition AuxAngle.hpp:47
Math::real radians() const
Definition AuxAngle.hpp:228
AuxAngle normalized() const
Definition AuxAngle.cpp:28
Math::real degrees() const
Definition AuxAngle.hpp:224
Math::real lam() const
Definition AuxAngle.hpp:232
Math::real tan() const
Definition AuxAngle.hpp:113
AuxAngle Convert(int auxin, int auxout, const AuxAngle &zeta, bool exact=false) const
static Math::real Dlam(real x, real y)
Math::real DParametric(const AuxAngle &phi1, const AuxAngle &phi2) const
Math::real DConvert(int auxin, int auxout, const AuxAngle &zeta1, const AuxAngle &zeta2) const
Math::real DIsometric(const AuxAngle &phi1, const AuxAngle &phi2) const
static Math::real Dp0Dpsi(real x, real y)
static Math::real DClenshaw(bool sinp, real Delta, real szeta1, real czeta1, real szeta2, real czeta2, const real c[], int K)
Math::real DRectifying(const AuxAngle &phi1, const AuxAngle &phi2) const
Discrete sine transforms.
Definition DST.hpp:63
Mathematical functions needed by GeographicLib.
Definition Math.hpp:77
static T degree()
Definition Math.hpp:209
static void sincosd(T x, T &sinx, T &cosx)
Definition Math.cpp:101
static T atan2d(T y, T x)
Definition Math.cpp:199
static constexpr int qd
degrees per quarter turn
Definition Math.hpp:145
static T AngNormalize(T x)
Definition Math.cpp:66
static int digits()
Definition Math.cpp:21
static T NaN()
Definition Math.cpp:277
static T polyval(int N, const T p[], T x)
Definition Math.hpp:280
static T AngDiff(T x, T y, T &e)
Definition Math.cpp:77
static constexpr int hd
degrees per half turn
Definition Math.hpp:148
Find a sequence of points on a single rhumb line.
Definition Rhumb.hpp:382
void GenPosition(real s12, unsigned outmask, real &lat2, real &lon2, real &S12) const
Definition Rhumb.cpp:283
Solve of the direct and inverse rhumb problems.
Definition Rhumb.hpp:80
RhumbLine Line(real lat1, real lon1, real azi12) const
Definition Rhumb.cpp:240
friend class RhumbLine
Definition Rhumb.hpp:83
Rhumb(real a, real f, bool exact=false)
Definition Rhumb.cpp:18
static const Rhumb & WGS84()
Definition Rhumb.cpp:32
Namespace for GeographicLib.