Coverage for pygeodesy/rhumbx.py: 97%

458 statements  

« prev     ^ index     » next       coverage.py v7.2.2, created at 2023-06-07 08:37 -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>}. 

7 

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

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

10 

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

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

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

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

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

16 

17Copyright (C) U{Charles Karney<mailto:Charles@Karney.com>} (2014-2022) 

18and licensed under the MIT/X11 License. For more information, see the 

19U{GeographicLib<https://GeographicLib.SourceForge.io>} documentation. 

20''' 

21# make sure int/int division yields float quotient 

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

23 

24from pygeodesy.basics import copysign0, neg, _xinstanceof, _zip 

25from pygeodesy.constants import INT0, _EPSqrt as _TOL, NAN, PI_2, isnan, _0_0s, \ 

26 _0_0, _0_5, _1_0, _2_0, _4_0, _90_0, _180_0, _720_0 

27# from pygeodesy.datums import _spherical_datum # in Rhumb.ellipsoid.setter 

28from pygeodesy.errors import IntersectionError, itemsorted, _ValueError, \ 

29 _xdatum, _xkwds 

30# from pygeodesy.etm import ExactTransverseMercator # _MODS in ._RhumbLine.xTM 

31from pygeodesy.fmath import euclid, favg, hypot, hypot1 

32# from pygeodesy.fsums import fsum1f_ # _MODS 

33from pygeodesy.interns import NN, _azi12_, _coincident_, _COMMASPACE_, \ 

34 _intersection_, _lat1_, _lat2_, _lon1_, _lon2_, \ 

35 _no_, _s12_, _S12_, _UNDER 

36from pygeodesy.karney import _a12_, _atan2d, Caps, _CapsBase as _RhumbBase, \ 

37 _diff182, Direct9Tuple, _EWGS84, _fix90, GDict, \ 

38 _GTuple, Inverse10Tuple, _norm180 

39from pygeodesy.ktm import KTransverseMercator, _Xorder, _Xs, \ 

40 _AlpCoeffs, _BetCoeffs # PYCHOK used! 

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

42from pygeodesy.namedTuples import Distance2Tuple, LatLon2Tuple, NearestOn4Tuple 

43from pygeodesy.props import deprecated_method, Property, Property_RO, property_RO, \ 

44 _update_all 

45from pygeodesy.streprs import Fmt, pairs, unstr 

46from pygeodesy.units import Bearing as _Azi, Degrees as _Deg, Int, Lat, Lon, \ 

47 Meter as _M, Meter2 as _M2 

48from pygeodesy.utily import sincos2_, sincos2d, _unrollon, _Wrap 

49from pygeodesy.vector3d import _intersect3d3, Vector3d # in .intersection2 below 

50 

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

52 

53__all__ = _ALL_LAZY.rhumbx 

54__version__ = '23.05.26' 

55 

56_rls = [] # instances of C{RbumbLine} to be updated 

57_TRIPS = 65 # .intersection2, 18+ 

58 

59 

60class _Lat(Lat): 

61 '''(INTERNAL) Latitude B{C{lat}}. 

62 ''' 

63 def __init__(self, *lat, **Error_name): 

64 kwds = _xkwds(Error_name, clip=0, Error=RhumbError) 

65 Lat.__new__(_Lat, *lat, **kwds) 

66 

67 

68class _Lon(Lon): 

69 '''(INTERNAL) Longitude B{C{lon}}. 

70 ''' 

71 def __init__(self, *lon, **Error_name): 

72 kwds = _xkwds(Error_name, clip=0, Error=RhumbError) 

73 Lon.__new__(_Lon, *lon, **kwds) 

74 

75 

76def _update_all_rls(r): 

77 '''(INTERNAL) Zap cached/memoized C{Property[_RO]}s 

78 of any L{RhumbLine} instances tied to the given 

79 L{Rhumb} instance B{C{r}}. 

80 ''' 

81 _xinstanceof(r, Rhumb) 

82 _update_all(r) 

83 for rl in _rls: # PYCHOK use weakref? 

84 if rl._rhumb is r: 

85 _update_all(rl) 

86 

87 

88class Rhumb(_RhumbBase): 

89 '''Class to solve of the I{direct} and I{inverse rhumb} problems, accurately. 

90 

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

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

93 ''' 

94 _E = _EWGS84 

95 _exact = True 

96 _mRA = 6 

97 _mTM = 6 

98 

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

100 '''New L{Rhumb}. 

101 

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

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

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

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

106 a C{scalar}, ignored otherwise. 

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

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

109 series expansion (C{bool}), see also property C{exact}. 

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

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

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

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

114 

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

116 ''' 

117 if f is not None: 

118 self.ellipsoid = a_earth, f 

119 elif a_earth not in (_EWGS84, None): 

120 self.ellipsoid = a_earth 

121 if not exact: 

122 self._exact = False 

123 if name: 

124 self.name = name 

125 if RA_TMorder: 

126 self.orders(**RA_TMorder) 

127 

128 @Property_RO 

129 def _A2(self): # Conformal2RectifyingCoeffs 

130 m = self.TMorder 

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

132 

133 @Property_RO 

134 def _B2(self): # Rectifying2ConformalCoeffs 

135 m = self.TMorder 

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

137 

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

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

140 

141 def _Direct(self, ll1, azi12, s12, **outmask): 

142 '''(INTERNAL) Short-cut version, see .latlonBase. 

143 ''' 

144 return self.Direct(ll1.lat, ll1.lon, azi12, s12, **outmask) 

145 

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

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

148 

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

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

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

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

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

154 

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

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

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

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

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

160 

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

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

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

164 

165 If the given point is a pole, the cosine of its latitude is 

166 taken to be C{epsilon}**-2 (where C{epsilon} is 2.0**-52. 

167 This position is extremely close to the actual pole and 

168 allows the calculation to be carried out in finite terms. 

169 ''' 

170 rl = _RhumbLine(self, lat1, lon1, azi12, caps=Caps.LINE_OFF, 

171 name=self.name) 

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

173 

174 @deprecated_method 

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

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

177 

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

179 ''' 

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

181 

182 def Direct8(self, lat1, lon1, azi12, s12, outmask=Caps.LATITUDE_LONGITUDE_AREA): 

183 '''Like method L{Rhumb.Direct} but returning a L{Rhumb8Tuple} with area C{S12}. 

184 ''' 

185 return self.Direct(lat1, lon1, azi12, s12, outmask=outmask).toRhumb8Tuple() 

186 

187 def _DirectLine(self, ll1, azi12, **name_caps): 

188 '''(INTERNAL) Short-cut version, see .latlonBase. 

189 ''' 

190 return self.DirectLine(ll1.lat, ll1.lon, azi12, **name_caps) 

191 

192 def DirectLine(self, lat1, lon1, azi12, name=NN, **caps): # caps=Caps.STANDARD 

193 '''Define a L{RhumbLine} in terms of the I{direct} rhumb problem. 

194 

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

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

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

198 @kwarg caps: Optional C{caps}, see L{RhumbLine} C{B{caps}}. 

199 

200 @return: A L{RhumbLine} instance and invoke its method 

201 L{RhumbLine.Position} to compute each point. 

202 

203 @note: Updates to this rhumb are reflected in the returned 

204 rhumb line. 

205 ''' 

206 return RhumbLine(self, lat1=lat1, lon1=lon1, azi12=azi12, 

207 name=name or self.name, **caps) 

208 

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

210 E = self.ellipsoid 

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

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

213 

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

215 if self.exact: 

216 E = self.ellipsoid 

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

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

219 r = self._DRectifyingt( tphix, tphiy, t) / \ 

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

221 else: 

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

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

224 return r 

225 

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

227 E = self.ellipsoid 

228 tbetx = E.f1 * tphix 

229 tbety = E.f1 * tphiy 

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

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

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

233 

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

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

236 

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

238 E = self.ellipsoid 

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

240 if self.exact: 

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

242 r = self._DIsometrict(phix, phiy, tphix, tphiy, t) / \ 

243 self._DRectifyingt( tphix, tphiy, t) 

244 else: 

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

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

247 return r 

248 

249 @Property_RO 

250 def _eF(self): 

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

252 ''' 

253 # .k2 = 0.006739496742276434 

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

255 

256 @Property 

257 def ellipsoid(self): 

258 '''Get this rhumb's ellipsoid (L{Ellipsoid}). 

259 ''' 

260 return self._E 

261 

262 @ellipsoid.setter # PYCHOK setter! 

263 def ellipsoid(self, a_earth_f): 

264 '''Set this rhumb's ellipsoid (L{Ellipsoid}, L{Ellipsoid2}, L{Datum}, 

265 L{a_f2Tuple}, 2-tuple C{(a, f)}) or the (equatorial) radius (C{scalar}). 

266 ''' 

267 E = _MODS.datums._spherical_datum(a_earth_f, Error=RhumbError).ellipsoid 

268 if self._E != E: 

269 _update_all_rls(self) 

270 self._E = E 

271 

272 @property_RO 

273 def equatoradius(self): 

274 '''Get the C{ellipsoid}'s equatorial radius, semi-axis (C{meter}). 

275 ''' 

276 return self.ellipsoid.a 

277 

278 a = equatoradius 

279 

280 @Property 

281 def exact(self): 

282 '''Get the I{exact} option (C{bool}). 

283 ''' 

284 return self._exact 

285 

286 @exact.setter # PYCHOK setter! 

287 def exact(self, exact): 

288 '''Set the I{exact} option (C{bool}). If C{True}, use I{exact} rhumb 

289 calculations, if C{False} results are less precise for more oblate 

290 or more prolate ellipsoids with M{abs(flattening) > 0.01} (C{bool}). 

291 

292 @see: Option U{B{-s}<https://GeographicLib.SourceForge.io/C++/doc/RhumbSolve.1.html>} 

293 and U{ACCURACY<https://GeographicLib.SourceForge.io/C++/doc/RhumbSolve.1.html#ACCURACY>}. 

294 ''' 

295 x = bool(exact) 

296 if self._exact != x: 

297 _update_all_rls(self) 

298 self._exact = x 

299 

300 def flattening(self): 

301 '''Get the C{ellipsoid}'s flattening (C{float}). 

302 ''' 

303 return self.ellipsoid.f 

304 

305 f = flattening 

306 

307 def _Inverse(self, ll1, ll2, wrap, **outmask): 

308 '''(INTERNAL) Short-cut version, see .latlonBase. 

309 ''' 

310 if wrap: 

311 ll2 = _unrollon(ll1, _Wrap.point(ll2)) 

312 return self.Inverse(ll1.lat, ll1.lon, ll2.lat, ll2.lon, **outmask) 

313 

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

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

316 

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

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

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

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

321 

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

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

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

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

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

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

328 

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

330 on opposite meridians, there are two shortest rhumb lines 

331 and the East-going one is chosen. 

332 

333 If either point is a pole, the cosine of its latitude is 

334 taken to be C{epsilon}**-2 (where C{epsilon} is 2.0**-52). 

335 This position is extremely close to the actual pole and 

336 allows the calculation to be carried out in finite terms. 

337 ''' 

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

339 if (outmask & Cs.AZIMUTH_DISTANCE_AREA): 

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

341 E = self.ellipsoid 

342 psi1 = E.auxIsometric(lat1) 

343 psi2 = E.auxIsometric(lat2) 

344 psi12 = psi2 - psi1 

345 lon12, _ = _diff182(lon1, lon2) 

346 if (outmask & Cs.AZIMUTH): 

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

348 if (outmask & Cs.DISTANCE): 

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

350 s12 = a12 * E._L_90 

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

352 if (outmask & Cs.AREA): 

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

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

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

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

357 lon12=lon12, psi1=psi1, exact=self.exact, 

358 psi12=psi12, psi2=psi2) 

359 return r 

360 

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

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

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

364# 

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

366# ''' 

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

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

369 

370 @deprecated_method 

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

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

373 

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

375 ''' 

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

377 

378 def Inverse8(self, lat1, lon1, azi12, s12, outmask=Caps.AZIMUTH_DISTANCE_AREA): 

379 '''Like method L{Rhumb.Inverse} but returning a L{Rhumb8Tuple} with area C{S12}. 

380 ''' 

381 return self.Inverse(lat1, lon1, azi12, s12, outmask=outmask).toRhumb8Tuple() 

382 

383 def _InverseLine(self, ll1, ll2, wrap, **name_caps): 

384 '''(INTERNAL) Short-cut version, see .latlonBase. 

385 ''' 

386 if wrap: 

387 ll2 = _unrollon(ll1, _Wrap.point(ll2)) 

388 return self.InverseLine(ll1.lat, ll1.lon, ll2.lat, ll2.lon, **name_caps) 

389 

390 def InverseLine(self, lat1, lon1, lat2, lon2, name=NN, **caps): # caps=Caps.STANDARD 

391 '''Define a L{RhumbLine} in terms of the I{inverse} rhumb problem. 

392 

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

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

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

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

397 @kwarg caps: Optional C{caps}, see L{RhumbLine} C{B{caps}}. 

398 

399 @return: A L{RhumbLine} instance and invoke its method 

400 L{RhumbLine.Position} to compute each point. 

401 

402 @note: Updates to this rhumb are reflected in the returned 

403 rhumb line. 

404 ''' 

405 r = self.Inverse(lat1, lon1, lat2, lon2, outmask=Caps.AZIMUTH) 

406 return RhumbLine(self, lat1=lat1, lon1=lon1, azi12=r.azi12, 

407 name=name or self.name, **caps) 

408 

409 Line = DirectLine # synonyms 

410 

411 def _MeanSinXi(self, x, y): # radians 

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

413 if self.f: 

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

415 return s 

416 

417 def orders(self, RAorder=None, TMorder=None): 

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

419 

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

421 or 8). 

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

423 5, 6, 7 or 8). 

424 

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

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

427 ''' 

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

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

430 self.RAorder = RAorder 

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

432 self.TMorder = TMorder 

433 return t 

434 

435 @Property_RO 

436 def _RA2(self): 

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

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

439 m = self.RAorder 

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

441 

442 @Property 

443 def RAorder(self): 

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

445 ''' 

446 return self._mRA 

447 

448 @RAorder.setter # PYCHOK setter! 

449 def RAorder(self, order): 

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

451 ''' 

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

453 if self._mRA != n: 

454 _update_all_rls(self) 

455 self._mRA = n 

456 

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

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

459 ''' 

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

461 self.ellipsoid.area) * lon12 / _720_0 

462 r *= self._MeanSinXi(radians(psi2), radians(psi1)) 

463 return r 

464 

465 @Property 

466 def TMorder(self): 

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

468 ''' 

469 return self._mTM 

470 

471 @TMorder.setter # PYCHOK setter! 

472 def TMorder(self, order): 

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

474 

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

476 ''' 

477 n = _Xorder(_AlpCoeffs, RhumbError, TMorder=order) 

478 if self._mTM != n: 

479 _update_all_rls(self) 

480 self._mTM = n 

481 self.exact = False 

482 

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

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

485 

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

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

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

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

490 

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

492 ''' 

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

494 exact=self.exact, TMorder=self.TMorder) 

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

496 

497 

498class RhumbError(_ValueError): 

499 '''Raised for a L{Rhumb} or L{RhumbLine} issue. 

500 ''' 

501 pass 

502 

503 

504class _RhumbLine(_RhumbBase): 

505 '''(INTERNAL) Class L{RhumbLine} 

506 ''' 

507 _azi12 = _0_0 

508# _lat1 = _0_0 

509# _lon1 = _0_0 

510 _salp = _0_0 

511 _calp = _1_0 

512 _rhumb = None # L{Rhumb} instance 

513 

514 def __init__(self, rhumb, lat1, lon1, azi12, caps=0, name=NN): # case=Caps.? 

515 '''New C{RhumbLine}. 

516 ''' 

517 _xinstanceof(Rhumb, rhumb=rhumb) 

518 self._lat1 = _Lat(lat1=_fix90(lat1)) 

519 self._lon1 = _Lon(lon1= lon1) 

520 self._debug |= (caps | rhumb._debug) & Caps._DEBUG_DIRECT_LINE 

521 if azi12: # non-zero 

522 self.azi12 = azi12 

523 self._caps = caps 

524 if not (caps & Caps.LINE_OFF): 

525 _rls.append(self) 

526 n = name or rhumb.name 

527 if n: 

528 self.name=n 

529 self._rhumb = rhumb # last 

530 

531 def __del__(self): # XXX use weakref? 

532 if _rls: # may be empty or None 

533 try: # PYCHOK no cover 

534 _rls.remove(self) 

535 except (TypeError, ValueError): 

536 pass 

537 self._rhumb = None 

538 # _update_all(self) # throws TypeError during Python 2 cleanup 

539 

540 @Property 

541 def azi12(self): 

542 '''Get this rhumb line's I{azimuth} (compass C{degrees}). 

543 ''' 

544 return self._azi12 

545 

546 @azi12.setter # PYCHOK setter! 

547 def azi12(self, azi12): 

548 '''Set this rhumb line's I{azimuth} (compass C{degrees}). 

549 ''' 

550 z = _norm180(azi12) 

551 if self._azi12 != z: 

552 if self._rhumb: 

553 _update_all(self) 

554 self._azi12 = z 

555 self._salp, self._calp = sincos2d(z) # no NEG0 

556 

557 def distance2(self, lat, lon): 

558 '''Return the distance from and (initial) bearing at the given 

559 point to this rhumb line's start point. 

560 

561 @arg lat: Latitude of the point (C{degrees}). 

562 @arg lon: Longitude of the points (C{degrees}). 

563 

564 @return: A L{Distance2Tuple}C{(distance, initial)} with the C{distance} 

565 in C{meter} and C{initial} bearing in C{degrees}. 

566 

567 @see: Methods L{RhumbLine.intersection2} and L{RhumbLine.nearestOn4}. 

568 ''' 

569 r = self.rhumb.Inverse(self.lat1, self.lon1, lat, lon) 

570# outmask=Caps.AZIMUTH_DISTANCE) 

571 return Distance2Tuple(r.s12, r.azi12) 

572 

573 @Property_RO 

574 def ellipsoid(self): 

575 '''Get this rhumb line's ellipsoid (L{Ellipsoid}). 

576 ''' 

577 return self.rhumb.ellipsoid 

578 

579 @property_RO 

580 def exact(self): 

581 '''Get this rhumb line's I{exact} option (C{bool}). 

582 ''' 

583 return self.rhumb.exact 

584 

585 def intersection2(self, other, tol=_TOL, **eps): 

586 '''I{Iteratively} find the intersection of this and an other rhumb line. 

587 

588 @arg other: The other rhumb line (L{RhumbLine}). 

589 @kwarg tol: Tolerance for longitudinal convergence (C{degrees}). 

590 @kwarg eps: Tolerance for L{intersection3d3} (C{EPS}). 

591 

592 @return: A L{LatLon2Tuple}{(lat, lon)} with the C{lat}- and 

593 C{lon}gitude of the intersection point. 

594 

595 @raise IntersectionError: No convergence for this B{C{tol}} or 

596 no intersection for an other reason. 

597 

598 @see: Methods L{RhumbLine.distance2} and L{RhumbLine.nearestOn4} 

599 and function L{pygeodesy.intersection3d3}. 

600 

601 @note: Each iteration involves a round trip to this rhumb line's 

602 L{ExactTransverseMercator} or L{KTransverseMercator} 

603 projection and invoking function L{intersection3d3} in 

604 that domain. 

605 ''' 

606 _xinstanceof(other, _RhumbLine) 

607 _xdatum(self.rhumb, other.rhumb, Error=RhumbError) 

608 try: 

609 if other is self: 

610 raise ValueError(_coincident_) 

611 # make globals and invariants locals 

612 _diff = euclid # approximate length 

613 _i3d3 = _intersect3d3 # NOT .vector3d.intersection3d3 

614 _LL2T = LatLon2Tuple 

615 _xTMr = self.xTM.reverse # ellipsoidal or spherical 

616 _s_3d, s_az = self._xTM3d, self.azi12 

617 _o_3d, o_az = other._xTM3d, other.azi12 

618 # use halfway point as initial estimate 

619 p = _LL2T(favg(self.lat1, other.lat1), 

620 favg(self.lon1, other.lon1)) 

621 for i in range(1, _TRIPS): 

622 v = _i3d3(_s_3d(p), s_az, # point + bearing 

623 _o_3d(p), o_az, useZ=False, **eps)[0] 

624 t = _xTMr(v.x, v.y, lon0=p.lon) # PYCHOK Reverse4Tuple 

625 d = _diff(t.lon - p.lon, t.lat) # PYCHOK t.lat + p.lat - p.lat 

626 p = _LL2T(t.lat + p.lat, t.lon) # PYCHOK t.lon + p.lon = lon0 

627 if d < tol: 

628 return _LL2T(p.lat, p.lon, iteration=i, # PYCHOK p... 

629 name=self.intersection2.__name__) 

630 except Exception as x: 

631 raise IntersectionError(_no_(_intersection_), cause=x) 

632 t = unstr(self.intersection2, tol=tol, **eps) 

633 raise IntersectionError(Fmt.no_convergence(d), txt=t) 

634 

635 @property_RO 

636 def lat1(self): 

637 '''Get this rhumb line's latitude (C{degrees90}). 

638 ''' 

639 return self._lat1 

640 

641 @property_RO 

642 def lon1(self): 

643 '''Get this rhumb line's longitude (C{degrees180}). 

644 ''' 

645 return self._lon1 

646 

647 @Property_RO 

648 def latlon1(self): 

649 '''Get this rhumb line's lat- and longitude (L{LatLon2Tuple}C{(lat, lon)}). 

650 ''' 

651 return LatLon2Tuple(self.lat1, self.lon1) 

652 

653 @Property_RO 

654 def _mu1(self): 

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

656 ''' 

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

658 

659 def nearestOn4(self, lat, lon, tol=_TOL, **eps): 

660 '''I{Iteratively} locate the point on this rhumb line nearest to 

661 the given point. 

662 

663 @arg lat: Latitude of the point (C{degrees}). 

664 @arg lon: Longitude of the point (C{degrees}). 

665 @kwarg tol: Longitudinal convergence tolerance (C{degrees}). 

666 @kwarg eps: Tolerance for L{intersection3d3} (C{EPS}). 

667 

668 @return: A L{NearestOn4Tuple}C{(lat, lon, distance, normal)} with 

669 the C{lat}- and C{lon}gitude of the nearest point on and 

670 the C{distance} in C{meter} to this rhumb line and with the 

671 azimuth of the C{normal}, perpendicular to this rhumb line. 

672 

673 @raise IntersectionError: No convergence for this B{C{eps}} or 

674 no intersection for an other reason. 

675 

676 @see: Methods L{RhumbLine.distance2} and L{RhumbLine.intersection2} 

677 and function L{intersection3d3}. 

678 ''' 

679 z = _norm180(self.azi12 + _90_0) # perpendicular 

680 r = _RhumbLine(self.rhumb, lat, lon, z, caps=Caps.LINE_OFF) 

681 p = self.intersection2(r, tol=tol, **eps) 

682 t = r.distance2(p.lat, p.lon) 

683 return NearestOn4Tuple(p.lat, p.lon, t.distance, z, 

684 iteration=p.iteration) 

685 

686 @Property_RO 

687 def _psi1(self): 

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

689 ''' 

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

691 

692 @property_RO 

693 def RAorder(self): 

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

695 ''' 

696 return self.rhumb.RAorder 

697 

698 @Property_RO 

699 def _r1rad(self): # PYCHOK no cover 

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

701 ''' 

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

703 

704 @Property_RO 

705 def rhumb(self): 

706 '''Get this rhumb line's rhumb (L{Rhumb}). 

707 ''' 

708 return self._rhumb 

709 

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

711 '''Compute a position at a distance on this rhumb line. 

712 

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

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

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

716 the quantities to be returned. 

717 

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

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

720 C{lon2} of the other point in C{degrees}, the rhumb angle 

721 C{a12} between both points in C{degrees} and the area C{S12} 

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

723 

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

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

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

727 

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

729 taken to be C{epsilon}**-2 (where C{epsilon} is 2**-52). 

730 This position is extremely close to the actual pole and 

731 allows the calculation to be carried out in finite terms. 

732 ''' 

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

734 if (outmask & Cs.LATITUDE_LONGITUDE_AREA): 

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

736 mu12 = s12 * self._calp / E._L_90 

737 mu2 = mu12 + self._mu1 

738 if fabs(mu2) > 90: # PYCHOK no cover 

739 mu2 = _norm180(mu2) # reduce to [-180, 180) 

740 if fabs(mu2) > 90: # point on anti-meridian 

741 mu2 = _norm180(_180_0 - mu2) 

742 lat2x = E.auxRectifying(mu2, inverse=True) 

743 lon2x = NAN 

744 if (outmask & Cs.AREA): 

745 r.set_(S12=NAN) 

746 else: 

747 psi2 = self._psi1 

748 if self._calp: 

749 lat2x = E.auxRectifying(mu2, inverse=True) 

750 psi12 = R._DRectifying2Isometricd(mu2, 

751 self._mu1) * mu12 

752 lon2x = psi12 * self._salp / self._calp 

753 psi2 += psi12 

754 else: # PYCHOK no cover 

755 lat2x = self.lat1 

756 lon2x = self._salp * s12 / self._r1rad 

757 if (outmask & Cs.AREA): 

758 r.set_(S12=R._S12d(lon2x, self._psi1, psi2)) 

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

760 if (outmask & Cs.LATITUDE): 

761 r.set_(lat2=lat2x, lat1=self.lat1) 

762 if (outmask & Cs.LONGITUDE): 

763 if (outmask & Cs.LONG_UNROLL) and not isnan(lat2x): 

764 lon2x += self.lon1 

765 else: 

766 lon2x = _norm180(_norm180(self.lon1) + lon2x) 

767 r.set_(lon2=lon2x, lon1=self.lon1) 

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

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

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

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

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

773 return r 

774 

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

776 '''Return this C{RhumbLine} as string. 

777 

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

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

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

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

782 

783 @return: C{RhumbLine} (C{str}). 

784 ''' 

785 d = dict(rhumb=self.rhumb, lat1=self.lat1, lon1=self.lon1, 

786 azi12=self.azi12, exact=self.exact, 

787 TMorder=self.TMorder, xTM=self.xTM) 

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

789 

790 @property_RO 

791 def TMorder(self): 

792 '''Get this rhumb line's I{Transverse Mercator} order (C{int}, 4, 5, 6, 7 or 8). 

793 ''' 

794 return self.rhumb.TMorder 

795 

796 @Property_RO 

797 def xTM(self): 

798 '''Get this rhumb line's I{Transverse Mercator} projection (L{ExactTransverseMercator} 

799 if I{exact} and I{ellipsoidal}, otherwise L{KTransverseMercator}). 

800 ''' 

801 E = self.ellipsoid 

802 # ExactTransverseMercator doesn't handle spherical earth models 

803 return _MODS.etm.ExactTransverseMercator(E) if self.exact and E.isEllipsoidal else \ 

804 KTransverseMercator(E, TMorder=self.TMorder) 

805 

806 def _xTM3d(self, latlon0, z=INT0, V3d=Vector3d): 

807 '''(INTERNAL) C{xTM.forward} this C{latlon1} to C{V3d} with B{C{latlon0}} 

808 as current intersection estimate and central meridian. 

809 ''' 

810 t = self.xTM.forward(self.lat1 - latlon0.lat, self.lon1, lon0=latlon0.lon) 

811 return V3d(t.easting, t.northing, z) 

812 

813 

814class RhumbLine(_RhumbLine): 

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

816 

817 Class L{RhumbLine} facilitates the determination of points on 

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

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

820 

821 Method L{RhumbLine.Position} returns the location of an other 

822 point and optionally the distance C{s12} along the corresponding 

823 area C{S12} under the rhumb line. 

824 

825 Method L{RhumbLine.intersection2} finds the intersection between 

826 two rhumb lines. 

827 

828 Method L{RhumbLine.nearestOn4} computes the nearest point on and 

829 its distance to a rhumb line. 

830 ''' 

831 def __init__(self, rhumb, lat1=0, lon1=0, azi12=None, caps=0, name=NN): # case=Caps.? 

832 '''New L{RhumbLine}. 

833 

834 @arg rhumb: The rhumb reference (L{Rhumb}). 

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

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

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

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

839 the capabilities. Include C{Caps.LINE_OFF} if 

840 updates to B{C{rhumb}} should I{not} be reflected 

841 in this rhumb line. 

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

843 ''' 

844 if (caps & Caps.LINE_OFF): # copy to avoid updates 

845 rhumb = rhumb.copy(deep=False, name=_UNDER(rhumb.name)) 

846 _RhumbLine.__init__(self, rhumb, lat1, lon1, azi12, caps=caps, name=name) 

847 

848 

849class RhumbOrder2Tuple(_GTuple): 

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

851 I{Transverse Mercator} order, both C{int}. 

852 ''' 

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

854 _Units_ = ( Int, Int) 

855 

856 

857class Rhumb8Tuple(_GTuple): 

858 '''8-Tuple C{(lat1, lon1, lat2, lon2, azi12, s12, S12, a12)} with lat- C{lat1}, 

859 C{lat2} and longitudes C{lon1}, C{lon2} of both points, the azimuth of the 

860 rhumb line C{azi12}, the distance C{s12}, the area C{S12} under the rhumb 

861 line and the angular distance C{a12} between both points. 

862 ''' 

863 _Names_ = (_lat1_, _lon1_, _lat2_, _lon2_, _azi12_, _s12_, _S12_, _a12_) 

864 _Units_ = (_Lat, _Lon, _Lat, _Lon, _Azi, _M, _M2, _Deg) 

865 

866 def toDirect9Tuple(self, dflt=NAN, **a12_azi1_azi2_m12_M12_M21): 

867 '''Convert this L{Rhumb8Tuple} result to a 9-tuple, like I{Karney}'s 

868 method C{geographiclib.geodesic.Geodesic._GenDirect}. 

869 

870 @kwarg dflt: Default value for missing items (C{any}). 

871 @kwarg a12_azi1_azi2_m12_M12_M21: Optional keyword arguments 

872 to specify or override L{Inverse10Tuple} items. 

873 

874 @return: L{Direct9Tuple}C{(a12, lat2, lon2, azi2, s12, 

875 m12, M12, M21, S12)} 

876 ''' 

877 d = dict(azi1=self.azi12, M12=_1_0, m12=self.s12, # PYCHOK attr 

878 azi2=self.azi12, M21=_1_0) # PYCHOK attr 

879 if a12_azi1_azi2_m12_M12_M21: 

880 d.update(a12_azi1_azi2_m12_M12_M21) 

881 return self._toTuple(Direct9Tuple, dflt, d) 

882 

883 def toInverse10Tuple(self, dflt=NAN, **a12_m12_M12_M21_salp1_calp1_salp2_calp2): 

884 '''Convert this L{Rhumb8Tuple} to a 10-tuple, like I{Karney}'s 

885 method C{geographiclib.geodesic.Geodesic._GenInverse}. 

886 

887 @kwarg dflt: Default value for missing items (C{any}). 

888 @kwarg a12_m12_M12_M21_salp1_calp1_salp2_calp2: Optional keyword 

889 arguments to specify or override L{Inverse10Tuple} items. 

890 

891 @return: L{Inverse10Tuple}C{(a12, s12, salp1, calp1, salp2, calp2, 

892 m12, M12, M21, S12)}. 

893 ''' 

894 s, c = sincos2d(self.azi12) # PYCHOK attr 

895 d = dict(salp1=s, calp1=c, M12=_1_0, m12=self.s12, # PYCHOK attr 

896 salp2=s, calp2=c, M21=_1_0) 

897 if a12_m12_M12_M21_salp1_calp1_salp2_calp2: 

898 d.update(a12_m12_M12_M21_salp1_calp1_salp2_calp2) 

899 return self._toTuple(Inverse10Tuple, dflt, d) 

900 

901 def _toTuple(self, nTuple, dflt, updates={}): 

902 '''(INTERNAL) Convert this C{Rhumb8Tuple} to an B{C{nTuple}}. 

903 ''' 

904 _g = self.toGDict(**updates).get 

905 t = tuple(_g(n, dflt) for n in nTuple._Names_) 

906 return nTuple(t, name=self.name) 

907 

908 @deprecated_method 

909 def _to7Tuple(self): 

910 '''DEPRECATED, do not use! 

911 ''' 

912 return _MODS.deprecated.Rhumb7Tuple(self[:-1]) 

913 

914 

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

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

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

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

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

920 

921def _Dasinh(x, y): 

922 hx = hypot1(x) 

923 d = x - y 

924 if d: 

925 hx *= y 

926 hy = x * hypot1(y) 

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

928 r = asinh(t) / d 

929 else: 

930 r = _1_0 / hx 

931 return r 

932 

933 

934def _Datan(x, y): 

935 xy = x * y 

936 r = xy + _1_0 

937 d = x - y 

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

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

940 else: 

941 r = _1_0 / r 

942 return r 

943 

944 

945def _Dcosh(x, y): 

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

947 

948 

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

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

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

952 # assert not isnear0(e) 

953 d = x - y 

954 return (E._es_atanh(d / e) / d) if d else (E.e2 / e) 

955 

956 

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

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

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

960 d = x - y 

961 if (x * y) > 0: 

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

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

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

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

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

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

968 # = t = d * Dt 

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

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

971 # holds if x*y > 0): 

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

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

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

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

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

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

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

979 t = D * d 

980 t2 = t**2 + _1_0 

981 D *= _2_0 / t2 

982 s = D * d 

983 if s: 

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

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

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

987 elif d: 

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

989 return r 

990 

991 

992def _Dgd(x, y): 

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

994 

995 

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

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

998 

999 

1000def _Dlog(x, y): 

1001 d = (x - y) * _0_5 

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

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

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

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

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

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

1008 

1009 

1010def _Dsin(x, y): 

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

1012 

1013 

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

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

1016 d = (x - y) * _0_5 

1017 if d: 

1018 r *= sin_(d) / d 

1019 return r 

1020 

1021 

1022def _Dsinh(x, y): 

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

1024 

1025 

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

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

1028 

1029 

1030def _Dtant(dxy, tx, ty): 

1031 txy = tx * ty 

1032 r = txy + _1_0 

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

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

1035 return r 

1036 

1037 

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

1039 # get inverse auxiliary lats in radians and tangents 

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

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

1042 return phix, phiy, tan(phix), tan(phiy) 

1043 

1044 

1045def _gd(x): 

1046 return atan(sinh(x)) 

1047 

1048 

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

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

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

1052 # Use Clenshaw summation to evaluate 

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

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

1055 # where 

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

1057 # SC = sinp ? sin : cos 

1058 # CS = sinp ? cos : sin 

1059 # ... 

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

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

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

1063 s = _neg(sp * sd) # negative 

1064 # 2x2 matrices in row-major order 

1065 a1 = s * d**2 

1066 a2 = s * _4_0 

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

1068 b2 = b1 = _0_0s(4) 

1069 if n > 0: 

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

1071 

1072 _fsum1f_ = _MODS.fsums.fsum1f_ 

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

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

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

1076 m0, m1, m2, m3 = b2 

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

1078 b1 = (_fsum1f_(a0 * m0, a1 * m2, n0, Cj), 

1079 _fsum1f_(a0 * m1, a1 * m3, n1), 

1080 _fsum1f_(a2 * m0, a3 * m2, n2), 

1081 _fsum1f_(a2 * m1, a3 * m3, n3, Cj)) 

1082 # Here are the full expressions for m and s 

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

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

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

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

1087 cd *= b1[2] 

1088 sd *= b1[3] * _2_0 

1089 s = _fsum1f_(cd * sp, sd * cp) if sinp else \ 

1090 _fsum1f_(cd * cp, _neg(sd * sp), _neg(b2[2])) 

1091 return s 

1092 

1093 

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

1095 4: ( # GEOGRAPHICLIB_RHUMBAREA_ORDER == 4 

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

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

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

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

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

1101 5: ( # GEOGRAPHICLIB_RHUMBAREA_ORDER == 5 

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

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

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

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

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

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

1108 6: ( # GEOGRAPHICLIB_RHUMBAREA_ORDER == 6 

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

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

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

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

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

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

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

1116 7: ( # GEOGRAPHICLIB_RHUMBAREA_ORDER == 7 

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

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

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

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

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

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

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

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

1125 8: ( # GEOGRAPHICLIB_RHUMBAREA_ORDER == 8 

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

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

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

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

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

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

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

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

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

1135} 

1136 

1137__all__ += _ALL_DOCS(Caps, _RhumbLine) 

1138 

1139if __name__ == '__main__': 

1140 

1141 from pygeodesy.lazily import printf 

1142 

1143 def _re(fmt, r3, x3): 

1144 e3 = [] 

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

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

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

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

1149 

1150 # <https://GeographicLib.SourceForge.io/cgi-bin/RhumbSolve> 

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

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

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

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

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

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

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

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

1159 

1160# % python3 -m pygeodesy.rhumbx 

1161 

1162# 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, 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) 

1163 

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

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

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

1167 

1168# **) MIT License 

1169# 

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

1171# 

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

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

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

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

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

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

1178# 

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

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

1181# 

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

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

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

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

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

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

1188# OTHER DEALINGS IN THE SOFTWARE.