1

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

3

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

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

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

10

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

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

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

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

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

16

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

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

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

20'''

21# make sure int/int division yields float quotient

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

23

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

25 from 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

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

29 _xdatum, _xkwds

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

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

32# from pygeodesy.fsums import fsum1_ # from .fmath

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

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

35 _no_, _s12_, _S12_, _UNDER

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

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

38 _GTuple, Inverse10Tuple, _norm180

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

40 _AlpCoeffs, _BetCoeffs # PYCHOK used!

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

42 from pygeodesy.namedTuples import Distance2Tuple, LatLon2Tuple, NearestOn4Tuple

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

44 _update_all

45 from pygeodesy.streprs import Fmt, pairs, unstr

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

47 Meter as _M, Meter2 as _M2

48 from pygeodesy.utily import sincos2_, sincos2d

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

50

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

52

53__all__ = _ALL_LAZY.rhumbx

54__version__ = '23.04.10'

55

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

57_TRIPS = 65 # .intersection2, 18+

58

59

60 class _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

68 class _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

76 def _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

88 class 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, lat1, lon1, azi12, s12, outmask=Caps.LATITUDE_LONGITUDE):

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

143

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

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

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

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

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

149

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

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

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

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

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

155

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

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

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

159

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

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

162 This position is extremely close to the actual pole and

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

164 '''

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

166 name=self.name)

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

168

169 @deprecated_method

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

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

172

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

174 '''

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

176

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

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

179 '''

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

181

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

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

184

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

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

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

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

189

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

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

192

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

194 rhumb line.

195 '''

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

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

198

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

200 E = self.ellipsoid

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

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

203

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

205 if self.exact:

206 E = self.ellipsoid

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

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

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

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

211 else:

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

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

214 return r

215

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

217 E = self.ellipsoid

218 tbetx = E.f1 * tphix

219 tbety = E.f1 * tphiy

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

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

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

223

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

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

226

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

228 E = self.ellipsoid

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

230 if self.exact:

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

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

233 self._DRectifyingt( tphix, tphiy, t)

234 else:

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

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

237 return r

238

239 @Property_RO

240 def _eF(self):

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

242 '''

243 # .k2 = 0.006739496742276434

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

245

246 @Property

247 def ellipsoid(self):

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

249 '''

250 return self._E

251

252 @ellipsoid.setter # PYCHOK setter!

253 def ellipsoid(self, a_earth_f):

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

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

256 '''

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

258 if self._E != E:

259 _update_all_rls(self)

260 self._E = E

261

262 @property_RO

263 def equatoradius(self):

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

265 '''

266 return self.ellipsoid.a

267

268 a = equatoradius

269

270 @Property

271 def exact(self):

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

273 '''

274 return self._exact

275

276 @exact.setter # PYCHOK setter!

277 def exact(self, exact):

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

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

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

281

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

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

284 '''

285 x = bool(exact)

286 if self._exact != x:

287 _update_all_rls(self)

288 self._exact = x

289

290 def flattening(self):

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

292 '''

293 return self.ellipsoid.f

294

295 f = flattening

296

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

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

299

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

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

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

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

304

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

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

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

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

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

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

311

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

313 on opposite meridians, there are two shortest rhumb lines

314 and the East-going one is chosen.

315

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

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

318 This position is extremely close to the actual pole and

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

320 '''

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

322 if (outmask & Cs.AZIMUTH_DISTANCE_AREA):

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

324 E = self.ellipsoid

325 psi1 = E.auxIsometric(lat1)

326 psi2 = E.auxIsometric(lat2)

327 psi12 = psi2 - psi1

328 lon12, _ = _diff182(lon1, lon2)

329 if (outmask & Cs.AZIMUTH):

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

331 if (outmask & Cs.DISTANCE):

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

333 s12 = a12 * E._L_90

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

335 if (outmask & Cs.AREA):

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

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

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

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

340 lon12=lon12, psi1=psi1, exact=self.exact,

341 psi12=psi12, psi2=psi2)

342 return r

343

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

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

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

347#

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

349# '''

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

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

352

353 @deprecated_method

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

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

356

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

358 '''

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

360

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

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

363 '''

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

365

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

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

368

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

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

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

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

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

374

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

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

377

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

379 rhumb line.

380 '''

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

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

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

384

385 Line = DirectLine # synonyms

386

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

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

389 if self.f:

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

391 return s

392

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

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

395

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

397 or 8).

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

399 5, 6, 7 or 8).

400

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

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

403 '''

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

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

406 self.RAorder = RAorder

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

408 self.TMorder = TMorder

409 return t

410

411 @Property_RO

412 def _RA2(self):

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

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

415 m = self.RAorder

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

417

418 @Property

419 def RAorder(self):

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

421 '''

422 return self._mRA

423

424 @RAorder.setter # PYCHOK setter!

425 def RAorder(self, order):

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

427 '''

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

429 if self._mRA != n:

430 _update_all_rls(self)

431 self._mRA = n

432

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

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

435 '''

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

437 self.ellipsoid.area) * lon12 / _720_0

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

439 return r

440

441 @Property

442 def TMorder(self):

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

444 '''

445 return self._mTM

446

447 @TMorder.setter # PYCHOK setter!

448 def TMorder(self, order):

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

450

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

452 '''

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

454 if self._mTM != n:

455 _update_all_rls(self)

456 self._mTM = n

457 self.exact = False

458

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

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

461

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

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

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

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

466

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

468 '''

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

470 exact=self.exact, TMorder=self.TMorder)

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

472

473

474 class RhumbError(_ValueError):

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

476 '''

477 pass

478

479

480 class _RhumbLine(_RhumbBase):

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

482 '''

483 _azi12 = _0_0

484# _lat1 = _0_0

485# _lon1 = _0_0

486 _salp = _0_0

487 _calp = _1_0

488 _rhumb = None # L{Rhumb} instance

489

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

491 '''New C{RhumbLine}.

492 '''

493 _xinstanceof(Rhumb, rhumb=rhumb)

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

495 self._lon1 = _Lon(lon1= lon1)

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

497 if azi12: # non-zero

498 self.azi12 = azi12

499 self._caps = caps

500 if not (caps & Caps.LINE_OFF):

501 _rls.append(self)

502 n = name or rhumb.name

503 if n:

504 self.name=n

505 self._rhumb = rhumb # last

506

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

508 if _rls: # may be empty or None

509 try: # PYCHOK no cover

510 _rls.remove(self)

511 except (TypeError, ValueError):

512 pass

513 self._rhumb = None

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

515

516 @Property

517 def azi12(self):

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

519 '''

520 return self._azi12

521

522 @azi12.setter # PYCHOK setter!

523 def azi12(self, azi12):

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

525 '''

526 z = _norm180(azi12)

527 if self._azi12 != z:

528 if self._rhumb:

529 _update_all(self)

530 self._azi12 = z

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

532

533 def distance2(self, lat, lon):

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

535 point to this rhumb line's start point.

536

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

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

539

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

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

542

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

544 '''

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

546# outmask=Caps.AZIMUTH_DISTANCE)

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

548

549 @Property_RO

550 def ellipsoid(self):

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

552 '''

553 return self.rhumb.ellipsoid

554

555 @property_RO

556 def exact(self):

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

558 '''

559 return self.rhumb.exact

560

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

562 '''Iteratively find the intersection of this and an other rhumb line.

563

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

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

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

567

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

569 C{lon}gitude of the intersection point.

570

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

572 no intersection for an other reason.

573

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

575 and function L{pygeodesy.intersection3d3}.

576

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

578 L{ExactTransverseMercator} or L{KTransverseMercator}

579 projection and invoking function L{intersection3d3} in

580 that domain.

581 '''

582 _xinstanceof(other, _RhumbLine)

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

584 try:

585 if other is self:

586 raise ValueError(_coincident_)

587 # make globals and invariants locals

588 _diff = euclid # approximate length

589 _i3d3 = _intersect3d3 # NOT .vector3d.intersection3d3

590 _LL2T = LatLon2Tuple

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

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

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

594 # use halfway point as initial estimate

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

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

597 for i in range(1, _TRIPS):

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

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

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

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

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

603 if d < tol:

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

605 name=self.intersection2.__name__)

606 except Exception as x:

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

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

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

610

611 @property_RO

612 def lat1(self):

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

614 '''

615 return self._lat1

616

617 @property_RO

618 def lon1(self):

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

620 '''

621 return self._lon1

622

623 @Property_RO

624 def latlon1(self):

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

626 '''

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

628

629 @Property_RO

630 def _mu1(self):

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

632 '''

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

634

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

636 '''Iteratively locate the point on this rhumb line nearest to

637 the given point.

638

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

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

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

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

643

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

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

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

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

648

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

650 no intersection for an other reason.

651

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

653 and function L{intersection3d3}.

654 '''

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

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

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

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

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

660 iteration=p.iteration)

661

662 @Property_RO

663 def _psi1(self):

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

665 '''

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

667

668 @property_RO

669 def RAorder(self):

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

671 '''

672 return self.rhumb.RAorder

673

674 @Property_RO

675 def _r1rad(self): # PYCHOK no cover

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

677 '''

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

679

680 @Property_RO

681 def rhumb(self):

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

683 '''

684 return self._rhumb

685

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

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

688

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

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

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

692 the quantities to be returned.

693

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

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

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

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

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

699

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

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

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

703

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

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

706 This position is extremely close to the actual pole and

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

708 '''

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

710 if (outmask & Cs.LATITUDE_LONGITUDE_AREA):

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

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

713 mu2 = mu12 + self._mu1

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

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

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

717 mu2 = _norm180(_180_0 - mu2)

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

719 lon2x = NAN

720 if (outmask & Cs.AREA):

721 r.set_(S12=NAN)

722 else:

723 psi2 = self._psi1

724 if self._calp:

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

726 psi12 = R._DRectifying2Isometricd(mu2,

727 self._mu1) * mu12

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

729 psi2 += psi12

730 else: # PYCHOK no cover

731 lat2x = self.lat1

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

733 if (outmask & Cs.AREA):

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

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

736 if (outmask & Cs.LATITUDE):

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

738 if (outmask & Cs.LONGITUDE):

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

740 lon2x += self.lon1

741 else:

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

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

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

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

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

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

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

749 return r

750

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

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

753

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

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

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

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

758

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

760 '''

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

762 azi12=self.azi12, exact=self.exact,

763 TMorder=self.TMorder, xTM=self.xTM)

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

765

766 @property_RO

767 def TMorder(self):

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

769 '''

770 return self.rhumb.TMorder

771

772 @Property_RO

773 def xTM(self):

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

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

776 '''

777 E = self.ellipsoid

778 # ExactTransverseMercator doesn't handle spherical earth models

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

780 KTransverseMercator(E, TMorder=self.TMorder)

781

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

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

784 as current intersection estimate and central meridian.

785 '''

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

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

788

789

790 class RhumbLine(_RhumbLine):

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

792

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

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

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

796

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

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

799 area C{S12} under the rhumb line.

800

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

802 two rhumb lines.

803

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

805 its distance to a rhumb line.

806 '''

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

808 '''New L{RhumbLine}.

809

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

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

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

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

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

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

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

817 in this rhumb line.

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

819 '''

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

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

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

823

824

825 class RhumbOrder2Tuple(_GTuple):

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

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

828 '''

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

830 _Units_ = ( Int, Int)

831

832

833 class Rhumb8Tuple(_GTuple):

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

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

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

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

838 '''

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

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

841

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

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

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

845

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

847 @kwarg a12_azi1_azi2_m12_M12_M21: Optional keyword arguments

848 to specify or override L{Inverse10Tuple} items.

849

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

851 m12, M12, M21, S12)}

852 '''

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

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

855 if a12_azi1_azi2_m12_M12_M21:

856 d.update(a12_azi1_azi2_m12_M12_M21)

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

858

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

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

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

862

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

864 @kwarg a12_m12_M12_M21_salp1_calp1_salp2_calp2: Optional keyword

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

866

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

868 m12, M12, M21, S12)}.

869 '''

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

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

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

873 if a12_m12_M12_M21_salp1_calp1_salp2_calp2:

874 d.update(a12_m12_M12_M21_salp1_calp1_salp2_calp2)

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

876

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

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

879 '''

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

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

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

883

884 @deprecated_method

885 def _to7Tuple(self):

886 '''DEPRECATED, do not use!

887 '''

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

889

890

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

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

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

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

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

896

897 def _Dasinh(x, y):

898 hx = hypot1(x)

899 d = x - y

900 if d:

901 hx *= y

902 hy = x * hypot1(y)

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

904 r = asinh(t) / d

905 else:

906 r = _1_0 / hx

907 return r

908

909

910 def _Datan(x, y):

911 xy = x * y

912 r = xy + _1_0

913 d = x - y

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

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

916 else:

917 r = _1_0 / r

918 return r

919

920

921 def _Dcosh(x, y):

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

923

924

925 def _DeatanhE(x, y, E):

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

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

928 # assert not isnear0(e)

929 d = x - y

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

931

932

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

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

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

936 d = x - y

937 if (x * y) > 0:

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

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

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

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

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

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

944 # = t = d * Dt

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

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

947 # holds if x*y > 0):

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

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

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

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

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

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

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

955 t = D * d

956 t2 = t**2 + _1_0

957 D *= _2_0 / t2

958 s = D * d

959 if s:

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

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

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

963 elif d:

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

965 return r

966

967

968 def _Dgd(x, y):

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

970

971

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

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

974

975

976 def _Dlog(x, y):

977 d = (x - y) * _0_5

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

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

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

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

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

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

984

985

986 def _Dsin(x, y):

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

988

989

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

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

992 d = (x - y) * _0_5

993 if d:

994 r *= sin_(d) / d

995 return r

996

997

998 def _Dsinh(x, y):

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

1000

1001

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

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

1004

1005

1006 def _Dtant(dxy, tx, ty):

1007 txy = tx * ty

1008 r = txy + _1_0

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

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

1011 return r

1012

1013

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

1015 # get inverse auxiliary lats in radians and tangents

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

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

1018 return phix, phiy, tan(phix), tan(phiy)

1019

1020

1021 def _gd(x):

1022 return atan(sinh(x))

1023

1024

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

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

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

1028 # Use Clenshaw summation to evaluate

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

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

1031 # where

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

1033 # SC = sinp ? sin : cos

1034 # CS = sinp ? cos : sin

1035 # ...

1036 d = x - y

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

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

1039 m = cp * cd * _2_0

1040 s = neg(sp * sd) # negative

1041 # 2x2 matrices in row-major order

1042 a0, a1 = m, (s * d**2)

1043 a2, a3 = (s * _4_0), m

1044 b2 = b1 = _0_0s(4)

1045 if n > 0:

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

1047 _fsum1_, _neg = fsum1_, neg

1048 for j in range(n - 1, 0, -1):

1049 b1, b2, Cj = b2, b1, C[j] # C[0] unused

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

1051 m0, m1, m2, m3 = b2

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

1053 b1 = (_fsum1_(a0 * m0, a1 * m2, n0, Cj, floats=True),

1054 _fsum1_(a0 * m1, a1 * m3, n1, floats=True),

1055 _fsum1_(a2 * m0, a3 * m2, n2, floats=True),

1056 _fsum1_(a2 * m1, a3 * m3, n3, Cj, floats=True))

1057 # Here are the full expressions for m and s

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

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

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

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

1062 cd *= b1[2]

1063 sd *= b1[3] * _2_0

1064 s = _fsum1_(cd * sp, sd * cp, floats=True) if sinp else \

1065 _fsum1_(cd * cp, _neg(sd * sp), _neg(b2[2]), floats=True)

1066 return s

1067

1068

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

1070 4: ( # GEOGRAPHICLIB_RHUMBAREA_ORDER == 4

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

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

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

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

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

1076 5: ( # GEOGRAPHICLIB_RHUMBAREA_ORDER == 5

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

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

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

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

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

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

1083 6: ( # GEOGRAPHICLIB_RHUMBAREA_ORDER == 6

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

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

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

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

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

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

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

1091 7: ( # GEOGRAPHICLIB_RHUMBAREA_ORDER == 7

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

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

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

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

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

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

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

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

1100 8: ( # GEOGRAPHICLIB_RHUMBAREA_ORDER == 8

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

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

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

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

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

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

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

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

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

1110}

1111

1112__all__ += _ALL_DOCS(Caps, _RhumbLine)

1113

1114 if __name__ == '__main__':

1115

1116 def _re(fmt, r3, x3):

1117 e3 = []

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

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

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

1121 print((fmt % r3) + ' rel errors: ' + ', '.join(e3))

1122

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

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

1125 print('# %r\n' % rhumb)

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

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

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

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

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

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

1132

1133# % python3 -m pygeodesy.rhumbx

1134

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

1136

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

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

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

1140

1141# **) MIT License

1142#

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

1144#

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

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

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

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

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

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

1151#

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

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

1154#

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

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

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

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

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

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

1161# OTHER DEALINGS IN THE SOFTWARE.