1# -*- coding: utf-8 -*-

2

3 u'''A pure Python version of I{Karney}'s C++ classes U{Rhumb

4<https://GeographicLib.SourceForge.io/C++/doc/classGeographicLib_1_1Rhumb.html>} and U{RhumbLine

5<https://GeographicLib.SourceForge.io/C++/doc/classGeographicLib_1_1RhumbLine.html>} from

6 I{GeographicLib version 2.2+}.

7

8 Class L{RhumbLine} has been enhanced with methods C{intersection2} and C{nearestOn4} to iteratively

9 find the intersection of two rhumb lines, respectively the nearest point on a rumb line along a

10 geodesic or perpendicular rhumb line.

11

12 For more details, see the I{2.2} U{GeographicLib<https://GeographicLib.SourceForge.io/C++/doc/index.html>}

13 documentation, especially the U{Class List<https://GeographicLib.SourceForge.io/C++/doc/annotated.html>},

14 the background information on U{Rhumb lines<https://GeographicLib.SourceForge.io/C++/doc/rhumb.html>},

15 utility U{RhumbSolve<https://GeographicLib.SourceForge.io/C++/doc/RhumbSolve.1.html>} and U{Online rhumb

16 line calculations<https://GeographicLib.SourceForge.io/cgi-bin/RhumbSolve>}.

17

18 Copyright (C) U{Charles Karney<mailto:Karney@Alum.MIT.edu>} (2022-2023) and licensed under the MIT/X11

19 License. For more information, see the U{GeographicLib<https://GeographicLib.SourceForge.io>} documentation.

20

21@note: Class L{AuxDST} requires package U{numpy<https://PyPI.org/project/numpy>} to be installed,

22 version 1.16 or newer and needed for I{exact} area calculations.

23'''

24# make sure int/int division yields float quotient

25 from __future__ import division as _; del _ # PYCHOK semicolon

26

27 from pygeodesy.auxilats.auxAngle import AuxMu, AuxPhi, atan2d, hypot

28 from pygeodesy.auxilats.auxDLat import AuxDLat, _DClenshaw

29# from pygeodesy.auxilats.auxDST import AuxDST # _MODS

30 from pygeodesy.auxilats.auxily import _Dlam, _Dp0Dpsi, _Ufloats

31 from pygeodesy.basics import _reverange, unsigned0, _zip, _xkwds_get

32 from pygeodesy.constants import EPS_2, MANT_DIG, NAN, PI4, isinf, \

33 _0_0, _4_0, _720_0, _log2, _over

34# from pygeodesy.ellipsoids import _EWGS84 # from .karney

35# from pygeodesy.errors import itemsorted, _xkwds_get # from .basics, ...

36 from pygeodesy.karney import Caps, _diff182, GDict, _norm180, _EWGS84

37# from pygeodesy.fmath import hypot # from .auxilats.auxAngle

38 from pygeodesy.interns import NN, _COMMASPACE_

39 from pygeodesy.lazily import _ALL_LAZY, _ALL_DOCS, _ALL_MODS as _MODS

40# from pygeodesy.props import Property, Property_RO # from .rhumbBase

41 from pygeodesy.rhumbBase import RhumbBase, RhumbLineBase, itemsorted, \

42 pairs, Property, Property_RO

43# from pygeodesy.streprs import pairs # from .rhumbBase

44# from pygeodesy.utily import atan2d # from .auxilats.auxAngle

45

46 from math import ceil as _ceil, fabs, radians

47

48__all__ = _ALL_LAZY.rhumbaux

49__version__ = '23.09.15'

50

51# DIGITS = (sizeof(real) * 8) bits

52# = (ctypes.sizeof(ctypes.c_double(1.0)) * 8) bits

53# For |n| <= 0.99, actual max for doubles is 2163. This scales

54# as DIGITS and for long doubles (GEOGRAPHICLIB_PRECISION = 3,

55# DIGITS = 64), this becomes 2163 * 64 / 53 = 2612. Round this

56# up to 2^12 = 4096 and scale this by DIGITS//64 if DIGITS > 64.

57#

58# 64 = DIGITS for long double, 6 = 12 - _log2(64)

59_Lbits = 1 << (int(_ceil(_log2(max(MANT_DIG, 64)))) + 6)

60

61

62 class RhumbAux(RhumbBase):

63 '''Class to solve the I{direct} and I{inverse rhumb} problems, based

64 on I{Auxiliary Latitudes} for accuracy near the poles.

65

66 @note: Package U{numpy<https://PyPI.org/project/numpy>} must be

67 installed, version 1.16 or later.

68 '''

69

70 def __init__(self, a_earth=_EWGS84, f=None, exact=True, name=NN, **TMorder): # PYCHOK signature

71 '''New C{rhumbaux.RhumbAux}.

72

73 @kwarg a_earth: This rhumb's earth model (L{Ellipsoid}, L{Ellipsoid2},

74 L{a_f2Tuple}, L{Datum}, 2-tuple C{(a, f)}) or the

75 (equatorial) radius (C{scalar}).

76 @kwarg f: The ellipsoid's flattening (C{scalar}), iff B{C{a_earth}} is

77 a C{scalar}, ignored otherwise.

78 @kwarg exact: If C{True}, use the exact expressions for the I{Auxiliary

79 Latitudes}, otherwise use the I{Fourier} series expansion

80 (C{bool}), see also property C{exact}.

81 @kwarg name: Optional name (C{str}).

82 @kwarg TMorder: Optional keyword argument B{C{TMorder}}, see property

83 C{TMorder}.

84

85 @raise ImportError: Package C{numpy} not found or not installed, only

86 required for area C{S12} when C{B{exact} is True}.

87

88 @raise RhumbError: Invalid B{C{a_earth}}, B{C{f}} or B{C{RA_TMorder}}.

89 '''

90 RhumbBase.__init__(self, a_earth, f, exact, name)

91 if TMorder:

92 self.Tmorder = _xkwds_get(TMorder, TMorder=RhumbBase._mTM)

93

94 def areaux(self, **exact):

95 '''Get this ellipsoid's B{C{exact}} surface area (C{meter} I{squared}).

96

97 @kwarg exact: Optional C{exact} (C{bool}), overriding this rhumb's

98 C{exact} setting, if C{True}, use the exact expression

99 for the authalic radius otherwise the I{Taylor} series.

100

101 @return: The (signed?) surface area (C{meter} I{squared}).

102

103 @raise AuxError: If C{B{exact}=False} and C{abs(flattening)} exceeds

104 property C{f_max}.

105

106 @note: The area of a polygon encircling a pole can be found by adding

107 C{areaux / 2} to the sum of C{S12} for each side of the polygon.

108

109 @see: U{The area of rhumb polygons<https://ArXiv.org/pdf/2303.03219.pdf>}

110 and method L{auxilats.AuxLat.AuthalicRadius2}.

111 '''

112 x = _xkwds_get(exact, exact=self.exact)

113 a = (self._c2 * _720_0) if bool(x) is self.exact else \

114 (self._auxD.AuthalicRadius2(exact=x, f_max=self.f_max) * PI4)

115 return a

116

117 @Property_RO

118 def _auxD(self):

119 return AuxDLat(self.ellipsoid)

120

121 @Property_RO

122 def _c2(self): # radians makes _c2 a factor per degree

123 return radians(self._auxD.AuthalicRadius2(exact=self.exact, f_max=self.f_max))

124

125 def Direct(self, lat1, lon1, azi12, s12, outmask=Caps.LATITUDE_LONGITUDE):

126 '''Solve the I{direct rhumb} problem, optionally with the area.

127

128 @arg lat1: Latitude of the first point (C{degrees90}).

129 @arg lon1: Longitude of the first point (C{degrees180}).

130 @arg azi12: Azimuth of the rhumb line (compass C{degrees}).

131 @arg s12: Distance along the rhumb line from the given to

132 the destination point (C{meter}), can be negative.

133

134 @return: L{GDict} with 2 up to 8 items C{lat2, lon2, a12, S12,

135 lat1, lon1, azi12, s12} with the destination point's

136 latitude C{lat2} and longitude C{lon2} in C{degrees},

137 the rhumb angle C{a12} in C{degrees} and area C{S12}

138 under the rhumb line in C{meter} I{squared}.

139

140 @raise ImportError: Package C{numpy} not found or not installed,

141 only required for area C{S12} when C{B{exact}

142 is True}.

143

144 @note: If B{C{s12}} is large enough that the rhumb line crosses

145 a pole, the longitude of the second point is indeterminate

146 and C{NAN} is returned for C{lon2} and area C{S12}.

147

148 @note: If the given point is a pole, the cosine of its latitude is

149 taken to be C{sqrt(L{EPS})}. This position is extremely

150 close to the actual pole and allows the calculation to be

151 carried out in finite terms.

152 '''

153 rl = RhumbLineAux(self, lat1, lon1, azi12, caps=Caps.LINE_OFF,

154 name=self.name)

155 return rl.Position(s12, outmask) # lat2, lon2, S12

156

157 def _DMu_DPsi(self, Phi1, Phi2, Chi1, Chi2):

158 xD = self._auxD

159 return _over(xD.DRectifying(Phi1, Phi2),

160 xD.DIsometric( Phi1, Phi2)) if self.exact else \

161 _over(xD.CRectifying(Chi1, Chi2),

162 _Dlam(Chi1.tan, Chi2.tan)) # not Lambertian!

163

164 def Inverse(self, lat1, lon1, lat2, lon2, outmask=Caps.AZIMUTH_DISTANCE):

165 '''Solve the I{inverse rhumb} problem.

166

167 @arg lat1: Latitude of the first point (C{degrees90}).

168 @arg lon1: Longitude of the first point (C{degrees180}).

169 @arg lat2: Latitude of the second point (C{degrees90}).

170 @arg lon2: Longitude of the second point (C{degrees180}).

171

172 @return: L{GDict} with 5 to 8 items C{azi12, s12, a12, S12,

173 lat1, lon1, lat2, lon2}, the rhumb line's azimuth C{azi12}

174 in compass C{degrees} between C{-180} and C{+180}, the

175 distance C{s12} and rhumb angle C{a12} between both points

176 in C{meter} respectively C{degrees} and the area C{S12}

177 under the rhumb line in C{meter} I{squared}.

178

179 @raise ImportError: Package C{numpy} not found or not installed,

180 only required for area C{S12} when C{B{exact}

181 is True}.

182

183 @note: The shortest rhumb line is found. If the end points are

184 on opposite meridians, there are two shortest rhumb lines

185 and the East-going one is chosen.

186

187 @note: If either point is a pole, the cosine of its latitude is

188 taken to be C{sqrt(L{EPS})}. This position is extremely

189 close to the actual pole and allows the calculation to be

190 carried out in finite terms.

191 '''

192 r, Cs = GDict(name=self.name), Caps

193 if (outmask & Cs.AZIMUTH_DISTANCE_AREA):

194 psi1, Chi1, Phi1 = self._psiChiPhi3(lat1)

195 psi2, Chi2, Phi2 = self._psiChiPhi3(lat2)

196

197 psi12 = psi2 - psi1

198 lon12, _ = _diff182(lon1, lon2, K_2_0=True)

199 lam12 = radians(lon12)

200 if (outmask & Cs.AZIMUTH):

201 r.set_(azi12=atan2d(lam12, psi12))

202 if (outmask & Cs.DISTANCE):

203 if isinf(psi1) or isinf(psi2): # PYCHOK no cover

204 d = Phi2.toMu(self).toRadians

205 d -= Phi1.toMu(self).toRadians

206 s = fabs(d)

207 else: # dmu/dpsi = dmu/dchi/dpsi/dchi

208 s = self._DMu_DPsi(Phi1, Phi2, Chi1, Chi2)

209 s *= hypot(lam12, psi12)

210 r.set_(s12=self._rrm * s)

211 if (outmask & Cs.AREA):

212 S = self._c2SinXi(Chi1, Chi2)

213 r.set_(S12=unsigned0(S * lon12)) # like .gx

214 return r

215

216 def _c2SinXi(self, Chix, Chiy):

217 pP, xD = self._RA, self._auxD

218

219 tx, Phix = Chix.tan, Chix.toPhi(self)

220 ty, Phiy = Chiy.tan, Chiy.toPhi(self)

221 dD = _DClenshaw(False, Phix.toBeta(self).normalized,

222 Phiy.toBeta(self).normalized,

223 pP, min(len(pP), 8)) # Fsum

224 dD *= _over(xD.DParametric(Phix, Phiy),

225 xD.DIsometric( Phix, Phiy)) if self.exact else \

226 _over(xD.CParametric(Chix, Chiy), _Dlam(tx, ty)) # not Lambertian!

227 dD += _Dp0Dpsi(tx, ty)

228 dD *= self._c2

229 return float(dD)

230

231 def _psiChiPhi3(self, lat):

232 Phi = AuxPhi.fromDegrees(lat)

233 Chi = Phi.toChi(self)

234 psi = Chi.toLambertianRadians

235 return psi, Chi, Phi

236

237 @Property_RO

238 def _RA(self): # get the coefficients for area calculation

239 return tuple(_RAintegrate(self._auxD) if self.exact else

240 _RAseries(self._auxD))

241

242 @Property_RO

243 def _RhumbLine(self):

244 '''(INTERNAL) Get this module's C{RhumbLineAux} class.

245 '''

246 return RhumbLineAux

247

248 @Property_RO

249 def _rrm(self):

250 return self._auxD.RectifyingRadius(exact=self.exact)

251

252 @Property

253 def TMorder(self):

254 '''Get the I{Transverse Mercator} order (C{int}, 4, 5, 6, 7 or 8).

255 '''

256 return self._mTM

257

258 @TMorder.setter # PYCHOK setter!

259 def TMorder(self, order):

260 '''Set the I{Transverse Mercator} order (C{int}, 4, 5, 6, 7 or 8).

261 '''

262 self._TMorder(order)

263

264 def toStr(self, prec=6, sep=_COMMASPACE_, **unused): # PYCHOK signature

265 '''Return this C{Rhumb} as string.

266

267 @kwarg prec: The C{float} precision, number of decimal digits (0..9).

268 Trailing zero decimals are stripped for B{C{prec}} values

269 of 1 and above, but kept for negative B{C{prec}} values.

270 @kwarg sep: Separator to join (C{str}).

271

272 @return: Tuple items (C{str}).

273 '''

274 d = dict(ellipsoid=self.ellipsoid, exact=self.exact,

275 TMorder=self.TMorder)

276 return sep.join(pairs(itemsorted(d, asorted=False), prec=prec))

277

278

279 class RhumbLineAux(RhumbLineBase):

280 '''Compute one or several points on a single rhumb line.

281

282 Class C{RhumbLineAux} facilitates the determination of points

283 on a single rhumb line. The starting point (C{lat1}, C{lon1})

284 and the azimuth C{azi12} are specified once.

285

286 Method C{RhumbLineAux.Position} returns the location of an

287 other point and optionally the distance C{s12} along and the

288 area C{S12} under the rhumb line.

289

290 Method C{RhumbLineAux.intersection2} finds the intersection

291 between two rhumb lines.

292

293 Method C{RhumbLineAux.nearestOn4} computes the nearest point

294 on and the distance to a rhumb line in different ways.

295 '''

296 _Rhumb = RhumbAux # rhumbaux.RhumbAux

297

298 def __init__(self, rhumb, lat1=0, lon1=0, azi12=None, **caps_name): # PYCHOK signature

299 '''New C{rhumbaux.RhumbLineAux}.

300

301 @arg rhumb: The rhumb reference (C{rhumbaux.RhumbAux}).

302 @kwarg lat1: Latitude of the start point (C{degrees90}).

303 @kwarg lon1: Longitude of the start point (C{degrees180}).

304 @kwarg azi12: Azimuth of this rhumb line (compass C{degrees}).

305 @kwarg caps_name: Optional keyword arguments C{B{name}=NN} and

306 C{B{caps}=0}, a bit-or'ed combination of L{Caps}

307 values specifying the required capabilities. Include

308 C{Caps.LINE_OFF} if updates to the B{C{rhumb}} should

309 I{not} be reflected in this rhumb line.

310 '''

311 RhumbLineBase.__init__(self, rhumb, lat1, lon1, azi12, **caps_name)

312

313 @Property_RO

314 def _Chi1(self):

315 return self._Phi1.toChi(self.rhumb)

316

317 @Property_RO

318 def _mu1(self):

319 return self._Phi1.toMu(self.rhumb).toDegrees

320

321 @Property_RO

322 def _Phi1(self):

323 return AuxPhi.fromDegrees(self.lat1)

324

325 def Position(self, s12, outmask=Caps.LATITUDE_LONGITUDE):

326 '''Compute a point at a given distance on this rhumb line.

327

328 @arg s12: The distance along this rhumb line between its origin

329 and the point (C{meters}), can be negative.

330 @kwarg outmask: Bit-or'ed combination of L{Caps} values specifying

331 the quantities to be returned.

332

333 @return: L{GDict} with 4 to 8 items C{azi12, a12, s12, S12, lat2,

334 lon2, lat1, lon1} with latitude C{lat2} and longitude

335 C{lon2} of the point in C{degrees}, the rhumb angle C{a12}

336 in C{degrees} from the start point of and the area C{S12}

337 under this rhumb line in C{meter} I{squared}.

338

339 @raise ImportError: Package C{numpy} not found or not installed,

340 only required for area C{S12} when C{B{exact}

341 is True}.

342

343 @note: If B{C{s12}} is large enough that the rhumb line crosses a

344 pole, the longitude of the second point is indeterminate and

345 C{NAN} is returned for C{lon2} and area C{S12}.

346

347 If the first point is a pole, the cosine of its latitude is

348 taken to be C{sqrt(L{EPS})}. This position is extremely

349 close to the actual pole and allows the calculation to be

350 carried out in finite terms.

351 '''

352 r, Cs = GDict(name=self.name), Caps

353 if (outmask & Cs.LATITUDE_LONGITUDE_AREA):

354 E, R = self.ellipsoid, self.rhumb

355 r12 = _over(s12, radians(R._rrm))

356 mu2, x90 = self._mu22(self._calp * r12, self._mu1)

357 Mu2 = AuxMu.fromDegrees(mu2)

358 Phi2 = Mu2.toPhi(R)

359 lat2 = Phi2.toDegrees

360 if x90: # PYCHOK no cover

361 lon2 = NAN

362 if (outmask & Cs.AREA):

363 r.set_(S12=NAN)

364 else:

365 Chi2 = Phi2.toChi(R)

366 Chi1 = self._Chi1

367 lon2 = R._DMu_DPsi(self._Phi1, Phi2, Chi1, Chi2)

368 lon2 = _over(self._salp * r12, lon2)

369 if (outmask & Cs.AREA):

370 S = R._c2SinXi(Chi1, Chi2)

371 r.set_(S12=unsigned0(S * lon2)) # like .gx

372 if (outmask & Cs.LONGITUDE):

373 if (outmask & Cs.LONG_UNROLL):

374 lon2 += self.lon1

375 else:

376 lon2 = _norm180(self._lon12 + lon2)

377 r.set_(azi12=self.azi12, s12=s12, a12=s12 / E._L_90)

378 if (outmask & Cs.LATITUDE):

379 r.set_(lat2=lat2, lat1=self.lat1)

380 if (outmask & Cs.LONGITUDE):

381 r.set_(lon2=lon2, lon1=self.lon1)

382 return r

383

384# @Property_RO

385# def _psi1(self):

386# return self._Chi1.toLambertianRadians

387

388

389 def _RAintegrate(auxD):

390 # Compute coefficients by Fourier transform of integrand

391 L = 2

392 fft = _MODS.auxilats.auxDST.AuxDST(L)

393 f = auxD._qIntegrand

394 c = fft.transform(f)

395 # assert L < _Lbits

396 while L < _Lbits:

397 fft.reset(L)

398 c = fft.refine(f, c)

399 L *= 2 # == len(c)

400 # assert len(c) == L

401 pP, k = [], -1

402 for j in range(1, L + 1):

403 # Compute Fourier coefficients of integral

404 p = -(c[j - 1] + (c[j] if j < L else _0_0)) / (_4_0 * j)

405 if fabs(p) > EPS_2:

406 k = -1 # run interrupted

407 else:

408 if k < 0:

409 k = 1 # mark as first small value

410 if (j - k) >= ((j + 7) // 8):

411 # run of at least (j - 1) // 8 small values

412 return pP[:j]

413 pP.append(p)

414 return pP # no convergence, use pP as-is

415

416

417 def _RAseries(auxD):

418 # Series expansions in n for Fourier coeffients of the integral

419 # @see: U{"Series expansions for computing rhumb areas"

420 # <https:#DOI.org/10.5281/zenodo.7685484>}.

421 d = n = auxD._n

422 i = 0

423 pP = []

424 aL = auxD.ALorder

425 Cs = _RACoeffs[aL]

426 # assert len(Cs) == (aL * (aL + 1)) // 2

427 _p = _MODS.karney._polynomial

428 for m in _reverange(aL): # order

429 j = i + m + 1

430 pP.append(_p(n, Cs, i, j) * d)

431 d *= n

432 i = j

433 # assert i == len(pP)

434 return pP

435

436

437_f, _u = float, _Ufloats()

438_RACoeffs = { # Rhumb Area Coefficients in matrix Q

439 4: _u( # GEOGRAPHICLIB_RHUMBAREA_ORDER == 4

440 596 / _f(2025), -398 / _f(945), 22 / _f(45), -1 / _f(3),

441 1543 / _f(4725), -118 / _f(315), 1 / _f(5),

442 152 / _f(945), -17 / _f(315),

443 5 / _f(252)),

444 5: _u( # GEOGRAPHICLIB_RHUMBAREA_ORDER == 5

445 -102614 / _f(467775), 596 / _f(2025), -398 / _f(945), 22 / _f(45),

446 -1 / _f(3),

447 -24562 / _f(155925), 1543 / _f(4725), -118 / _f(315), 1 / _f(5),

448 -38068 / _f(155925), 152 / _f(945), -17 / _f(315),

449 -752 / _f(10395), 5 / _f(252),

450 -101 / _f(17325)),

451 6: _u( # GEOGRAPHICLIB_RHUMBAREA_ORDER == 6

452 138734126 / _f(638512875), -102614 / _f(467775), 596 / _f(2025),

453 -398 / _f(945), 22 / _f(45), -1 / _f(3),

454 17749373 / _f(425675250), -24562 / _f(155925), 1543 / _f(4725),

455 -118 / _f(315), 1 / _f(5),

456 1882432 / _f(8513505), -38068 / _f(155925), 152 / _f(945),

457 -17 / _f(315),

458 268864 / _f(2027025), -752 / _f(10395), 5 / _f(252),

459 62464 / _f(2027025), -101 / _f(17325),

460 11537 / _f(4054050)),

461 7: _u( # GEOGRAPHICLIB_RHUMBAREA_ORDER == 7

462 -565017322 / _f(1915538625), 138734126 / _f(638512875),

463 -102614 / _f(467775), 596 / _f(2025), -398 / _f(945), 22 / _f(45),

464 -1 / _f(3),

465 -1969276 / _f(58046625), 17749373 / _f(425675250), -24562 / _f(155925),

466 1543 / _f(4725), -118 / _f(315), 1 / _f(5),

467 -58573784 / _f(638512875), 1882432 / _f(8513505), -38068 / _f(155925),

468 152 / _f(945), -17 / _f(315),

469 -6975184 / _f(42567525), 268864 / _f(2027025), -752 / _f(10395),

470 5 / _f(252),

471 -112832 / _f(1447875), 62464 / _f(2027025), -101 / _f(17325),

472 -4096 / _f(289575), 11537 / _f(4054050),

473 -311 / _f(525525)),

474 8: _u( # GEOGRAPHICLIB_RHUMBAREA_ORDER == 8

475 188270561816 / _f(488462349375), -565017322 / _f(1915538625),

476 138734126 / _f(638512875), -102614 / _f(467775), 596 / _f(2025),

477 -398 / _f(945), 22 / _f(45), -1 / _f(3),

478 2332829602 / _f(23260111875), -1969276 / _f(58046625),

479 17749373 / _f(425675250), -24562 / _f(155925), 1543 / _f(4725),

480 -118 / _f(315), 1 / _f(5),

481 -41570288 / _f(930404475), -58573784 / _f(638512875),

482 1882432 / _f(8513505), -38068 / _f(155925), 152 / _f(945),

483 -17 / _f(315),

484 1538774036 / _f(10854718875), -6975184 / _f(42567525),

485 268864 / _f(2027025), -752 / _f(10395), 5 / _f(252),

486 436821248 / _f(3618239625), -112832 / _f(1447875),

487 62464 / _f(2027025), -101 / _f(17325),

488 3059776 / _f(80405325), -4096 / _f(289575), 11537 / _f(4054050),

489 4193792 / _f(723647925), -311 / _f(525525),

490 1097653 / _f(1929727800))

491}

492 del _f, _u, _Ufloats

493

494

495__all__ += _ALL_DOCS(Caps, RhumbAux, RhumbLineAux)

496

497 if __name__ == '__main__':

498

499 from pygeodesy.lazily import printf

500 from pygeodesy import RhumbAux # PYCHOK RhumbAux is __main__.RhumbAux

501

502 def _re(fmt, r3, x3):

503 e3 = []

504 for r, x in _zip(r3, x3): # strict=True

505 e = fabs(r - x) / fabs(x)

506 e3.append('%.g' % (e,))

507 printf((fmt % r3) + ' rel errors: ' + ', '.join(e3))

508

509 # <https://GeographicLib.SourceForge.io/cgi-bin/RhumbSolveå -p 9> version 2.2

510 rhumb = RhumbAux(exact=True) # WGS84 default

511 printf('# %r\n', rhumb)

512 r = rhumb.Direct8(40.6, -73.8, 51, 5.5e6) # from JFK about NE

513 _re('# JFK NE lat2=%.8f, lon2=%.8f, S12=%.1f', (r.lat2, r.lon2, r.S12), (71.688899882813, 0.2555198244234, 44095641862956.11))

514 r = rhumb.Inverse8(40.6, -73.8, 51.6, -0.5) # JFK to LHR

515 _re('# JFK-LHR azi12=%.8f, s12=%.3f S12=%.1f', (r.azi12, r.s12, r.S12), (77.7683897102557, 5771083.38332803, 37395209100030.39))

516 r = rhumb.Inverse8(40.6, -73.8, 35.8, 140.3) # JFK to Tokyo Narita

517 _re('# JFK-NRT azi12=%.8f, s12=%.3f S12=%.1f', (r.azi12, r.s12, r.S12), (-92.38888798169965, 12782581.067684170, -63760642939072.50))

518

519# % python3 -m pygeodesy.rhumbaux

520

521# RhumbAux(TMorder=6, ellipsoid=Ellipsoid(name='WGS84', a=6378137, b=6356752.31424518, f_=298.25722356, f=0.00335281, f2=0.00336409, n=0.00167922, e=0.08181919, e2=0.00669438, e21=0.99330562, e22=0.0067395, e32=0.00335843, A=6367449.14582341, L=10001965.72931272, R1=6371008.77141506, R2=6371007.18091847, R3=6371000.79000916, Rbiaxial=6367453.63451633, Rtriaxial=6372797.5559594), exact=True)

522

523# JFK NE lat2=71.68889988, lon2=0.25551982, S12=44095641862956.1 rel errors: 4e-11, 2e-08, 5e-16

524# JFK-LHR azi12=77.76838971, s12=5771083.383 S12=37395209100030.3 rel errors: 3e-12, 5e-15, 6e-16

525# JFK-NRT azi12=-92.38888798, s12=12782581.068 S12=-63760642939072.5 rel errors: 2e-16, 3e-16, 6e-16