Coverage for pygeodesy/rhumbx.py: 98%

269 statements  

« prev     ^ index     » next       coverage.py v7.2.2, created at 2023-09-15 09:43 -0400

1 

2# -*- coding: utf-8 -*- 

3 

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

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

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

7I{GeographicLib version 2.0}. 

8 

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

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

11geodesic or perpendicular rhumb line. 

12 

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

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

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

16the utily U{RhumbSolve<https://GeographicLib.SourceForge.io/C++/doc/RhumbSolve.1.html>} and U{Online 

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

18 

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

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

21''' 

22# make sure int/int division yields float quotient 

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

24 

25from pygeodesy.basics import copysign0, neg, unsigned0, _zip 

26from pygeodesy.constants import NAN, PI_2, _0_0s, _0_0, _0_5, \ 

27 _1_0, _2_0, _4_0, _720_0, _over 

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

29from pygeodesy.errors import itemsorted, RhumbError, _Xorder 

30from pygeodesy.fmath import hypot, hypot1 

31# from pygeodesy.fsums import fsum1f_ # _MODS 

32from pygeodesy.interns import NN, _COMMASPACE_ 

33from pygeodesy.karney import _atan2d, Caps, _diff182, GDict, _GTuple, \ 

34 _norm180, _EWGS84 

35from pygeodesy.ktm import KTransverseMercator, _Xs, \ 

36 _AlpCoeffs, _BetCoeffs # PYCHOK used! 

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

38from pygeodesy.props import deprecated_method, Property, Property_RO, property_RO 

39from pygeodesy.rhumbBase import RhumbBase, RhumbLineBase, Int, pairs, \ 

40 sincos2_, _update_all_rls 

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

42# from pygeodesy.units import Int # from .rhumbBase 

43# from pygeodesy.utily import sincos2_ # from .rhumbBase 

44 

45from math import asinh, atan, cos, cosh, fabs, radians, sin, sinh, sqrt, tan 

46 

47__all__ = _ALL_LAZY.rhumbx 

48__version__ = '23.09.15' 

49 

50 

51class Rhumb(RhumbBase): 

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

53 I{elliptic functions} or I{Krüger} series expansion. 

54 

55 @see: The U{Detailed Description<https://GeographicLib.SourceForge.io/C++/doc/ 

56 classGeographicLib_1_1Rhumb.html>} of I{Karney}'s C++ C{Rhumb Class}. 

57 ''' 

58 _mRA = 6 # see .RAorder 

59 

60 def __init__(self, a_earth=_EWGS84, f=None, exact=True, name=NN, **RA_TMorder): 

61 '''New C{rhumbx.Rhumb}. 

62 

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

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

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

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

67 a C{scalar}, ignored otherwise. 

68 @kwarg exact: If C{True}, use an addition theorem for elliptic integrals 

69 to compute I{Divided differences}, otherwise use the I{Krüger} 

70 series expansion (C{bool} or C{None}), see also properties 

71 C{exact} and C{TMorder}. 

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

73 @kwarg RA_TMorder: Optional keyword arguments B{C{RAorder}} and B{C{TMorder}} 

74 to set the respective C{order}, see properties C{RAorder} 

75 and C{TMorder} and method C{orders}. 

76 

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

78 ''' 

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

80 if RA_TMorder: 

81 self.orders(**RA_TMorder) 

82 

83 @Property_RO 

84 def _A2(self): # Conformal2RectifyingCoeffs 

85 m = self.TMorder 

86 return _Xs(_AlpCoeffs, m, self.ellipsoid), m 

87 

88 @Property_RO 

89 def _B2(self): # Rectifying2ConformalCoeffs 

90 m = self.TMorder 

91 return _Xs(_BetCoeffs, m, self.ellipsoid), m 

92 

93 def _DConformal2Rectifying(self, x, y): # radians 

94 return _1_0 + (_sincosSeries(True, x, y, *self._A2) if self.f else _0_0) 

95 

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

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

98 

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

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

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

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

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

104 

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

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

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

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

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

110 

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

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

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

114 

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

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

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

118 carried out in finite terms. 

119 ''' 

120 rl = RhumbLine(self, lat1, lon1, azi12, caps=Caps.LINE_OFF, 

121 name=self.name) 

122 return rl.Position(s12, outmask | self._debug) # lat2, lon2, S12 

123 

124 @deprecated_method 

125 def Direct7(self, lat1, lon1, azi12, s12, outmask=Caps.LATITUDE_LONGITUDE_AREA): 

126 '''DEPRECATED, use method L{Rhumb.Direct8}. 

127 

128 @return: A I{DEPRECATED} L{Rhumb7Tuple}. 

129 ''' 

130 return self.Direct8(lat1, lon1, azi12, s12, outmask=outmask)._to7Tuple() 

131 

132 def _DIsometrict(self, phix, phiy, tphix, tphiy, _Dtan_phix_phiy): 

133 E = self.ellipsoid 

134 return _Dtan_phix_phiy * _Dasinh(tphix, tphiy) - \ 

135 _Dsin(phix, phiy) * _DeatanhE(sin(phix), sin(phiy), E) 

136 

137 def _DIsometric2Rectifyingd(self, psix, psiy): # degrees 

138 if self.exact: 

139 E = self.ellipsoid 

140 phix, phiy, tphix, tphiy = _Eaux4(E.auxIsometric, psix, psiy) 

141 t = _Dtant(phix - phiy, tphix, tphiy) 

142 r = _over(self._DRectifyingt( tphix, tphiy, t), 

143 self._DIsometrict(phix, phiy, tphix, tphiy, t)) 

144 else: 

145 x, y = radians(psix), radians(psiy) 

146 r = self._DConformal2Rectifying(_gd(x), _gd(y)) * _Dgd(x, y) 

147 return r 

148 

149 def _DRectifyingt(self, tphix, tphiy, _Dtan_phix_phiy): 

150 E = self.ellipsoid 

151 tbetx = E.f1 * tphix 

152 tbety = E.f1 * tphiy 

153 return (E.f1 * _Dtan_phix_phiy * E.b * PI_2 

154 * _DfEt( tbetx, tbety, self._eF) 

155 * _Datan(tbetx, tbety)) / E.L 

156 

157 def _DRectifying2Conformal(self, x, y): # radians 

158 return _1_0 - (_sincosSeries(True, x, y, *self._B2) if self.f else _0_0) 

159 

160 def _DRectifying2Isometricd(self, mux, muy): # degrees 

161 E = self.ellipsoid 

162 phix, phiy, tphix, tphiy = _Eaux4(E.auxRectifying, mux, muy) 

163 if self.exact: 

164 t = _Dtant(phix - phiy, tphix, tphiy) 

165 r = _over(self._DIsometrict(phix, phiy, tphix, tphiy, t), 

166 self._DRectifyingt( tphix, tphiy, t)) 

167 else: 

168 r = self._DRectifying2Conformal(radians(mux), radians(muy)) * \ 

169 _Dgdinv(E.es_taupf(tphix), E.es_taupf(tphiy)) 

170 return r 

171 

172 @Property_RO 

173 def _eF(self): 

174 '''(INTERNAL) Get the ellipsoid's elliptic function. 

175 ''' 

176 # .k2 = 0.006739496742276434 

177 return self._E._elliptic_e12 # _MODS.elliptic.Elliptic(-self._E._e12) 

178 

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

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

181 

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

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

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

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

186 

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

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

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

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

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

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

193 

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

195 on opposite meridians, there are two shortest rhumb lines 

196 and the East-going one is chosen. 

197 

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

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

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

201 carried out in finite terms. 

202 ''' 

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

204 if (outmask & Cs.AZIMUTH_DISTANCE_AREA): 

205 r.set_(lat1=lat1, lon1=lon1, lat2=lat2, lon2=lon2) 

206 E = self.ellipsoid 

207 psi1 = E.auxIsometric(lat1) 

208 psi2 = E.auxIsometric(lat2) 

209 psi12 = psi2 - psi1 

210 lon12, _ = _diff182(lon1, lon2) 

211 if (outmask & Cs.AZIMUTH): 

212 r.set_(azi12=_atan2d(lon12, psi12)) 

213 if (outmask & Cs.DISTANCE): 

214 a12 = hypot(lon12, psi12) * self._DIsometric2Rectifyingd(psi2, psi1) 

215 s12 = a12 * E._L_90 

216 r.set_(s12=s12, a12=copysign0(a12, s12)) 

217 if (outmask & Cs.AREA): 

218 r.set_(S12=self._S12d(lon12, psi2, psi1)) 

219 if ((outmask | self._debug) & Cs._DEBUG_INVERSE): # PYCHOK no cover 

220 r.set_(a=E.a, f=E.f, f1=E.f1, L=E.L, 

221 b=E.b, e=E.e, e2=E.e2, k2=self._eF.k2, 

222 lon12=lon12, psi1=psi1, exact=self.exact, 

223 psi12=psi12, psi2=psi2) 

224 return r 

225 

226# def Inverse3(self, lat1, lon1, lat2, lon2): # PYCHOK outmask 

227# '''Return the distance in C{meter} and the forward and 

228# reverse azimuths (initial and final bearing) in C{degrees}. 

229# 

230# @return: L{Distance3Tuple}C{(distance, initial, final)}. 

231# ''' 

232# r = self.Inverse(lat1, lon1, lat2, lon2) 

233# return Distance3Tuple(r.s12, r.azi12, r.azi12) 

234 

235 @deprecated_method 

236 def Inverse7(self, lat1, lon1, azi12, s12, outmask=Caps.AZIMUTH_DISTANCE_AREA): 

237 '''DEPRECATED, use method L{Rhumb.Inverse8}. 

238 

239 @return: A I{DEPRECATED} L{Rhumb7Tuple}. 

240 ''' 

241 return self.Inverse8(lat1, lon1, azi12, s12, outmask=outmask)._to7Tuple() 

242 

243 def _meanSinXi(self, x, y): # radians 

244 s = _Dlog(cosh(x), cosh(y)) * _Dcosh(x, y) 

245 if self.f: 

246 s += _sincosSeries(False, _gd(x), _gd(y), *self._RA2) * _Dgd(x, y) 

247 return s 

248 

249 @deprecated_method 

250 def orders(self, RAorder=None, TMorder=None): # PYCHOK expected 

251 '''DEPRECATED, use properties C{RAorder} and/or C{TMorder}. 

252 

253 Get and set the I{RAorder} and/or I{TMorder}. 

254 

255 @kwarg RAorder: I{Rhumb Area} order (C{int}, 4, 5, 6, 7 

256 or 8). 

257 @kwarg TMorder: I{Transverse Mercator} order (C{int}, 4, 

258 5, 6, 7 or 8). 

259 

260 @return: L{RhumbOrder2Tuple}C{(RAorder, TMorder)} with 

261 the previous C{RAorder} and C{TMorder} setting. 

262 ''' 

263 t = RhumbOrder2Tuple(self.RAorder, self.TMorder) 

264 if RAorder not in (None, t.RAorder): # PYCHOK attr 

265 self.RAorder = RAorder 

266 if TMorder not in (None, t.TMorder): # PYCHOK attr 

267 self.TMorder = TMorder 

268 return t 

269 

270 @Property_RO 

271 def _RA2(self): 

272 # for WGS84: (0, -0.0005583633519275459, -3.743803759172812e-07, -4.633682270824446e-10, 

273 # RAorder 6: -7.709197397676237e-13, -1.5323287106694307e-15, -3.462875359099873e-18) 

274 m = self.RAorder 

275 return _Xs(_RACoeffs, m, self.ellipsoid, RA=True), m 

276 

277 @Property 

278 def RAorder(self): 

279 '''Get the I{Rhumb Area} order (C{int}, 4, 5, 6, 7 or 8). 

280 ''' 

281 return self._mRA 

282 

283 @RAorder.setter # PYCHOK setter! 

284 def RAorder(self, order): 

285 '''Set the I{Rhumb Area} order (C{int}, 4, 5, 6, 7 or 8). 

286 ''' 

287 n = _Xorder(_RACoeffs, RhumbError, RAorder=order) 

288 if self._mRA != n: 

289 _update_all_rls(self) 

290 self._mRA = n 

291 

292 @Property_RO 

293 def _RhumbLine(self): 

294 '''(INTERNAL) Get this module's C{RhumbLine} class. 

295 ''' 

296 return RhumbLine 

297 

298 def _S12d(self, lon12, psi2, psi1): # degrees 

299 '''(INTERNAL) Compute the area C{S12}. 

300 ''' 

301 r = (self.ellipsoid.areax if self.exact else 

302 self.ellipsoid.area) * lon12 / _720_0 

303 r *= self._meanSinXi(radians(psi2), radians(psi1)) 

304 return r 

305 

306 @Property 

307 def TMorder(self): 

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

309 ''' 

310 return self._mTM 

311 

312 @TMorder.setter # PYCHOK setter! 

313 def TMorder(self, order): 

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

315 

316 @note: Setting C{TMorder} turns property C{exact} off. 

317 ''' 

318 self.exact = self._TMorder(order) 

319 

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

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

322 

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

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

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

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

327 

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

329 ''' 

330 d = dict(ellipsoid=self.ellipsoid, RAorder=self.RAorder, 

331 exact=self.exact, TMorder=self.TMorder) 

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

333 

334 

335class RhumbLine(RhumbLineBase): 

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

337 

338 Class C{RhumbLine} facilitates the determination of points on 

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

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

341 

342 Method C{RhumbLine.Position} returns the location of an other 

343 point at distance C{s12} along and the area C{S12} under the 

344 rhumb line. 

345 

346 Method C{RhumbLine.intersection2} finds the intersection between 

347 two rhumb lines. 

348 

349 Method C{RhumbLine.nearestOn4} computes the nearest point on and 

350 the distance to a rhumb line in different ways. 

351 ''' 

352 _Rhumb = Rhumb # rhumbx.Rhumb 

353 

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

355 '''New C{rhumbx.RhumbLine}. 

356 

357 @arg rhumb: The rhumb reference (C{rhumbx.Rhumb}). 

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

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

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

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

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

363 values specifying the required capabilities. Include 

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

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

366 ''' 

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

368 

369 @Property_RO 

370 def _mu1(self): 

371 '''(INTERNAL) Get the I{rectifying auxiliary} latitude C{mu} (C{degrees}). 

372 ''' 

373 return self.ellipsoid.auxRectifying(self.lat1) 

374 

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

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

377 

378 @arg s12: The distance along this rhumb between its point and 

379 the other point (C{meters}), can be negative. 

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

381 the quantities to be returned. 

382 

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

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

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

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

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

388 

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

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

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

392 

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

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

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

396 carried out in finite terms. 

397 ''' 

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

399 if (outmask & Cs.LATITUDE_LONGITUDE_AREA): 

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

401 a12 = s12 / E._L_90 

402 mu12 = self._calp * a12 

403 mu2, x90 = self._mu22(mu12, self._mu1) 

404 if x90: # PYCHOK no cover 

405 lat2 = E.auxRectifying(mu2, inverse=True) 

406 lon2 = NAN 

407 if (outmask & Cs.AREA): 

408 r.set_(S12=NAN) 

409 else: 

410 psi2 = self._psi1 

411 if self._calp: 

412 lat2 = E.auxRectifying(mu2, inverse=True) 

413 psi12 = R._DRectifying2Isometricd(mu2, 

414 self._mu1) * mu12 

415 lon2 = psi12 * self._salp / self._calp 

416 psi2 += psi12 

417 else: # PYCHOK no cover 

418 lat2 = self.lat1 

419 lon2 = self._salp * s12 / self._r1rad 

420 if (outmask & Cs.AREA): 

421 S12 = R._S12d(lon2, self._psi1, psi2) 

422 r.set_(S12=unsigned0(S12)) # like .gx 

423 if (outmask & Cs.LONGITUDE): 

424 if (outmask & Cs.LONG_UNROLL): 

425 lon2 += self.lon1 

426 else: 

427 lon2 = _norm180(self._lon12 + lon2) 

428 r.set_(azi12=self.azi12, s12=s12, a12=a12) 

429 if (outmask & Cs.LATITUDE): 

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

431 if (outmask & Cs.LONGITUDE): 

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

433 if ((outmask | self._debug) & Cs._DEBUG_DIRECT_LINE): # PYCHOK no cover 

434 r.set_(a=E.a, f=E.f, f1=E.f1, L=E.L, exact=R.exact, 

435 b=E.b, e=E.e, e2=E.e2, k2=R._eF.k2, 

436 calp=self._calp, mu1 =self._mu1, mu12=mu12, 

437 salp=self._salp, psi1=self._psi1, mu2=mu2) 

438 return r 

439 

440 @Property_RO 

441 def _psi1(self): 

442 '''(INTERNAL) Get the I{isometric auxiliary} latitude C{psi} (C{degrees}). 

443 ''' 

444 return self.ellipsoid.auxIsometric(self.lat1) 

445 

446 @property_RO 

447 def RAorder(self): 

448 '''Get this rhumb line's I{Rhumb Area} order (C{int}, 4, 5, 6, 7 or 8). 

449 ''' 

450 return self.rhumb.RAorder 

451 

452 @Property_RO 

453 def _r1rad(self): # PYCHOK no cover 

454 '''(INTERNAL) Get this rhumb line's parallel I{circle radius} (C{meter}). 

455 ''' 

456 return radians(self.ellipsoid.circle4(self.lat1).radius) 

457 

458 

459class RhumbOrder2Tuple(_GTuple): 

460 '''2-Tuple C{(RAorder, TMorder)} with a I{Rhumb Area} and 

461 I{Transverse Mercator} order, both C{int}, DEPRECATED. 

462 ''' 

463 _Names_ = (Rhumb.RAorder.name, Rhumb.TMorder.name) 

464 _Units_ = ( Int, Int) 

465 

466 

467# Use I{Divided Differences} to determine (mu2 - mu1) / (psi2 - psi1) accurately. 

468# Definition: _Df(x,y,d) = (f(x) - f(y)) / (x - y), @see W. M. Kahan & R. J. 

469# Fateman, "Symbolic computation of Divided Differences", SIGSAM Bull. 33(3), 

470# 7-28 (1999). U{ACM<https://DL.ACM.org/doi/pdf/10.1145/334714.334716> and @see 

471# U{UCB<https://www.CS.Berkeley.edu/~fateman/papers/divdiff.pdf>}, Dec 8, 1999. 

472 

473def _Dasinh(x, y): 

474 hx = hypot1(x) 

475 d = x - y 

476 if d: 

477 hx *= y 

478 hy = x * hypot1(y) 

479 t = (d * (x + y) / (hy + hx)) if (x * y) > 0 else (hy - hx) 

480 r = asinh(t) / d 

481 else: 

482 r = _1_0 / hx 

483 return r 

484 

485 

486def _Datan(x, y): 

487 xy = x * y 

488 r = xy + _1_0 

489 d = x - y 

490 if d: # 2 * xy > -1 == 2 * xy + 1 > 0 == xy + r > 0 == xy > -r 

491 r = (atan(d / r) if xy > -r else (atan(x) - atan(y))) / d 

492 else: 

493 r = _1_0 / r 

494 return r 

495 

496 

497def _Dcosh(x, y): 

498 return _Dsincos(x, y, sinh, sinh) 

499 

500 

501def _DeatanhE(x, y, E): # see .albers._Datanhee 

502 # Deatanhe(x, y) = eatanhe((x - y) / (1 - e^2 * x * y)) / (x - y) 

503 e = _1_0 - E.e2 * x * y 

504 if e: # assert not isnear0(e) 

505 d = x - y 

506 e = (E._es_atanh(d / e) / d) if d else (E.e2 / e) 

507 return e 

508 

509 

510def _DfEt(tx, ty, eF): # tangents 

511 # eF = Elliptic(-E.e12) # -E.e2 / (1 - E.e2) 

512 r, x, y, = _1_0, atan(tx), atan(ty) 

513 d = x - y 

514 if (x * y) > 0: 

515 # See U{DLMF<https://DLMF.NIST.gov/19.11>}: 19.11.2 and 19.11.4 

516 # letting theta -> x, phi -> -y, psi -> z 

517 # (E(x) - E(y)) / d = E(z)/d - k2 * sin(x) * sin(y) * sin(z)/d 

518 # tan(z/2) = (sin(x)*Delta(y) - sin(y)*Delta(x)) / (cos(x) + cos(y)) 

519 # = d * Dsin(x,y) * (sin(x) + sin(y))/(cos(x) + cos(y)) / 

520 # (sin(x)*Delta(y) + sin(y)*Delta(x)) 

521 # = t = d * Dt 

522 # sin(z) = 2*t/(1+t^2); cos(z) = (1-t^2)/(1+t^2) 

523 # Alt (this only works for |z| <= pi/2 -- however, this conditions 

524 # holds if x*y > 0): 

525 # sin(z) = d * Dsin(x,y) * (sin(x) + sin(y)) / 

526 # (sin(x)*cos(y)*Delta(y) + sin(y)*cos(x)*Delta(x)) 

527 # cos(z) = sqrt((1-sin(z))*(1+sin(z))) 

528 sx, cx, sy, cy = sincos2_(x, y) 

529 D = (cx + cy) * (eF.fDelta(sy, cy) * sx + 

530 eF.fDelta(sx, cx) * sy) 

531 D = (sx + sy) * _Dsin(x, y) / D 

532 t = D * d 

533 t2 = _1_0 + t**2 

534 D *= _2_0 / t2 

535 s = D * d 

536 if s: 

537 c = (t + _1_0) * (_1_0 - t) / t2 

538 r = eF.fE(s, c, eF.fDelta(s, c)) / s 

539 r = D * (r - eF.k2 * sx * sy) 

540 elif d: 

541 r = (eF.fE(x) - eF.fE(y)) / d 

542 return r 

543 

544 

545def _Dgd(x, y): 

546 return _Datan(sinh(x), sinh(y)) * _Dsinh(x, y) 

547 

548 

549def _Dgdinv(x, y): # x, y are tangents 

550 return _Dasinh(x, y) / _Datan(x, y) 

551 

552 

553def _Dlog(x, y): 

554 d = (x - y) * _0_5 

555 # Changed atanh(t / (x + y)) to asinh(t / (2 * sqrt(x*y))) to 

556 # avoid taking atanh(1) when x is large and y is 1. This also 

557 # fixes bogus results being returned for the area when an endpoint 

558 # is at a pole. N.B. this routine is invoked with positive x 

559 # and y, so the sqrt is always taken of a positive quantity. 

560 return (asinh(d / sqrt(x * y)) / d) if d else (_1_0 / x) 

561 

562 

563def _Dsin(x, y): 

564 return _Dsincos(x, y, sin, cos) 

565 

566 

567def _Dsincos(x, y, sin_, cos_): 

568 r = cos_((x + y) * _0_5) 

569 d = (x - y) * _0_5 

570 if d: 

571 r *= sin_(d) / d 

572 return r 

573 

574 

575def _Dsinh(x, y): 

576 return _Dsincos(x, y, sinh, cosh) 

577 

578 

579def _Dtan(x, y): # PYCHOK no cover 

580 return _Dtant(x - y, tan(x), tan(y)) 

581 

582 

583def _Dtant(dxy, tx, ty): 

584 txy = tx * ty 

585 r = txy + _1_0 

586 if dxy: # 2 * txy > -1 == 2 * txy + 1 > 0 == txy + r > 0 == txy > -r 

587 r = ((tan(dxy) * r) if txy > -r else (tx - ty)) / dxy 

588 return r 

589 

590 

591def _Eaux4(E_aux, mu_psi_x, mu_psi_y): # degrees 

592 # get inverse auxiliary lats in radians and tangents 

593 phix = radians(E_aux(mu_psi_x, inverse=True)) 

594 phiy = radians(E_aux(mu_psi_y, inverse=True)) 

595 return phix, phiy, tan(phix), tan(phiy) 

596 

597 

598def _gd(x): 

599 return atan(sinh(x)) 

600 

601 

602def _sincosSeries(sinp, x, y, C, n): 

603 # N.B. C[] has n+1 elements of which 

604 # C[0] is ignored and n >= 0 

605 # Use Clenshaw summation to evaluate 

606 # m = (g(x) + g(y)) / 2 -- mean value 

607 # s = (g(x) - g(y)) / (x - y) -- average slope 

608 # where 

609 # g(x) = sum(C[j] * SC(2 * j * x), j = 1..n) 

610 # SC = sinp ? sin : cos 

611 # CS = sinp ? cos : sin 

612 # ... 

613 d, _neg = (x - y), neg 

614 sp, cp, sd, cd = sincos2_(x + y, d) 

615 sd = (sd / d) if d else _1_0 

616 s = _neg(sp * sd) # negative 

617 # 2x2 matrices in row-major order 

618 a1 = s * d**2 

619 a2 = s * _4_0 

620 a0 = a3 = _2_0 * cp * cd # m 

621 b2 = b1 = _0_0s(4) 

622 if n > 0: 

623 b1 = C[n], _0_0, _0_0, C[n] 

624 

625 _fsum = _MODS.fsums.fsum1f_ 

626 for j in range(n - 1, 0, -1): # C[0] unused 

627 b1, b2, Cj = b2, b1, C[j] 

628 # b1 = a * b2 - b1 + C[j] * I 

629 m0, m1, m2, m3 = b2 

630 n0, n1, n2, n3 = map(_neg, b1) 

631 b1 = (_fsum(a0 * m0, a1 * m2, n0, Cj), 

632 _fsum(a0 * m1, a1 * m3, n1), 

633 _fsum(a2 * m0, a3 * m2, n2), 

634 _fsum(a2 * m1, a3 * m3, n3, Cj)) 

635 # Here are the full expressions for m and s 

636 # f01, f02, f11, f12 = (0, 0, cd * sp, 2 * sd * cp) if sinp else \ 

637 # (1, 0, cd * cp, -2 * sd * sp) 

638 # m = -b2[1] * f02 + (C[0] - b2[0]) * f01 + b1[0] * f11 + b1[1] * f12 

639 # s = -b2[2] * f01 + (C[0] - b2[3]) * f02 + b1[2] * f11 + b1[3] * f12 

640 cd *= b1[2] 

641 sd *= b1[3] * _2_0 

642 s = _fsum(cd * sp, sd * cp) if sinp else \ 

643 _fsum(cd * cp, _neg(sd * sp), _neg(b2[2])) 

644 return s 

645 

646 

647_RACoeffs = { # Generated by Maxima on 2015-05-15 08:24:04-04:00 

648 4: ( # GEOGRAPHICLIB_RHUMBAREA_ORDER == 4 

649 691, 7860, -20160, 18900, 0, 56700, # R[0]/n^0, polynomial(n), order 4 

650 1772, -5340, 6930, -4725, 14175, # R[1]/n^1, polynomial(n), order 3 

651 -1747, 1590, -630, 4725, # PYCHOK R[2]/n^2, polynomial(n), order 2 

652 104, -31, 315, # R[3]/n^3, polynomial(n), order 1 

653 -41, 420), # PYCHOK R[4]/n^4, polynomial(n), order 0, count = 20 

654 5: ( # GEOGRAPHICLIB_RHUMBAREA_ORDER == 5 

655 -79036, 22803, 259380, -665280, 623700, 0, 1871100, # PYCHOK R[0]/n^0, polynomial(n), order 5 

656 41662, 58476, -176220, 228690, -155925, 467775, # PYCHOK R[1]/n^1, polynomial(n), order 4 

657 18118, -57651, 52470, -20790, 155925, # PYCHOK R[2]/n^2, polynomial(n), order 3 

658 -23011, 17160, -5115, 51975, # PYCHOK R[3]/n^3, polynomial(n), order 2 

659 5480, -1353, 13860, # PYCHOK R[4]/n^4, polynomial(n), order 1 

660 -668, 5775), # PYCHOK R[5]/n^5, polynomial(n), order 0, count = 27 

661 6: ( # GEOGRAPHICLIB_RHUMBAREA_ORDER == 6 

662 128346268, -107884140, 31126095, 354053700, -908107200, 851350500, 0, 2554051500, # R[0]/n^0, polynomial(n), order 6 

663 -114456994, 56868630, 79819740, -240540300, 312161850, -212837625, 638512875, # PYCHOK R[1]/n^1, polynomial(n), order 5 

664 51304574, 24731070, -78693615, 71621550, -28378350, 212837625, # R[2]/n^2, polynomial(n), order 4 

665 1554472, -6282003, 4684680, -1396395, 14189175, # R[3]/n^3, polynomial(n), order 3 

666 -4913956, 3205800, -791505, 8108100, # PYCHOK R[4]/n^4, polynomial(n), order 2 

667 1092376, -234468, 2027025, # R[5]/n^5, polynomial(n), order 1 

668 -313076, 2027025), # PYCHOK R[6]/n^6, polynomial(n), order 0, count = 35 

669 7: ( # GEOGRAPHICLIB_RHUMBAREA_ORDER == 7 

670 -317195588, 385038804, -323652420, 93378285, 1062161100, -2724321600, 2554051500, 0, 7662154500, # PYCHOK R[0]/n^0, polynomial(n), order 7 

671 258618446, -343370982, 170605890, 239459220, -721620900, 936485550, -638512875, 1915538625, # PYCHOK R[1]/n^1, polynomial(n), order 6 

672 -248174686, 153913722, 74193210, -236080845, 214864650, -85135050, 638512875, # PYCHOK R[2]/n^2, polynomial(n), order 5 

673 114450437, 23317080, -94230045, 70270200, -20945925, 212837625, # PYCHOK R[3]/n^3, polynomial(n), order 4 

674 15445736, -103193076, 67321800, -16621605, 170270100, # PYCHOK R[4]/n^4, polynomial(n), order 3 

675 -27766753, 16385640, -3517020, 30405375, # PYCHOK R[4]/n^4, polynomial(n), order 3 

676 4892722, -939228, 6081075, # PYCHOK R[4]/n^4, polynomial(n), order 3 

677 -3189007, 14189175), # PYCHOK R[7]/n^7, polynomial(n), order 0, count = 44 

678 8: ( # GEOGRAPHICLIB_RHUMBAREA_ORDER == 8 

679 71374704821, -161769749880, 196369790040, -165062734200, 47622925350, 541702161000, -1389404016000, 1302566265000, 0, 3907698795000, # R[0]/n^0, polynomial(n), order 8 

680 -13691187484, 65947703730, -87559600410, 43504501950, 61062101100, -184013329500, 238803815250, -162820783125, 488462349375, # PYCHOK R[1]/n^1, polynomial(n), order 7 

681 30802104839, -63284544930, 39247999110, 18919268550, -60200615475, 54790485750, -21709437750, 162820783125, # R[2]/n^2, polynomial(n), order 6 

682 -8934064508, 5836972287, 1189171080, -4805732295, 3583780200, -1068242175, 10854718875, # PYCHOK R[3]/n^3, polynomial(n), order 5 

683 50072287748, 3938662680, -26314234380, 17167059000, -4238509275, 43418875500, # R[4]/n^4, polynomial(n), order 4 

684 359094172, -9912730821, 5849673480, -1255576140, 10854718875, # R[5]/n^5, polynomial(n), order 3 

685 -16053944387, 8733508770, -1676521980, 10854718875, # PYCHOK R[6]/n^6, polynomial(n), order 2 

686 930092876, -162639357, 723647925, # R[7]/n^7, polynomial(n), order 1 

687 -673429061, 1929727800) # PYCHOK R[8]/n^8, polynomial(n), order 0, count = 54 

688} 

689 

690__all__ += _ALL_DOCS(Caps, Rhumb, RhumbLine) 

691 

692if __name__ == '__main__': 

693 

694 from pygeodesy.lazily import printf 

695 

696 def _re(fmt, r3, x3): 

697 e3 = [] 

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

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

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

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

702 

703 # <https://GeographicLib.SourceForge.io/cgi-bin/RhumbSolve> version 2.0 

704 rhumb = Rhumb(exact=True) # WGS84 default 

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

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

707 _re('# JFK NE lat2=%.8f, lon2=%.8f, S12=%.1f', (r.lat2, r.lon2, r.S12), (71.68889988, 0.25551982, 44095641862956.148438)) 

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

709 _re('# JFK-LHR azi12=%.8f, s12=%.3f S12=%.1f', (r.azi12, r.s12, r.S12), (77.76838971, 5771083.383328, 37395209100030.367188)) 

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

711 _re('# JFK-NRT azi12=%.8f, s12=%.3f S12=%.1f', (r.azi12, r.s12, r.S12), (-92.388887981699639, 12782581.0676841792, -63760642939072.492)) 

712 

713# % python3 -m pygeodesy.rhumbx 

714 

715# Rhumb(RAorder=6, 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) 

716 

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

718# JFK-LHR azi12=77.76838971, s12=5771083.383 S12=37395209100030.4 rel errors: 3e-12, 5e-15, 0 

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

720 

721# **) MIT License 

722# 

723# Copyright (C) 2022-2023 -- mrJean1 at Gmail -- All Rights Reserved. 

724# 

725# Permission is hereby granted, free of charge, to any person obtaining a 

726# copy of this software and associated documentation files (the "Software"), 

727# to deal in the Software without restriction, including without limitation 

728# the rights to use, copy, modify, merge, publish, distribute, sublicense, 

729# and/or sell copies of the Software, and to permit persons to whom the 

730# Software is furnished to do so, subject to the following conditions: 

731# 

732# The above copyright notice and this permission notice shall be included 

733# in all copies or substantial portions of the Software. 

734# 

735# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS 

736# OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, 

737# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL 

738# THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR 

739# OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, 

740# ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR 

741# OTHER DEALINGS IN THE SOFTWARE.