Coverage for pygeodesy/rhumbaux.py: 96%

163 statements  

« prev     ^ index     » next       coverage.py v7.2.2, created at 2023-08-12 12:31 -0400

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

2 

3u'''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 

6I{GeographicLib version 2.2+}. 

7 

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

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

10geodesic or perpendicular rhumb line. 

11 

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

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

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

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

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

17 

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

19License. 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 

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

26 

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

28from pygeodesy.auxilats.auxDLat import AuxDLat, _DClenshaw 

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

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

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

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

33 _0_0, _4_0, _720_0, _log2, _over 

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

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

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

37from pygeodesy.interns import NN, _COMMASPACE_ 

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

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

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

41 pairs, Property, Property_RO 

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

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

44 

45from math import ceil as _ceil, fabs, radians 

46 

47__all__ = _ALL_LAZY.rhumbaux 

48__version__ = '23.08.09' 

49 

50# DIGITS = (sizeof(real) * 8) bits # assert DIGITS > MANT_DIG 

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

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

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

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

55# 

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

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

58_rls = [] # instances of C{RbumbLineAux} to be updated 

59 

60 

61class RhumbAux(RhumbBase): 

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

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

64 

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

66 installed, version 1.16 or later. 

67 ''' 

68 

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

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

71 

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

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

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

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

76 a C{scalar}, ignored otherwise. 

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

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

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

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

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

82 C{TMorder}. 

83 

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

85 required when C{B{exact} is True}. 

86 

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

88 ''' 

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

90 if TMorder: 

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

92 

93 def areaux(self, **exact): 

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

95 

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

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

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

99 

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

101 

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

103 property C{f_max}. 

104 

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

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

107 

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

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

110 ''' 

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

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

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

114 return a 

115 

116 @Property_RO 

117 def _auxD(self): 

118 return AuxDLat(self.ellipsoid) 

119 

120 @Property_RO 

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

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

123 

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

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

126 

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

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

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

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

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

132 

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

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

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

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

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

138 

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

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

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

142 

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

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

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

146 carried out in finite terms. 

147 ''' 

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

149 name=self.name) 

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

151 

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

153 xD = self._auxD 

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

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

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

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

158 

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

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

161 

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

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

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

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

166 

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

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

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

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

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

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

173 

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

175 on opposite meridians, there are two shortest rhumb lines 

176 and the East-going one is chosen. 

177 

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

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

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

181 carried out in finite terms. 

182 ''' 

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

184 if (outmask & Cs.AZIMUTH_DISTANCE_AREA): 

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

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

187 

188 psi12 = psi2 - psi1 

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

190 lam12 = radians(lon12) 

191 if (outmask & Cs.AZIMUTH): 

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

193 if (outmask & Cs.DISTANCE): 

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

195 d = Phi2.toMu(self).toRadians 

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

197 s = fabs(d) 

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

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

200 s *= hypot(lam12, psi12) 

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

202 if (outmask & Cs.AREA): 

203 S = self._c2SinXi(Chi1, Chi2) 

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

205 return r 

206 

207 def _c2SinXi(self, Chix, Chiy): 

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

209 

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

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

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

213 Phiy.toBeta(self).normalized, 

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

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

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

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

218 dD += _Dp0Dpsi(tx, ty) 

219 dD *= self._c2 

220 return float(dD) 

221 

222 def _psiChiPhi3(self, lat): 

223 Phi = AuxPhi.fromDegrees(lat) 

224 Chi = Phi.toChi(self) 

225 psi = Chi.toLambertianRadians 

226 return psi, Chi, Phi 

227 

228 @Property_RO 

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

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

231 _RAseries(self._auxD)) 

232 

233 @Property_RO 

234 def _RhumbLine(self): 

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

236 ''' 

237 return RhumbLineAux 

238 

239 @Property_RO 

240 def _rrm(self): 

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

242 

243 @Property 

244 def TMorder(self): 

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

246 ''' 

247 return self._mTM 

248 

249 @TMorder.setter # PYCHOK setter! 

250 def TMorder(self, order): 

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

252 ''' 

253 self._TMorder(order) 

254 

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

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

257 

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

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

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

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

262 

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

264 ''' 

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

266 TMorder=self.TMorder) 

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

268 

269 

270class RhumbLineAux(RhumbLineBase): 

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

272 

273 Class C{RhumbLineAux} facilitates the determination of points 

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

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

276 

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

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

279 area C{S12} under the rhumb line. 

280 

281 Method C{RhumbLineAux.intersection2} finds the intersection 

282 between two rhumb lines. 

283 

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

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

286 ''' 

287 _Rhumb = RhumbAux # rhumbaux.RhumbAux 

288 

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

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

291 

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

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

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

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

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

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

298 values specifying the required capabilities. Include 

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

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

301 ''' 

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

303 

304 @Property_RO 

305 def _Chi1(self): 

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

307 

308 @Property_RO 

309 def _mu1(self): 

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

311 

312 @Property_RO 

313 def _Phi1(self): 

314 return AuxPhi.fromDegrees(self.lat1) 

315 

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

317 '''Compute a point at a distance on this rhumb line. 

318 

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

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

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

322 the quantities to be returned. 

323 

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

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

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

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

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

329 

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

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

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

333 

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

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

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

337 carried out in finite terms. 

338 ''' 

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

340 if (outmask & Cs.LATITUDE_LONGITUDE_AREA): 

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

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

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

344 Mu2 = AuxMu.fromDegrees(mu2) 

345 Phi2 = Mu2.toPhi(R) 

346 lat2 = Phi2.toDegrees 

347 if x90: # PYCHOK no cover 

348 lon2 = NAN 

349 if (outmask & Cs.AREA): 

350 r.set_(S12=NAN) 

351 else: 

352 Chi2 = Phi2.toChi(R) 

353 Chi1 = self._Chi1 

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

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

356 if (outmask & Cs.AREA): 

357 S = R._c2SinXi(Chi1, Chi2) 

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

359 if (outmask & Cs.LONGITUDE): 

360 if (outmask & Cs.LONG_UNROLL): 

361 lon2 += self.lon1 

362 else: 

363 lon2 = _norm180(self._lon12 + lon2) 

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

365 if (outmask & Cs.LATITUDE): 

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

367 if (outmask & Cs.LONGITUDE): 

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

369 return r 

370 

371# @Property_RO 

372# def _psi1(self): 

373# return self._Chi1.toLambertianRadians 

374 

375 

376def _RAintegrate(auxD): 

377 # Compute coefficients by Fourier transform of integrand 

378 L = 2 

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

380 f = auxD._qIntegrand 

381 c = fft.transform(f) 

382 # assert L < _Lbits 

383 while L < _Lbits: 

384 fft.reset(L) 

385 c = fft.refine(f, c) 

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

387 # assert len(c) == L 

388 pP, k = [], -1 

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

390 # Compute Fourier coefficients of integral 

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

392 if fabs(p) > EPS_2: 

393 k = -1 # run interrupted 

394 else: 

395 if k < 0: 

396 k = 1 # mark as first small value 

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

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

399 return pP[:j] 

400 pP.append(p) 

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

402 

403 

404def _RAseries(auxD): 

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

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

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

408 d = n = auxD._n 

409 i = 0 

410 pP = [] 

411 aL = auxD.ALorder 

412 Cs = _RACoeffs[aL] 

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

414 _p = _MODS.karney._polynomial 

415 for m in _reverange(aL): # order 

416 j = i + m + 1 

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

418 d *= n 

419 i = j 

420 # assert i == len(pP) 

421 return pP 

422 

423 

424_f, _u = float, _Ufloats() 

425_RACoeffs = { # Rhumb Area Coefficients in matrix Q 

426 4: _u( # GEOGRAPHICLIB_RHUMBAREA_ORDER == 4 

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

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

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

430 5 / _f(252)), 

431 5: _u( # GEOGRAPHICLIB_RHUMBAREA_ORDER == 5 

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

433 -1 / _f(3), 

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

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

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

437 -101 / _f(17325)), 

438 6: _u( # GEOGRAPHICLIB_RHUMBAREA_ORDER == 6 

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

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

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

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

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

444 -17 / _f(315), 

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

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

447 11537 / _f(4054050)), 

448 7: _u( # GEOGRAPHICLIB_RHUMBAREA_ORDER == 7 

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

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

451 -1 / _f(3), 

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

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

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

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

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

457 5 / _f(252), 

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

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

460 -311 / _f(525525)), 

461 8: _u( # GEOGRAPHICLIB_RHUMBAREA_ORDER == 8 

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

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

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

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

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

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

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

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

470 -17 / _f(315), 

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

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

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

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

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

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

477 1097653 / _f(1929727800)) 

478} 

479del _f, _u, _Ufloats 

480 

481 

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

483 

484if __name__ == '__main__': 

485 

486 from pygeodesy.lazily import printf 

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

488 

489 def _re(fmt, r3, x3): 

490 e3 = [] 

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

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

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

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

495 

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

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

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

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

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

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

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

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

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

505 

506# % python3 -m pygeodesy.rhumbaux 

507 

508# 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) 

509 

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

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

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