1

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

3

4 u'''Base classes C{RhumbBase} and C{RhumbLineBase}, pure Python version of I{Karney}'s C++

5 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{RhumbLine} has been enhanced with method C{Intersection} to iteratively find the intersection

11 of two rhumb lines and C{PlumbTo} to find the I{perpendicular} intersection of a rumb line and a

12 geodesic or rhumb line from a given point.

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, 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, _spherical_datum # from .formy

30# from pygeodesy.datums import _earth_datum, _WGS84 # from .karney

31 from pygeodesy.errors import IntersectionError, itemsorted, RhumbError, \

32 _xdatum, _xkwds, _xkwds_pop, _Xorder

33# from pygeodesy.etm import ExactTransverseMercator # _MODS

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

35 from pygeodesy.formy import opposing, Datum, _spherical_datum

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

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

38 _dunder_nameof, _parallel_, _too_, _under

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

40 _norm180, GDict, _earth_datum, _WGS84

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

42 from pygeodesy.lazily import _ALL_DOCS, _ALL_MODS as _MODS

43# from pygeodesy.named import notOverloaded # _MODS

44 from pygeodesy.namedTuples import Distance2Tuple, LatLon2Tuple

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

46 property_RO, _update_all

47 from pygeodesy.streprs import Fmt, pairs

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

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

50 _unrollon, _Wrap

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

52

53 from math import cos, fabs

54

55__all__ = ()

56__version__ = '23.11.09'

57

58_anti_ = _Dash('anti')

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

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

61

62

63 class _Lat(Lat):

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

65 '''

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

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

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

69

70

71 class _Lon(Lon):

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

73 '''

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

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

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

77

78

79 def _update_all_rls(r):

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

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

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

83 '''

84 # _xinstanceof(_MODS.rhumbaux.RhumbAux, _MODS.rhumbx.Rhumb, r=r)

85 _update_all(r)

86 for rl in _rls: # PYCHOK use weakref?

87 if rl._rhumb is r:

88 _update_all(rl)

89

90

91 class RhumbBase(_CapsBase):

92 '''(INTERNAL) Base class for C{rhumbaux.RhumbAux} and C{rhumbx.Rhumb}.

93 '''

94 _datum = _WGS84

95 _exact = True

96 _f_max = _0_01

97 _mTM = 6 # see .TMorder

98

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

100 '''New C{rhumbaux.RhumbAux} or C{rhumbx.Rhum}.

101 '''

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

103 if not exact:

104 self.exact = False

105 if name:

106 self.name = name

107

108 @Property_RO

109 def a(self):

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

111 '''

112 return self.ellipsoid.a

113

114 equatoradius = a

115

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

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

118

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

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

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

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

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

124

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

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

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

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

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

130

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

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

133 is True} and L{RhumbAux}.

134

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

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

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

138

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

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

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

142 carried out in finite terms.

143 '''

144 s12 = a12 * self._mpd

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

146

147 @Property_RO

148 def b(self):

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

150 '''

151 return self.ellipsoid.b

152

153 polaradius = b

154

155 @property

156 def datum(self):

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

158 '''

159 return self._datum

160

161 @datum.setter # PYCHOK setter!

162 def datum(self, datum):

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

164

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

166 '''

167 _xinstanceof(Datum, datum=datum)

168 if self._datum != datum:

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

170 _update_all_rls(self)

171 self._datum = datum

172

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

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

175 '''

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

177

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

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

180

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

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

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

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

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

186

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

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

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

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

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

192

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

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

195 is True} and L{RhumbAux}.

196

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

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

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

200

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

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

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

204 carried out in finite terms.

205 '''

206 a12 = _over(s12, self._mpd)

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

208

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

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

211 '''

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

213

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

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

216 '''

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

218

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

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

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

222

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

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

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

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

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

228 L{Caps} values specifying the required capabilities.

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

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

231

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

233 C{.Position} to compute each point.

234

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

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

237 '''

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

239

240 Line = DirectLine # synonyms

241

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

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

244 '''

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

246 name=self.name)

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

248

249 @Property

250 def ellipsoid(self):

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

252 '''

253 return self.datum.ellipsoid

254

255 @ellipsoid.setter # PYCHOK setter!

256 def ellipsoid(self, a_earth_f):

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

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

259

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

261 '''

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

263

264 @Property

265 def exact(self):

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

267 '''

268 return self._exact

269

270 @exact.setter # PYCHOK setter!

271 def exact(self, exact):

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

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

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

275

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

277 exceeds non-zero C{f_max}.

278

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

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

281 '''

282 x = bool(exact)

283 if self._exact != x:

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

285 _update_all_rls(self)

286 self._exact = x

287

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

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

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

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

292

293 @Property_RO

294 def f(self):

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

296 '''

297 return self.ellipsoid.f

298

299 flattening = f

300

301 @property

302 def f_max(self):

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

304 '''

305 return self._f_max

306

307 @f_max.setter # PYCHOK setter!

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

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

310

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

312 exceeds non-zero C{f_max}.

313 '''

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

315 if self._f_max != f:

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

317 self._f_max = f

318

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

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

321 '''

322 if wrap:

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

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

325

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

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

328

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

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

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

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

333

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

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

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

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

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

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

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

341 I{squared}.

342

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

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

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

346

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

348 on opposite meridians, there are two shortest rhumb lines

349 and the East-going one is chosen.

350

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

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

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

354 carried out in finite terms.

355 '''

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

357 Cs = Caps

358 if (outmask & Cs.AZIMUTH_DISTANCE_AREA):

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

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

361 if (outmask & Cs.AZIMUTH):

362 z = _atan2d(y, x)

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

364 if (outmask & Cs.AREA):

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

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

367 return r

368

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

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

371 _MODS.named.notOverloaded(self, lon12, r, Caps.toStr(outmask))

372

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

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

375 '''

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

377

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

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

380 '''

381 if wrap:

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

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

384

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

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

387

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

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

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

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

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

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

394 L{Caps} values specifying the required capabilities.

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

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

397

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

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

400

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

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

403 '''

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

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

406

407 @Property_RO

408 def _mpd(self): # PYCHOK no cover

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

410 _MODS.named.notOverloaded(self)

411

412 @property_RO

413 def RAorder(self):

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

415 '''

416 return None

417

418 @property_RO

419 def _RhumbLine(self): # PYCHOK no cover

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

421 _MODS.named.notOverloaded(self, underOK=True)

422

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

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

425 _MODS.named.notOverloaded(self, s1, s2, lon)

426

427 @Property

428 def TMorder(self):

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

430 '''

431 return self._mTM

432

433 @TMorder.setter # PYCHOK setter!

434 def TMorder(self, order):

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

436

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

438 for L{Rhumb} instances.

439 '''

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

441 if self._mTM != m:

442 _update_all_rls(self)

443 self._mTM = m

444 if self.exact and isinstance(self, _MODS.rhumbx.Rhumb):

445 self.exact = False

446

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

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

449

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

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

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

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

454

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

456 '''

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

458 exact=self.exact, TMorder=self.TMorder)

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

460

461

462 class RhumbLineBase(_CapsBase):

463 '''(INTERNAL) Base class for C{rhumbaux.RhumbLineAux} and C{rhumbx.RhumbLine}.

464 '''

465 _azi12 = _0_0

466 _calp = _1_0

467# _caps = 0

468# _debug = 0

469# _lat1 = _0_0

470# _lon1 = _0_0

471# _lon12 = _0_0

472 _Rhumb = RhumbBase # compatible C{Rhumb} class

473 _rhumb = None # C{Rhumb} instance

474 _salp = _0_0

475 _talp = _0_0

476

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

478 '''New C{RhumbLine}.

479 '''

480 _xinstanceof(self._Rhumb, rhumb=rhumb)

481

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

483 self._lon1 = _Lon(lon1= lon1)

484 self._lon12 = _norm180(self._lon1)

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

486 self.azi12 = _norm180(azi12)

487

488 n = name or rhumb.name

489 if n:

490 self.name=n

491

492 self._caps = caps

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

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

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

496 else:

497 self._rhumb = rhumb

498 _rls.append(self)

499

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

501 if _rls: # may be empty or None

502 try: # PYCHOK no cover

503 _rls.remove(self)

504 except (TypeError, ValueError):

505 pass

506 self._rhumb = None

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

508

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

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

511

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

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

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

515 the quantities to be returned.

516

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

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

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

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

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

522

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

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

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

526

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

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

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

530

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

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

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

534 carried out in finite terms.

535 '''

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

537

538 @Property

539 def azi12(self):

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

541 '''

542 return self._azi12

543

544 @azi12.setter # PYCHOK setter!

545 def azi12(self, azi12):

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

547 '''

548 z = _norm180(azi12)

549 if self._azi12 != z:

550 if self._rhumb:

551 _update_all(self)

552 self._azi12 = z

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

554 self._talp = _over(*t)

555

556 @property_RO

557 def azi12_sincos2(self): # PYCHOK no cover

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

559 '''

560 return self._scalp, self._calp

561

562 @property_RO

563 def datum(self):

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

565 '''

566 return self.rhumb.datum

567

568 def degrees2m(self, angle):

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

570

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

572

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

574 '''

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

576

577 @deprecated_method

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

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

580

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

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

583 '''

584 r = self.Inverse(lat, lon)

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

586

587 @property_RO

588 def ellipsoid(self):

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

590 '''

591 return self.rhumb.ellipsoid

592

593 @property_RO

594 def exact(self):

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

596 '''

597 return self.rhumb.exact

598

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

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

601

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

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

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

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

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

607 @kwarg tol_eps: Optional keyword arguments, see method

608 method L{Intersection} for further details.

609

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

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

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

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

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

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

616 rhumb line is tangential to the circle, both points

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

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

619

620 @raise IntersectionError: The circle and this rhumb line

621 do not intersect.

622

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

624 '''

625 r = Radius_(radius)

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

627 a = q.s02

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

629 if a < r:

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

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

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

633 d, h = _0_0, r

634 else:

635 d = q.s12

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

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

638 m = self.rhumb._mpr

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

640 else:

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

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

643 if h:

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

645 elif a > r:

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

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

648 txt=t, **tol_eps)

649 else: # tangential

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

651 return p, q

652

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

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

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

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

657 r._iteration = p.iteration

658 return r

659

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

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

662

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

664 @kwarg tol: Tolerance for longitudinal convergence and parallel

665 error (C{degrees}).

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

667

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

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

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

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

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

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

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

675 See method L{Position} for further details.

676

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

678 no intersection for an other reason.

679

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

681 L{pygeodesy.intersection3d3}.

682

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

684 L{ExactTransverseMercator} or L{KTransverseMercator}

685 projection and function L{pygeodesy.intersection3d3} in

686 that domain.

687 '''

688 _xinstanceof(RhumbLineBase, other=other)

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

690 try:

691 if self.others(other) is self:

692 raise ValueError(_coincident_)

693 # make invariants and globals locals

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

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

696 p = opposing(s_az, o_az, margin=tol)

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

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

699 _diff = euclid # approximate length

700 _i3d3 = _intersect3d3 # NOT .vector3d.intersection3d3

701 _LL2T = LatLon2Tuple

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

703 # use halfway point as initial estimate

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

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

706 for i in range(1, _TRIPS):

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

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

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

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

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

712 if d < tol: # 19 trips

713 break

714 else:

715 raise ValueError(Fmt.no_convergence(d))

716

717 n = _dunder_nameof(self.Intersection, self.name)

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

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

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

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

722 azi12= self.azi12, a12=r.a12, s12=r.s12,

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

724 at=other.azi12 - self.azi12, name=n)

725 P._iteration = i # .set_(iteration=i, ...) only

726 except Exception as x:

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

728 eps=eps, cause=x)

729 return P

730

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

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

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

734

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

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

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

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

739

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

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

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

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

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

745 '''

746 if wrap:

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

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

749 return r

750

751 @Property_RO

752 def isLoxodrome(self):

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

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

755

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

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

758 '''

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

760

761 @Property_RO

762 def lat1(self):

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

764 '''

765 return self._lat1

766

767 @Property_RO

768 def lon1(self):

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

770 '''

771 return self._lon1

772

773 @Property_RO

774 def latlon1(self):

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

776 '''

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

778

779 def m2degrees(self, distance):

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

781

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

783

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

785 '''

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

787

788 @property_RO

789 def _mu1(self): # PYCHOK no cover

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

791 _MODS.named.notOverloaded(self, underOK=True)

792

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

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

795 _MODS.named.notOverloaded(self, mu2, underOK=True)

796

797 @deprecated_method

798 def nearestOn4(self, lat0, lon0, **exact_eps_est_tol):

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

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

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

802 name=self.nearestOn4.__name__)

803 r._iteration = P.iteration

804 return r

805

806 @deprecated_method

807 def NearestOn(self, lat0, lon0, **exact_eps_est_tol):

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

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

810

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

812 '''Compute the I{perpendicular} intersection of this rumb line with a geodesic

813 from the given point, in part transcoded from I{Karney}'s C++ U{rhumb-intercept

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

815

816 @arg lat0: Latitude of the point (C{degrees}).

817 @arg lon0: Longitude of the point (C{degrees}).

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

819 line, otherwise use an I{exact} C{Geodesic...} from the

820 given point perpendicular to this rhumb line (C{bool} or

821 C{Geodesic...}), see method L{Ellipsoid.geodesic_}.

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

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

824 @kwarg est: Optional, initial estimate for the distance C{s12} of the

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

826 if C{B{exact} is not None}.

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

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

829 C{not None}.

830

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

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

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

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

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

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

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

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

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

840 method L{Position} for further details.

841

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

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

844 package not found or not installed.

845

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

847 intersection for some other reason.

848

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

850 and function L{pygeodesy.intersection3d3}.

851 '''

852 Cs = Caps

853 if exact is None:

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

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

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

857

858 else: # C{rhumb-intercept}

859 E = self.ellipsoid

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

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

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

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

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

865 _, s12 = sincos2d(d)

866 s12 *= r.s12 # signed

867 else:

868 s12 = Meter(est=est)

869 try:

870 tol = Float_(tol=tol, low=EPS, high=None)

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

872 # return (p / (q or _copysign(tol, q))) if isfinite(q) else NAN

873

874 _ErT = E.rocPrimeVertical # aka rocTransverse

875 _S12 = Fsum(s12).fsum2_

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

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

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

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

880 s, c, s2, c2 = sincos2d_(d, r.lat2)

881 c2 *= _ErT(r.lat2)

882 s *= _over(s2 * self._salp, c2) - _over(s * r.M21, r.m12)

883 s12, t = _S12(c / s) # XXX _over?

884 if fabs(t) < tol: # or fabs(c) < EPS

885 break

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

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

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

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

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

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

892

893 return P

894

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

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

897

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

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

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

901 the quantities to be returned.

902

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

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

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

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

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

908

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

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

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

912

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

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

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

916

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

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

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

920 carried out in finite terms.

921 '''

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

923

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

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

926 '''

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

928 Cs = Caps

929 if (outmask & Cs.LATITUDE_LONGITUDE_AREA):

930 if a12 or s12:

931 mu12 = self._calp * a12

932 mu2 = self._mu1 + mu12

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

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

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

936 mu2 = _norm180(_loneg(mu2))

937 lat2 = self._mu2lat(mu2)

938 lon2 = S12 = NAN

939 else:

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

941 if (outmask & Cs.AREA):

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

943 S12 = unsigned0(S12) # like .gx

944# else:

945# S12 = None # unused

946 if (outmask & Cs.LONGITUDE):

947 if (outmask & Cs.LONG_UNROLL):

948 lon2 += self.lon1

949 else:

950 lon2 = _norm180(self._lon12 + lon2)

951 else: # coincident

952 lat2, lon2 = self.latlon1

953 S12 = _0_0

954

955 if (outmask & Cs.AREA):

956 r.set_(S12=S12)

957 if (outmask & Cs.LATITUDE):

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

959 if (outmask & Cs.LONGITUDE):

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

961 return r

962

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

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

965 _MODS.named.notOverloaded(self, a12, s12, mu2, mu12)

966

967 @Property_RO

968 def rhumb(self):

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

970 '''

971 return self._rhumb

972

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

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

975

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

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

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

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

980

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

982 '''

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

984 azi12=self.azi12, exact=self.exact,

985 TMorder=self.TMorder, xTM=self.xTM)

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

987

988 @property_RO

989 def TMorder(self):

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

991 '''

992 return self.rhumb.TMorder

993

994 @Property_RO

995 def xTM(self):

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

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

998 '''

999 E = self.ellipsoid

1000 # ExactTransverseMercator doesn't handle spherical earth models

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

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

1003

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

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

1006 as current intersection estimate and central meridian.

1007 '''

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

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

1010

1011

1012 class _PseudoRhumbLine(RhumbLineBase):

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

1014 '''

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

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

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

1018 self._azi1 = self.azi12 = gl.azi1

1019 self._gl = gl

1020 self._gD = gl.geodesic.Direct

1021

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

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

1024 P.set_(azi1=self._gl.azi1, azi2=_xkwds_pop(P, azi12=None))

1025 return P # geodesic L{Position}

1026

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

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

1029 self._azi1 = r.azi1

1030 self.azi12 = z = r.azi2

1031 self._salp, _ = sincos2d(z)

1032 return r.set_(azi12=z)

1033

1034

1035__all__ += _ALL_DOCS(RhumbBase, RhumbLineBase)

1036

1037 if __name__ == '__main__':

1038

1039 from pygeodesy import printf, Rhumb as R, RhumbAux as A

1040 from pygeodesy.ellipsoids import _EWGS84

1041

1042 A = A(_EWGS84).Line(30, 0, 45)

1043 R = R(_EWGS84).Line(30, 0, 45)

1044

1045 for i in range(1, 10):

1046 s = .5e6 + 1e6 / i

1047 a = A.Position(s).lon2

1048 r = R.Position(s).lon2

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

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

1051

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

1053 for est in (None, 1e6):

1054 a = A.PlumbTo(60, 0, exact=exact, est=est)

1055 r = R.PlumbTo(60, 0, exact=exact, est=est)

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

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

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

1059

1060# % python3 -m pygeodesy.rhumbBase

1061

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

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

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

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

1066# Position.lon2 5.25385527435707 vs 5.25385527435707, diff 0

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

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

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

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

1071

1072# 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

1073# 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

1074

1075# 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

1076# 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

1077

1078# 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

1079# 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

1080

1081# 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

1082# 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

1083

1084# 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

1085# 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

1086

1087# 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

1088# 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

1089

1090# **) MIT License

1091#

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

1093#

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

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

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

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

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

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

1100#

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

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

1103#

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

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

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

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

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

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

1110# OTHER DEALINGS IN THE SOFTWARE.