1

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

3

4 u'''(INTERNAL) base classes C{RhumbBase} and C{RhumbLineBase}, pure Python version of I{Karney}'s

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

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

7 from I{GeographicLib versions 2.0} and I{2.2} and I{Karney}'s C++ example U{Rhumb intersect

8<https://SourceForge.net/p/geographiclib/discussion/1026620/thread/2ddc295e/>}.

9

10 Class L{RhumbLineBase} has been enhanced with methods C{Intersecant2}, C{Intersection} and C{PlumbTo}

11 to iteratively find the intersection of a rhumb line and a circle or an other rhumb line, respectively

12 a perpendicular geodesic or other rhumb line.

13

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

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

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

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

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

19

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

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

22'''

23# make sure int/int division yields float quotient

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

25

26 from pygeodesy.basics import _copysign, itemsorted, unsigned0, _xinstanceof

27 from pygeodesy.constants import EPS, EPS0, EPS1, INT0, NAN, _over, \

28 _EPSqrt as _TOL, _0_0, _0_01, _1_0, _90_0

29 from pygeodesy.datums import Datum, _earth_datum, _spherical_datum, _WGS84

30 from pygeodesy.errors import IntersectionError, RhumbError, _xdatum, \

31 _xkwds, _xkwds_pop2, _Xorder

32# from pygeodesy.etm import ExactTransverseMercator # _MODS

33 from pygeodesy.fmath import euclid, favg, sqrt_a, Fsum

34# from pygeodesy.formy import opposing # _MODS

35# from pygeodesy.fsums import Fsum # from .fmath

36 from pygeodesy.interns import NN, _coincident_, _COMMASPACE_, _Dash, \

37 _dunder_nameof, _parallel_, _too_, _under

38 from pygeodesy.karney import _atan2d, Caps, _CapsBase, _diff182, _fix90, \

39 _norm180, GDict

40# from pygeodesy.ktm import KTransverseMercator, _AlpCoeffs # _MODS

41 from pygeodesy.lazily import _ALL_DOCS, _ALL_MODS as _MODS

42 from pygeodesy.namedTuples import Distance2Tuple, LatLon2Tuple

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

44 property_RO, _update_all

45 from pygeodesy.streprs import Fmt, pairs

46 from pygeodesy.units import Float_, Lat, Lon, Meter, Radius_, Int # PYCHOK shared

47 from pygeodesy.utily import acos1, _azireversed, _loneg, sincos2d, sincos2d_, \

48 _unrollon, _Wrap

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

50

51 from math import cos, fabs

52

53__all__ = ()

54__version__ = '24.04.07'

55

56_anti_ = _Dash('anti')

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

58_TRIPS = 65 # .Intersection, .PlumbTo, 19+

59

60

61 class _Lat(Lat):

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

63 '''

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

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

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

67

68

69 class _Lon(Lon):

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

71 '''

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

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

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

75

76

77 def _update_all_rls(r):

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

79 of any C{RhumbLine} instances tied to the given

80 C{Rhumb} instance B{C{r}}.

81 '''

82 # _xinstanceof(_MODS.rhumb.aux_.RhumbAux, _MODS.rhumb.ekx.Rhumb, r=r)

83 _update_all(r)

84 for rl in _rls: # PYCHOK use weakref?

85 if rl._rhumb is r:

86 _update_all(rl)

87

88

89 class RhumbBase(_CapsBase):

90 '''(INTERNAL) Base class for C{rhumb.aux_.RhumbAux} and C{rhumb.ekx.Rhumb}.

91 '''

92 _datum = _WGS84

93 _exact = True

94 _f_max = _0_01

95 _mTM = 6 # see .TMorder

96

97 def __init__(self, a_earth, f, exact, name):

98 '''New C{RhumbAux} or C{Rhumb}.

99 '''

100 _earth_datum(self, a_earth, f=f, name=name)

101 if not exact:

102 self.exact = False

103 if name:

104 self.name = name

105

106 @Property_RO

107 def a(self):

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

109 '''

110 return self.ellipsoid.a

111

112 equatoradius = a

113

114 def ArcDirect(self, lat1, lon1, azi12, a12, outmask=Caps.LATITUDE_LONGITUDE):

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

116

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

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

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

120 @arg a12: Angle along the rhumb line from the given to the

121 destination point (C{degrees}), can be negative.

122

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

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

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

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

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

128

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

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

131 is True} and L{RhumbAux}.

132

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

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

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

136

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

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

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

140 carried out in finite terms.

141 '''

142 s12 = a12 * self._mpd

143 return self._DirectRhumb(lat1, lon1, azi12, a12, s12, outmask)

144

145 @Property_RO

146 def b(self):

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

148 '''

149 return self.ellipsoid.b

150

151 polaradius = b

152

153 @property

154 def datum(self):

155 '''Get this rhumb's datum (L{Datum}).

156 '''

157 return self._datum

158

159 @datum.setter # PYCHOK setter!

160 def datum(self, datum):

161 '''Set this rhumb's datum (L{Datum}).

162

163 @raise RhumbError: If C{abs(B{f}} exceeds non-zero C{f_max} and C{exact=False}.

164 '''

165 _xinstanceof(Datum, datum=datum)

166 if self._datum != datum:

167 self._exactest(self.exact, datum.ellipsoid, self.f_max)

168 _update_all_rls(self)

169 self._datum = datum

170

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

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

173 '''

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

175

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

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

178

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

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

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

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

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

184

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

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

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

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

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

190

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

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

193 is True} and L{RhumbAux}.

194

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

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

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

198

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

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

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

202 carried out in finite terms.

203 '''

204 a12 = _over(s12, self._mpd)

205 return self._DirectRhumb(lat1, lon1, azi12, a12, s12, outmask)

206

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

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

209 '''

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

211

212 def _DirectLine(self, ll1, azi12, **caps_name):

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

214 '''

215 return self.DirectLine(ll1.lat, ll1.lon, azi12, **caps_name)

216

217 def DirectLine(self, lat1, lon1, azi12, **caps_name):

218 '''Define a C{RhumbLine} in terms of the I{direct} rhumb

219 problem to compute several points on a single rhumb line.

220

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

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

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

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

225 C{B{caps}=Caps.STANDARD}, a bit-or'ed combination of

226 L{Caps} values specifying the required capabilities.

227 Include C{Caps.LINE_OFF} if updates to the B{C{rhumb}}

228 should I{not} be reflected in this rhumb line.

229

230 @return: A C{RhumbLine...} instance and invoke its method

231 C{.Position} to compute each point.

232

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

234 rhumb line, unless C{B{caps} |= Caps.LINE_OFF}.

235 '''

236 return self._RhumbLine(self, lat1, lon1, azi12, **caps_name)

237

238 Line = DirectLine # synonyms

239

240 def _DirectRhumb(self, lat1, lon1, azi12, a12, s12, outmask):

241 '''(INTERNAL) See methods C{.ArcDirect} and C{.Direct}.

242 '''

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

244 name=self.name)

245 return rl._Position(a12, s12, outmask | self._debug) # lat2, lon2, S12

246

247 @Property

248 def ellipsoid(self):

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

250 '''

251 return self.datum.ellipsoid

252

253 @ellipsoid.setter # PYCHOK setter!

254 def ellipsoid(self, a_earth_f):

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

256 L{a_f2Tuple}) or (equatorial) radius and flattening (2-tuple C{(a, f)}).

257

258 @raise RhumbError: If C{abs(B{f}} exceeds non-zero C{f_max} and C{exact=False}.

259 '''

260 self.datum = _spherical_datum(a_earth_f, Error=RhumbError)

261

262 @Property

263 def exact(self):

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

265 '''

266 return self._exact

267

268 @exact.setter # PYCHOK setter!

269 def exact(self, exact):

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

271 expressions, otherwise a series expansion (accurate for oblate or

272 prolate ellipsoids with C{abs(flattening)} below C{f_max}.

273

274 @raise RhumbError: If C{B{exact}=False} and C{abs(flattening})

275 exceeds non-zero C{f_max}.

276

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

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

279 '''

280 x = bool(exact)

281 if self._exact != x:

282 self._exactest(x, self.ellipsoid, self.f_max)

283 _update_all_rls(self)

284 self._exact = x

285

286 def _exactest(self, exact, ellipsoid, f_max):

287 # Helper for property setters C{ellipsoid}, C{exact} and C{f_max}

288 if fabs(ellipsoid.f) > f_max > 0 and not exact:

289 raise RhumbError(exact=exact, f=ellipsoid.f, f_max=f_max)

290

291 @Property_RO

292 def f(self):

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

294 '''

295 return self.ellipsoid.f

296

297 flattening = f

298

299 @property

300 def f_max(self):

301 '''Get the I{max.} flattening (C{float}).

302 '''

303 return self._f_max

304

305 @f_max.setter # PYCHOK setter!

306 def f_max(self, f_max): # PYCHOK no cover

307 '''Set the I{max.} flattening, not to exceed (C{float}).

308

309 @raise RhumbError: If C{exact=False} and C{abs(flattening})

310 exceeds non-zero C{f_max}.

311 '''

312 f = Float_(f_max=f_max, low=_0_0, high=EPS1)

313 if self._f_max != f:

314 self._exactest(self.exact, self.ellipsoid, f)

315 self._f_max = f

316

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

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

319 '''

320 if wrap:

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

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

323

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

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

326

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

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

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

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

331

332 @return: L{GDict} with 4 to 9 items C{lat1, lon1, lat2, lon2,

333 azi12, azi21, s12, a12, S12}, the rhumb line's azimuth

334 C{azi12} and I{reverse} azimuth C{azi21}, both in

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

336 rhumb distance C{s12} and rhumb angle C{a12} between

337 both points in C{meter} respectively C{degrees} and

338 the area C{S12} under the rhumb line in C{meter}

339 I{squared}.

340

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

342 only required for L{RhumbAux} area C{S12}

343 when C{B{exact} is True}.

344

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

346 on opposite meridians, there are two shortest rhumb lines

347 and the East-going one is chosen.

348

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

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

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

352 carried out in finite terms.

353 '''

354 r = GDict(lat1=lat1, lon1=lon1, lat2=lat2, lon2=lon2, name=self.name)

355 Cs = Caps

356 if (outmask & Cs.AZIMUTH_DISTANCE_AREA):

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

358 y, x, s1, s2 = self._Inverse4(lon12, r, outmask)

359 if (outmask & Cs.AZIMUTH):

360 z = _atan2d(y, x)

361 r.set_(azi12=z, azi21=_azireversed(z))

362 if (outmask & Cs.AREA):

363 S12 = self._S12d(s1, s2, lon12)

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

365 return r

366

367 def _Inverse4(self, lon12, r, outmask): # PYCHOK no cover

368 '''(INTERNAL) I{Must be overloaded}.'''

369 self._notOverloaded(lon12, r, Caps.toStr(outmask)) # underOK=True

370

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

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

373 '''

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

375

376 def _InverseLine(self, ll1, ll2, wrap, **caps_name):

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

378 '''

379 if wrap:

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

381 return self.InverseLine(ll1.lat, ll1.lon, ll2.lat, ll2.lon, **caps_name)

382

383 def InverseLine(self, lat1, lon1, lat2, lon2, **caps_name):

384 '''Define a C{RhumbLine} in terms of the I{inverse} rhumb problem.

385

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

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

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

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

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

391 C{B{caps}=Caps.STANDARD}, a bit-or'ed combination of

392 L{Caps} values specifying the required capabilities.

393 Include C{Caps.LINE_OFF} if updates to the B{C{rhumb}}

394 should I{not} be reflected in this rhumb line.

395

396 @return: A C{RhumbLine...} instance and invoke its method

397 C{ArcPosition} or C{Position} to compute points.

398

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

400 rhumb line, unless C{B{caps} |= Caps.LINE_OFF}.

401 '''

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

403 return self._RhumbLine(self, lat1, lon1, r.azi12, **caps_name)

404

405 @Property_RO

406 def _mpd(self): # PYCHOK no cover

407 '''(INTERNAL) I{Must be overloaded}.'''

408 _MODS.named.notOverloaded(self)

409

410 @property_RO

411 def RAorder(self):

412 '''Get the I{Rhumb Area} order, C{None} always.

413 '''

414 return None

415

416 @property_RO

417 def _RhumbLine(self): # PYCHOK no cover

418 '''(INTERNAL) I{Must be overloaded}.'''

419 self._notOverloaded(underOK=True)

420

421 def _S12d(self, s1, s2, lon): # PYCHOK no cover

422 '''(INTERNAL) I{Must be overloaded}.'''

423 self._notOverloaded(s1, s2, lon) # underOK=True

424

425 @Property

426 def TMorder(self):

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

428 '''

429 return self._mTM

430

431 @TMorder.setter # PYCHOK setter!

432 def TMorder(self, order):

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

434

435 @note: Setting C{TMorder} turns property C{exact} off, but only

436 for L{Rhumb} instances.

437 '''

438 m = _Xorder(_MODS.ktm._AlpCoeffs, RhumbError, TMorder=order)

439 if self._mTM != m:

440 _update_all_rls(self)

441 self._mTM = m

442 if self.exact and isinstance(self, _MODS.rhumb.ekx.Rhumb):

443 self.exact = False

444

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

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

447

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

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

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

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

452

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

454 '''

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

456 exact=self.exact, TMorder=self.TMorder)

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

458

459

460 class RhumbLineBase(_CapsBase):

461 '''(INTERNAL) Base class for C{rhumb.aux_.RhumbLineAux} and C{rhumb.ekx.RhumbLine}.

462 '''

463 _azi12 = _0_0

464 _calp = _1_0

465# _caps = \

466# _debug = 0

467# _lat1 = \

468# _lon1 = \

469# _lon12 = _0_0

470 _Rhumb = RhumbBase # compatible C{Rhumb} class

471 _rhumb = None # C{Rhumb} instance

472 _salp = \

473 _talp = _0_0

474

475 def __init__(self, rhumb, lat1, lon1, azi12, caps=Caps.STANDARD, name=NN):

476 '''New C{RhumbLine} or C{RhumbLineAux}.

477 '''

478 _xinstanceof(self._Rhumb, rhumb=rhumb)

479

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

481 self._lon1 = _Lon(lon1= lon1)

482 self._lon12 = _norm180(self._lon1)

483 if azi12: # non-zero, non-None

484 self.azi12 = _norm180(azi12)

485

486 n = name or rhumb.name

487 if n:

488 self.name=n

489

490 self._caps = caps

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

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

493 self._rhumb = rhumb.copy(deep=False, name=_under(rhumb.name))

494 else:

495 self._rhumb = rhumb

496 _rls.append(self)

497

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

499 if _rls: # may be empty or None

500 try: # PYCHOK no cover

501 _rls.remove(self)

502 except (TypeError, ValueError):

503 pass

504 self._rhumb = None

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

506

507 def ArcPosition(self, a12, outmask=Caps.LATITUDE_LONGITUDE):

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

509

510 @arg a12: The angle along this rhumb line from its origin to the

511 point (C{degrees}), can be negative.

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

513 the quantities to be returned.

514

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

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

517 C{lon2} of the point in C{degrees}, the rhumb distance

518 C{s12} in C{meter} from the start point of and the area

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

520

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

522 only required for L{RhumbLineAux} area C{S12}

523 when C{B{exact} is True}.

524

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

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

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

528

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

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

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

532 carried out in finite terms.

533 '''

534 return self._Position(a12, self.degrees2m(a12), outmask)

535

536 @Property

537 def azi12(self):

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

539 '''

540 return self._azi12

541

542 @azi12.setter # PYCHOK setter!

543 def azi12(self, azi12):

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

545 '''

546 z = _norm180(azi12)

547 if self._azi12 != z:

548 if self._rhumb:

549 _update_all(self)

550 self._azi12 = z

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

552 self._talp = _over(*t)

553

554 @property_RO

555 def azi12_sincos2(self): # PYCHOK no cover

556 '''Get the sine and cosine of this rhumb line's I{azimuth} (2-tuple C{(sin, cos)}).

557 '''

558 return self._scalp, self._calp

559

560 @property_RO

561 def datum(self):

562 '''Get this rhumb line's datum (L{Datum}).

563 '''

564 return self.rhumb.datum

565

566 def degrees2m(self, angle):

567 '''Convert an angular distance along this rhumb line to C{meter}.

568

569 @arg angle: Angular distance (C{degrees}).

570

571 @return: Distance (C{meter}).

572 '''

573 return float(angle) * self.rhumb._mpd

574

575 @deprecated_method

576 def distance2(self, lat, lon): # PYCHOK no cover

577 '''DEPRECATED on 23.09.23, use method L{RhumbLineAux.Inverse} or L{RhumbLine.Inverse}.

578

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

580 in C{meter} and C{initial} bearing (azimuth) in C{degrees}.

581 '''

582 r = self.Inverse(lat, lon)

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

584

585 @property_RO

586 def ellipsoid(self):

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

588 '''

589 return self.rhumb.ellipsoid

590

591 @property_RO

592 def exact(self):

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

594 '''

595 return self.rhumb.exact

596

597 def Intersecant2(self, lat0, lon0, radius, napier=True, **tol_eps):

598 '''Compute the intersection(s) of this rhumb line and a circle.

599

600 @arg lat0: Latitude of the circle center (C{degrees}).

601 @arg lon0: Longitude of the circle center (C{degrees}).

602 @arg radius: Radius of the circle (C{meter}, conventionally).

603 @kwarg napier: If C{True}, apply I{Napier}'s spherical triangle

604 instead of planar trigonometry (C{bool}).

605 @kwarg tol_eps: Optional keyword arguments, see method

606 method L{Intersection} for further details.

607

608 @return: 2-Tuple C{(P, Q)} with both intersections (representing

609 a rhumb chord), each a L{GDict} from method L{Intersection}

610 extended to 18 items by C{lat3, lon3, azi03, a03, s03}

611 with azimuth C{azi03} of, distance C{a03} in C{degrees}

612 and C{s03} in C{meter} along the rhumb line from the circle

613 C{lat0, lon0} to the chord center C{lat3, lon3}. If this

614 rhumb line is tangential to the circle, both points

615 are the same L{GDict} instance with distances C{s02} and

616 C{s03} near-equal to the B{C{radius}}.

617

618 @raise IntersectionError: The circle and this rhumb line

619 do not intersect.

620

621 @raise UnitError: Invalid B{C{radius}}.

622 '''

623 r = Radius_(radius)

624 p = q = self.PlumbTo(lat0, lon0, exact=None, **tol_eps)

625 a = q.s02

626 t = dict(lat3=q.lat2, lon3=q.lon2, azi03=q.azi02, a03=q.a02, s03=a)

627 if a < r:

628 t.update(iteration=q.iteration, lat0=q.lat1, lon0=q.lon1, # or lat0, lon0

629 name=_dunder_nameof(self.Intersecant2, self.name))

630 if fabs(a) < EPS0: # coincident centers

631 d, h = _0_0, r

632 else:

633 d = q.s12

634 if napier: # Napier rule (R1) cos(b) = cos(c) / cos(a)

635 # <https://WikiPedia.org/wiki/Spherical_trigonometry>

636 m = self.rhumb._mpr

637 h = (acos1(cos(r / m) / cos(a / m)) * m) if m else _0_0

638 else:

639 h = _copysign(sqrt_a(r, a), a)

640 p = q = self.Position(d + h).set_(**t)

641 if h:

642 q = self.Position(d - h).set_(**t)

643 elif a > r:

644 t = _too_(Fmt.distant(a))

645 raise IntersectionError(self, lat0, lon0, radius,

646 txt=t, **tol_eps)

647 else: # tangential

648 q.set_(**t) # == p.set(_**t)

649 return p, q

650

651 @deprecated_method

652 def intersection2(self, other, **tol_eps): # PYCHOK no cover

653 '''DEPRECATED on 23.10.10, use method L{Intersection}.'''

654 p = self.Intersection(other, **tol_eps)

655 r = LatLon2Tuple(p.lat2, p.lon2, name=self.intersection2.__name__)

656 r._iteration = p.iteration

657 return r

658

659 def Intersection(self, other, tol=_TOL, **eps):

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

661

662 @arg other: The other rhumb line (C{RhumbLine}).

663 @kwarg tol: Tolerance for longitudinal convergence and parallel

664 error (C{degrees}).

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

666

667 @return: The intersection point, a L{Position}-like L{GDict} with

668 13 items C{lat1, lon1, azi12, a12, s12, lat2, lon2, lat0,

669 lon0, azi02, a02, s02, at} with the rhumb angle C{a02}

670 and rhumb distance C{s02} between the start point C{lat0,

671 lon0} of the B{C{other}} rhumb line and the intersection

672 C{lat2, lon2}, the azimuth C{azi02} of the B{C{other}}

673 rhumb line and the angle C{at} between both rhumb lines.

674 See method L{Position} for further details.

675

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

677 no intersection for an other reason.

678

679 @see: Methods C{distance2} and C{PlumbTo} and function

680 L{pygeodesy.intersection3d3}.

681

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

683 L{ExactTransverseMercator} or L{KTransverseMercator}

684 projection and function L{pygeodesy.intersection3d3} in

685 that domain.

686 '''

687 _xinstanceof(RhumbLineBase, other=other)

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

689 try:

690 if self.others(other) is self:

691 raise ValueError(_coincident_)

692 # make invariants and globals locals

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

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

695 p = _MODS.formy.opposing(s_az, o_az, margin=tol)

696 if p is not None: # == p in (True, False)

697 raise ValueError(_anti_(_parallel_) if p else _parallel_)

698 _diff = euclid # approximate length

699 _i3d3 = _intersect3d3 # NOT .vector3d.intersection3d3

700 _LL2T = LatLon2Tuple

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

702 # use halfway point as initial estimate

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

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

705 for i in range(1, _TRIPS):

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

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

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

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

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

711 if d < tol: # 19 trips

712 break

713 else:

714 raise ValueError(Fmt.no_convergence(d, tol))

715

716 P = GDict(lat1=self.lat1, lat2=p.lat, lat0=other.lat1,

717 lon1=self.lon1, lon2=p.lon, lon0=other.lon1,

718 name=_dunder_nameof(self.Intersection, self.name))

719 r = self.Inverse( p.lat, p.lon, outmask=Caps.DISTANCE)

720 t = other.Inverse(p.lat, p.lon, outmask=Caps.DISTANCE)

721 P.set_(azi12= self.azi12, a12=r.a12, s12=r.s12,

722 azi02=other.azi12, a02=t.a12, s02=t.s12,

723 at=other.azi12 - self.azi12, iteration=i)

724 except Exception as x:

725 raise IntersectionError(self, other, tol=tol,

726 eps=eps, cause=x)

727 return P

728

729 def Inverse(self, lat2, lon2, wrap=False, **outmask):

730 '''Return the rhumb angle, distance, azimuth, I{reverse} azimuth, etc. of

731 a rhumb line between the given point and this rhumb line's start point.

732

733 @arg lat2: Latitude of the point (C{degrees}).

734 @arg lon2: Longitude of the points (C{degrees}).

735 @kwarg wrap: If C{True}, wrap or I{normalize} and unroll B{C{lat2}}

736 and B{C{lon2}} (C{bool}).

737

738 @return: L{GDict} with 8 items C{a12, s12, azi12, azi21, lat1, lon1,

739 lat2, lon2}, the rhumb angle C{a12} and rhumb distance C{s12}

740 between both points in C{degrees} respectively C{meter}, the

741 rhumb line's azimuth C{azi12} and I{reverse} azimuth C{azi21}

742 both in compass C{degrees} between C{-180} and C{+180}.

743 '''

744 if wrap:

745 _, lat2, lon2 = _Wrap.latlon3(self.lon1, _fix90(lat2), lon2, wrap)

746 r = self.rhumb.Inverse(self.lat1, self.lon1, lat2, lon2, **outmask)

747 return r

748

749 @Property_RO

750 def isLoxodrome(self):

751 '''Is this rhumb line a meridional (C{None}), a parallel

752 (C{False}) or a C{True} loxodrome?

753

754 @see: I{Osborne's} U{2.5 Rumb lines and loxodromes

755 <https://Zenodo.org/record/35392>}, page 37.

756 '''

757 return bool(self._salp) if self._calp else None

758

759 @Property_RO

760 def lat1(self):

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

762 '''

763 return self._lat1

764

765 @Property_RO

766 def lon1(self):

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

768 '''

769 return self._lon1

770

771 @Property_RO

772 def latlon1(self):

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

774 '''

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

776

777 def m2degrees(self, distance):

778 '''Convert a distance along this rhumb line to an angular distance.

779

780 @arg distance: Distance (C{meter}).

781

782 @return: Angular distance (C{degrees}).

783 '''

784 return _over(float(distance), self.rhumb._mpd)

785

786 @property_RO

787 def _mu1(self): # PYCHOK no cover

788 '''(INTERNAL) I{Must be overloaded}.'''

789 self._notOverloaded(underOK=True)

790

791 def _mu2lat(self, mu2): # PYCHOK no cover

792 '''(INTERNAL) I{Must be overloaded}.'''

793 self._notOverloaded(mu2) # underOK=True

794

795 @deprecated_method

796 def nearestOn4(self, lat0, lon0, **exact_eps_est_tol): # PYCHOK no cover

797 '''DEPRECATED on 23.10.10, use method L{PlumbTo}.'''

798 P = self.PlumbTo(lat0, lon0, **exact_eps_est_tol)

799 r = _MODS.deprecated.classes.NearestOn4Tuple(P.lat2, P.lon2, P.s12, P.azi02,

800 name=self.nearestOn4.__name__)

801 r._iteration = P.iteration

802 return r

803

804 @deprecated_method

805 def NearestOn(self, lat0, lon0, **exact_eps_est_tol): # PYCHOK no cover

806 '''DEPRECATED on 23.10.30, use method L{PlumbTo}.'''

807 return self.PlumbTo(lat0, lon0, **exact_eps_est_tol)

808

809 def PlumbTo(self, lat0, lon0, exact=None, eps=EPS, est=None, tol=_TOL):

810 '''Compute the I{perpendicular} intersection of this rhumb line with a geodesic

811 from the given point (transcoded from I{Karney}'s C++ U{rhumb-intercept

812 <https://SourceForge.net/p/geographiclib/discussion/1026620/thread/2ddc295e/>}).

813

814 @arg lat0: Latitude of the point on the geodesic (C{degrees}).

815 @arg lon0: Longitude of the point on the geodesic (C{degrees}).

816 @kwarg exact: If C{None}, use a rhumb line perpendicular to this rhumb line,

817 otherwise use an I{exact} C{Geodesic...} from the given point

818 perpendicular to this rhumb line (C{bool} or C{Geodesic...}),

819 see method L{Ellipsoid.geodesic_}.

820 @kwarg eps: Optional tolerance for L{pygeodesy.intersection3d3} (C{EPS}),

821 used only if C{B{exact} is None}.

822 @kwarg est: Optionally, an initial estimate for the distance C{s12} of the

823 intersection I{along} this rhumb line (C{meter}), used only if

824 C{B{exact} is not None}.

825 @kwarg tol: Longitudinal convergence tolerance (C{degrees}) or distance

826 tolerance (C(meter)) when C{B{exact} is None}, respectively

827 C{not None}.

828

829 @return: The intersection point on this rhumb line, a L{GDict} from method

830 L{Intersection} if B{C{exact}=None}. If C{B{exact} is not None},

831 a L{Position}-like L{GDict} of 13 items C{azi12, a12, s12, lat2,

832 lat1, lat0, lon2, lon1, lon0, azi0, a02, s02, at} with distance

833 C{a02} in C{degrees} and C{s02} in C{meter} between the given point

834 C{lat0, lon0} and the intersection C{lat2, lon2}, geodesic azimuth

835 C{azi0} at the given point and the (perpendicular) angle C{at}

836 between the geodesic and this rhumb line at the intersection. The

837 I{geodesic} azimuth at the intersection is C{(at + azi12)}. See

838 method L{Position} for further details.

839

840 @raise ImportError: I{Karney}'s U{geographiclib

841 <https://PyPI.org/project/geographiclib>}

842 package not found or not installed.

843

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

845 intersection for some other reason.

846

847 @see: Methods C{distance2}, C{Intersecant2} and C{Intersection}

848 and function L{pygeodesy.intersection3d3}.

849 '''

850 Cs, tol = Caps, Float_(tol=tol, low=EPS, high=None)

851

852# def _over(p, q): # see @note at method C{.Position}

853# if p:

854# p = (p / (q or _copysign(tol, q))) if isfinite(q) else NAN

855# return p

856

857 if exact is None:

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

859 rl = RhumbLineBase(self.rhumb, lat0, lon0, z, caps=Cs.LINE_OFF)

860 P = self.Intersection(rl, tol=tol, eps=eps)

861

862 else: # C{rhumb-intercept}

863 E = self.ellipsoid

864 _gI = E.geodesic_(exact=exact).Inverse

865 gm = Cs.STANDARD | Cs._REDUCEDLENGTH_GEODESICSCALE # ^ Cs.DISTANCE_IN

866 if est is None: # get an estimate from the "perpendicular" geodesic

867 r = _gI(self.lat1, self.lon1, lat0, lon0, outmask=Cs.AZIMUTH_DISTANCE)

868 d, _ = _diff182(r.azi2, self.azi12, K_2_0=True)

869 _, s12 = sincos2d(d)

870 s12 *= r.s12 # signed

871 else:

872 s12 = Meter(est=est)

873 try:

874 _abs = fabs

875 _d2 = _diff182

876 _ErT = E.rocPrimeVertical # aka rocTransverse

877 _ovr = _over

878 _S12 = Fsum(s12).fsum2f_

879 _scd = sincos2d_

880 for i in range(1, _TRIPS): # 9+, suffix 1 == C++ 2, 2 == C++ 3

881 P = self.Position(s12) # outmask=Cs.LATITUDE_LONGITUDE

882 r = _gI(lat0, lon0, P.lat2, P.lon2, outmask=gm)

883 d, _ = _d2(self.azi12, r.azi2, K_2_0=True)

884 s, c, s2, c2 = _scd(d, r.lat2)

885 c2 *= _ErT(r.lat2)

886 s *= _ovr(s2 * self._salp, c2) - _ovr(s * r.M21, r.m12)

887 s12, t = _S12(c / s) # XXX _ovr?

888 if _abs(t) < tol: # or _abs(c) < EPS

889 break

890 P.set_(azi0=r.azi1, a02=r.a12, s02=r.s12, # azi2=r.azi2,

891 lat0=lat0, lon0=lon0, iteration=i, at=r.azi2 - self.azi12,

892 name=_dunder_nameof(self.PlumbTo, self.name))

893 except Exception as x: # Fsum(NAN) Value-, ZeroDivisionError

894 raise IntersectionError(lat0, lon0, tol=tol, exact=exact,

895 eps=eps, est=est, iteration=i, cause=x)

896

897 return P

898

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

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

901

902 @arg s12: The distance along this rhumb line from its origin to

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

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

905 the quantities to be returned.

906

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

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

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

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

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

912

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

914 only required for L{RhumbLineAux} area C{S12}

915 when C{B{exact} is True}.

916

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

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

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

920

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

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

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

924 carried out in finite terms.

925 '''

926 return self._Position(self.m2degrees(s12), s12, outmask)

927

928 def _Position(self, a12, s12, outmask):

929 '''(INTERNAL) C{Arc-/Position} helper.

930 '''

931 r = GDict(azi12=self.azi12, a12=a12, s12=s12, name=self.name)

932 Cs = Caps

933 if (outmask & Cs.LATITUDE_LONGITUDE_AREA):

934 if a12 or s12:

935 mu12 = self._calp * a12

936 mu2 = self._mu1 + mu12

937 if fabs(mu2) > 90: # past pole

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

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

940 mu2 = _norm180(_loneg(mu2))

941 lat2 = self._mu2lat(mu2)

942 lon2 = S12 = NAN

943 else:

944 lat2, lon2, S1, S2 = self._Position4(a12, mu2, s12, mu12)

945 if (outmask & Cs.AREA):

946 S12 = self.rhumb._S12d(S1, S2, lon2)

947 S12 = unsigned0(S12) # like .gx

948# else:

949# S12 = None # unused

950 if (outmask & Cs.LONGITUDE):

951 if (outmask & Cs.LONG_UNROLL):

952 lon2 += self.lon1

953 else:

954 lon2 = _norm180(self._lon12 + lon2)

955 else: # coincident

956 lat2, lon2 = self.latlon1

957 S12 = _0_0

958

959 if (outmask & Cs.AREA):

960 r.set_(S12=S12)

961 if (outmask & Cs.LATITUDE):

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

963 if (outmask & Cs.LONGITUDE):

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

965 return r

966

967 def _Position4(self, a12, mu2, s12, mu12): # PYCHOK no cover

968 '''(INTERNAL) I{Must be overloaded}.'''

969 self._notOverloaded(a12, s12, mu2, mu12) # underOK=True

970

971 @Property_RO

972 def rhumb(self):

973 '''Get this rhumb line's rhumb (L{RhumbAux} or L{Rhumb}).

974 '''

975 return self._rhumb

976

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

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

979

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

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

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

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

984

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

986 '''

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

988 azi12=self.azi12, exact=self.exact,

989 TMorder=self.TMorder, xTM=self.xTM)

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

991

992 @property_RO

993 def TMorder(self):

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

995 '''

996 return self.rhumb.TMorder

997

998 @Property_RO

999 def xTM(self):

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

1001 if I{exact} and I{ellipsoidal}, otherwise L{KTransverseMercator} for C{TMorder}).

1002 '''

1003 E = self.ellipsoid

1004 # ExactTransverseMercator doesn't handle spherical earth models

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

1006 _MODS.ktm.KTransverseMercator(E, TMorder=self.TMorder)

1007

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

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

1010 as current intersection estimate and central meridian.

1011 '''

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

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

1014

1015

1016 class _PseudoRhumbLine(RhumbLineBase):

1017 '''(INTERNAL) Pseudo-rhumb line for a geodesic (line), see C{geodesicw._PlumbTo}.

1018 '''

1019 def __init__(self, gl, name=NN):

1020 R = RhumbBase(gl.geodesic.ellipsoid, None, True, name)

1021 RhumbLineBase.__init__(self, R, gl.lat1, gl.lon1, 0, caps=Caps.LINE_OFF)

1022 self._azi1 = self.azi12 = gl.azi1

1023 self._gl = gl

1024 self._gD = gl.geodesic.Direct

1025

1026 def PlumbTo(self, lat0, lon0, **exact_eps_est_tol): # PYCHOK signature

1027 P = RhumbLineBase.PlumbTo(self, lat0, lon0, **exact_eps_est_tol)

1028 z, P = _xkwds_pop2(P, azi12=None)

1029 P.set_(azi1=self._gl.azi1, azi2=z)

1030 return P # geodesic L{Position}

1031

1032 def Position(self, s12, **unused): # PYCHOK signature

1033 r = self._gD(self.lat1, self.lon1, self._azi1, s12)

1034 self._azi1 = r.azi1

1035 self.azi12 = z = r.azi2

1036 self._salp, _ = sincos2d(z)

1037 return r.set_(azi12=z)

1038

1039

1040__all__ += _ALL_DOCS(RhumbBase, RhumbLineBase)

1041

1042 if __name__ == '__main__':

1043

1044 from pygeodesy import printf, Rhumb as Rh, RhumbAux as Ah

1045 from pygeodesy.basics import _zip

1046 from pygeodesy.ellipsoids import _EWGS84

1047

1048 Al = Ah(_EWGS84).Line(30, 0, 45)

1049 Rl = Rh(_EWGS84).Line(30, 0, 45)

1050

1051 for i in range(1, 10):

1052 s = .5e6 + 1e6 / i

1053 a = Al.Position(s).lon2

1054 r = Rl.Position(s).lon2

1055 e = (fabs(a - r) / a) if a else 0

1056 printf('# Position.lon2 %.14f vs %.14f, diff %g', r, a, e)

1057

1058 for exact in (None, False, True):

1059 for est in (None, 1e6):

1060 a = Al.PlumbTo(60, 0, exact=exact, est=est)

1061 r = Rl.PlumbTo(60, 0, exact=exact, est=est)

1062 printf('# %s, iteration=%s, exact=%s, est=%s\n# %s, iteration=%s',

1063 a.toRepr(), a.iteration, exact, est,

1064 r.toRepr(), r.iteration, nl=1)

1065

1066 NE_=(71.688899882813, 0.2555198244234, 44095641862956.11)

1067 LHR=(77.7683897102557, 5771083.38332803, 37395209100030.39)

1068 NRT=(-92.38888798169965, 12782581.067684170, -63760642939072.50)

1069

1070 def _ref(fmt, r3, x3):

1071 e3 = []

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

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

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

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

1076

1077 for R in (Ah, Rh): # <https://GeographicLib.SourceForge.io/cgi-bin/RhumbSolve -p 9> version 2.2

1078 rh = R(exact=True) # WGS84 default

1079 printf('# %r', rh, nl=1)

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

1081 _ref('# JFK NE lat2=%.12f, lon2=%.12f, S12=%.1f', (r.lat2, r.lon2, r.S12), NE_)

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

1083 _ref('# JFK-LHR azi12=%.12f, s12=%.3f S12=%.1f', (r.azi12, r.s12, r.S12), LHR)

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

1085 _ref('# JFK-NRT azi12=%.12f, s12=%.3f S12=%.1f', (r.azi12, r.s12, r.S12), NRT)

1086

1087# % python3.10 -m pygeodesy3.rhumb.Bases

1088

1089# Position.lon2 11.61455846901637 vs 11.61455846901637, diff 3.05885e-16

1090# Position.lon2 7.58982302826842 vs 7.58982302826842, diff 2.34045e-16

1091# Position.lon2 6.28526067416369 vs 6.28526067416369, diff 2.82623e-16

1092# Position.lon2 5.63938995325146 vs 5.63938995325146, diff 1.57495e-16

1093# Position.lon2 5.25385527435707 vs 5.25385527435707, diff 0

1094# Position.lon2 4.99764604290380 vs 4.99764604290380, diff 8.88597e-16

1095# Position.lon2 4.81503363740473 vs 4.81503363740473, diff 1.84459e-16

1096# Position.lon2 4.67828821748836 vs 4.67828821748835, diff 5.69553e-16

1097# Position.lon2 4.57205667906283 vs 4.57205667906283, diff 5.82787e-16

1098

1099# Intersection(a02=17.798332, a12=19.521356, at=90.0, azi02=135.0, azi12=45.0, lat0=60.0, lat1=30.0, lat2=45.0, lon0=0.0, lon1=0.0, lon2=15.830286, name='Intersection', s02=1977981.142985, s12=2169465.957531), iteration=9, exact=None, est=None

1100# Intersection(a02=17.798332, a12=19.521356, at=90.0, azi02=135.0, azi12=45.0, lat0=60.0, lat1=30.0, lat2=45.0, lon0=0.0, lon1=0.0, lon2=15.830286, name='Intersection', s02=1977981.142985, s12=2169465.957531), iteration=9

1101

1102# Intersection(a02=17.798332, a12=19.521356, at=90.0, azi02=135.0, azi12=45.0, lat0=60.0, lat1=30.0, lat2=45.0, lon0=0.0, lon1=0.0, lon2=15.830286, name='Intersection', s02=1977981.142985, s12=2169465.957531), iteration=9, exact=None, est=1000000.0

1103# Intersection(a02=17.798332, a12=19.521356, at=90.0, azi02=135.0, azi12=45.0, lat0=60.0, lat1=30.0, lat2=45.0, lon0=0.0, lon1=0.0, lon2=15.830286, name='Intersection', s02=1977981.142985, s12=2169465.957531), iteration=9

1104

1105# PlumbTo(a02=17.967658, a12=27.74256, at=90.0, azi0=113.73626, azi12=45.0, lat0=60, lat1=30.0, lat2=49.634582, lon0=0, lon1=0.0, lon2=25.767876, name='PlumbTo', s02=1997960.116871, s12=3083112.636236), iteration=5, exact=False, est=None

1106# PlumbTo(a02=17.967658, a12=27.74256, at=90.0, azi0=113.73626, azi12=45.0, lat0=60, lat1=30.0, lat2=49.634582, lon0=0, lon1=0.0, lon2=25.767876, name='PlumbTo', s02=1997960.116871, s12=3083112.636236), iteration=5

1107

1108# PlumbTo(a02=17.967658, a12=27.74256, at=90.0, azi0=113.73626, azi12=45.0, lat0=60, lat1=30.0, lat2=49.634582, lon0=0, lon1=0.0, lon2=25.767876, name='PlumbTo', s02=1997960.116871, s12=3083112.636236), iteration=7, exact=False, est=1000000.0

1109# PlumbTo(a02=17.967658, a12=27.74256, at=90.0, azi0=113.73626, azi12=45.0, lat0=60, lat1=30.0, lat2=49.634582, lon0=0, lon1=0.0, lon2=25.767876, name='PlumbTo', s02=1997960.116871, s12=3083112.636236), iteration=7

1110

1111# PlumbTo(a02=17.967658, a12=27.74256, at=90.0, azi0=113.73626, azi12=45.0, lat0=60, lat1=30.0, lat2=49.634582, lon0=0, lon1=0.0, lon2=25.767876, name='PlumbTo', s02=1997960.116871, s12=3083112.636236), iteration=5, exact=True, est=None

1112# PlumbTo(a02=17.967658, a12=27.74256, at=90.0, azi0=113.73626, azi12=45.0, lat0=60, lat1=30.0, lat2=49.634582, lon0=0, lon1=0.0, lon2=25.767876, name='PlumbTo', s02=1997960.116871, s12=3083112.636236), iteration=5

1113

1114# PlumbTo(a02=17.967658, a12=27.74256, at=90.0, azi0=113.73626, azi12=45.0, lat0=60, lat1=30.0, lat2=49.634582, lon0=0, lon1=0.0, lon2=25.767876, name='PlumbTo', s02=1997960.116871, s12=3083112.636236), iteration=7, exact=True, est=1000000.0

1115# PlumbTo(a02=17.967658, a12=27.74256, at=90.0, azi0=113.73626, azi12=45.0, lat0=60, lat1=30.0, lat2=49.634582, lon0=0, lon1=0.0, lon2=25.767876, name='PlumbTo', s02=1997960.116871, s12=3083112.636236), iteration=7

1116

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

1118# JFK NE lat2=71.688899882813, lon2=0.255519824423, S12=44095641862956.1, rel errors: 4e-16, 2e-13, 4e-16

1119# JFK-LHR azi12=77.768389710256, s12=5771083.383 S12=37395209100030.4, rel errors: 5e-16, 3e-16, 8e-16

1120# JFK-NRT azi12=-92.388887981700, s12=12782581.068 S12=-63760642939072.5, rel errors: 0, 1e-16, 7e-16

1121

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

1123# JFK NE lat2=71.688899882813, lon2=0.255519824423, S12=44095641862956.1, rel errors: 2e-16, 1e-13, 5e-16

1124# JFK-LHR azi12=77.768389710256, s12=5771083.383 S12=37395209100030.4, rel errors: 4e-16, 3e-16, 6e-16

1125# JFK-NRT azi12=-92.388887981700, s12=12782581.068 S12=-63760642939072.5, rel errors: 0, 1e-16, 1e-16

1126

1127# **) MIT License

1128#

1129# Copyright (C) 2022-2024 -- mrJean1 at Gmail -- All Rights Reserved.

1130#

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

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

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

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

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

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

1137#

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

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

1140#

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

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

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

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

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

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

1147# OTHER DEALINGS IN THE SOFTWARE.