24 , _rm(_aux.RectifyingRadius(_exact))
25 , _c2(_aux.AuthalicRadiusSquared(_exact) *
Math::degree())
26 , _lL(_exact ? 8 : Lmax_)
38 void Rhumb::AreaCoeffs() {
42 static const real eps = numeric_limits<real>::epsilon()/2;
46 DST fft(L); fft.transform(f, c.data()); L *= 2;
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());
57 for (
int l = 0, k = -1; l < L; ++l) {
59 _pP[l] = (c[l] + (l+1 < L ? c[l+1] : 0)) / (-4 * (l+1));
60 if (fabs(_pP[l]) <= eps) {
64 if (k >= 0 && l - k + 1 >= (l + 1 + 7) / 8) {
66 _lL = l + 1; _pP.resize(_lL);
break;
72 _lL = int(_pP.size());
77#if GEOGRAPHICLIB_RHUMBAREA_ORDER == 4
78 static const real coeffs[] = {
85#elif GEOGRAPHICLIB_RHUMBAREA_ORDER == 5
86 static const real coeffs[] = {
95#elif GEOGRAPHICLIB_RHUMBAREA_ORDER == 6
96 static const real coeffs[] = {
98 138734126/
real(638512875), -102614/
real(467775), 596/
real(2025),
100 17749373/
real(425675250), -24562/
real(155925), 1543/
real(4725),
102 1882432/
real(8513505), -38068/
real(155925), 152/
real(945),
105 62464/
real(2027025), -101/
real(17325),
108#elif GEOGRAPHICLIB_RHUMBAREA_ORDER == 7
109 static const real coeffs[] = {
111 -565017322/
real(1915538625), 138734126/
real(638512875),
114 -1969276/
real(58046625), 17749373/
real(425675250), -24562/
real(155925),
116 -58573784/
real(638512875), 1882432/
real(8513505), -38068/
real(155925),
118 -6975184/
real(42567525), 268864/
real(2027025), -752/
real(10395),
120 -112832/
real(1447875), 62464/
real(2027025), -101/
real(17325),
121 -4096/
real(289575), 11537/
real(4054050),
124#elif GEOGRAPHICLIB_RHUMBAREA_ORDER == 8
125 static const real coeffs[] = {
127 188270561816LL/
real(488462349375LL), -565017322/
real(1915538625),
128 138734126/
real(638512875), -102614/
real(467775), 596/
real(2025),
130 2332829602LL/
real(23260111875LL), -1969276/
real(58046625),
131 17749373/
real(425675250), -24562/
real(155925), 1543/
real(4725),
133 -41570288/
real(930404475), -58573784/
real(638512875),
134 1882432/
real(8513505), -38068/
real(155925), 152/
real(945),
136 1538774036/
real(10854718875LL), -6975184/
real(42567525),
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),
145#error "Bad value for GEOGRAPHICLIB_RHUMBAREA_ORDER"
148 static_assert(
sizeof(coeffs) /
sizeof(
real) ==
149 (Lmax_ * (Lmax_ + 1))/2,
150 "Coefficient array size mismatch for Rhumb");
153 for (
int l = 0; l < Lmax_; ++l) {
154 int m = Lmax_ - l - 1;
163 Rhumb::qIntegrand::qIntegrand(
const AuxLatitude& aux)
180 betaa,
true).normalized()),
182 phia ,
true).normalized()),
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);
205 void Rhumb::GenInverse(real lat1, real lon1, real lat2, real lon2,
207 real& s12, real& azi12, real& S12)
const {
214 lam12 = lon12 * Math::degree<real>(),
221 if (isinf(psi1) || isinf(psi2)) {
225 phi1, _exact).
radians()) * _rm;
227 real h = hypot(lam12, psi12);
229 real dmudpsi = _exact ?
233 s12 = h * dmudpsi * _rm;
237 S12 = _c2 * lon12 * MeanSinXi(chi1, chi2);
241 {
return RhumbLine(*
this, lat1, lon1, azi12); }
243 void Rhumb::GenDirect(real lat1, real lon1, real azi12, real s12,
245 real& lat2, real& lon2, real& S12)
const
246 {
Line(lat1, lon1, azi12).
GenPosition(s12, outmask, lat2, lon2, S12); }
257 betay.radians() - betax.radians(),
258 betax.y(), betax.x(), betay.y(), betay.x(),
260 tx = chix.
tan(), ty = chiy.
tan(),
268 RhumbLine::RhumbLine(
const Rhumb& rh,
real lat1,
real lon1,
real azi12)
270 , _lat1(Math::LatFix(lat1))
272 , _azi12(Math::AngNormalize(azi12))
277 _phi1, _rh._exact).degrees();
284 real& lat2, real& lon2, real& S12)
const {
296 lat2x = phi2.degrees();
297 real dmudpsi = _rh._exact ?
301 lon2x = r12 * _salp / dmudpsi;
303 S12 = _rh._c2 * lon2x * _rh.MeanSinXi(_chi1, chi2);
317 if (outmask &
LATITUDE) lat2 = lat2x;
Header for GeographicLib::DST class.
GeographicLib::Math::real real
Header for GeographicLib::Rhumb and GeographicLib::RhumbLine classes.
An accurate representation of angles.
Math::real radians() const
AuxAngle normalized() const
Math::real degrees() const
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.
Mathematical functions needed by GeographicLib.
static void sincosd(T x, T &sinx, T &cosx)
static T atan2d(T y, T x)
static constexpr int qd
degrees per quarter turn
static T AngNormalize(T x)
static T polyval(int N, const T p[], T x)
static T AngDiff(T x, T y, T &e)
static constexpr int hd
degrees per half turn
Find a sequence of points on a single rhumb line.
void GenPosition(real s12, unsigned outmask, real &lat2, real &lon2, real &S12) const
Solve of the direct and inverse rhumb problems.
RhumbLine Line(real lat1, real lon1, real azi12) const
Rhumb(real a, real f, bool exact=false)
static const Rhumb & WGS84()
Namespace for GeographicLib.