1

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

3

4 u'''Wrapper around Python classes C{geodesic.Geodesic} and C{geodesicline.GeodesicLine} from

5 I{Karney}'s Python package U{geographiclib<https://PyPI.org/project/geographiclib>}, provided

6 that package is installed.

7

8 The I{wrapped} class methods return a L{GDict} instance offering access to the C{dict} items

9 either by C{key} or by C{attribute} name.

10

11 With env variable C{PYGEODESY_GEOGRAPHICLIB} left undefined or set to C{"2"}, this module,

12 L{pygeodesy.geodesicx} and L{pygeodesy.karney} will use U{GeographicLib 2.0

13<https://GeographicLib.SourceForge.io/C++/doc/>} transcoding, otherwise C{1.52} or older.

14'''

15

16 from pygeodesy.basics import _copysign, _xinstanceof

17 from pygeodesy.constants import EPS, NAN, _EPSqrt as _TOL, _0_5

18 from pygeodesy.datums import _earth_datum, _WGS84, _EWGS84

19# from pygeodesy.dms import F_D # from .latlonBase

20# from pygeodesy.ellipsoids import _EWGS84 # from .datums

21 from pygeodesy.errors import IntersectionError, GeodesicError

22 from pygeodesy.fsums import Fsum, Fmt, unstr

23 from pygeodesy.internals import _dunder_nameof, _under

24 from pygeodesy.interns import NN, _DOT_, _SPACE_, _to_, _too_

25 from pygeodesy.karney import _atan2d, Caps, Direct9Tuple, GDict, \

26 Inverse10Tuple, _kWrapped, _llz2line

27 from pygeodesy.latlonBase import LatLonBase as _LLB, F_D, Radius_

28 from pygeodesy.lazily import _ALL_LAZY, _ALL_MODS as _MODS

29 from pygeodesy.named import callername, classname, _name1__, _name2__

30 from pygeodesy.namedTuples import Destination3Tuple, Distance3Tuple

31 from pygeodesy.props import Property, Property_RO, property_RO

32# from pygeodesy.streprs import Fmt, unstr # from .fsums

33# from pygeodesy.units import Radius_ # from .latlonBase

34 from pygeodesy.utily import _unrollon, _Wrap, wrap360, fabs # PYCHOK used!

35

36 from contextlib import contextmanager

37# from math import fabs # from .utily

38

39__all__ = _ALL_LAZY.geodesicw

40__version__ = '24.07.09'

41

42_plumb_ = 'plumb'

43_TRIPS = 65

44

45

46 class _gWrapped(_kWrapped):

47 ''''(INTERNAL) Wrapper for some of I{Karney}'s U{geographiclib

48 <https://PyPI.org/project/geographiclib>} classes.

49 '''

50

51 @Property_RO # MCCABE 24

52 def Geodesic(self):

53 '''Get the I{wrapped} C{geodesic.Geodesic} class from I{Karney}'s Python

54 U{geographiclib<https://GitHub.com/geographiclib/geographiclib-python>},

55 provided the latter is installed.

56 '''

57 _Geodesic = self.geographiclib.Geodesic

58 # assert Caps._STD == _Geodesic.STANDARD

59

60 class Geodesic(_Geodesic):

61 '''I{Wrapper} for I{Karney}'s Python U{geodesic.Geodesic

62 <https://PyPI.org/project/geographiclib>} class.

63 '''

64 _datum = _WGS84

65 _debug = 0 # like .geodesicx.bases._GeodesicBase

66 LINE_OFF = 0 # in .azimuthal._GnomonicBase and .css.CassiniSoldner

67 _name = NN

68

69 def __init__(self, a_ellipsoid=_EWGS84, f=None, **name): # PYCHOK signature

70 '''New I{wrapped} C{geodesic.Geodesic} instance.

71

72 @arg a_ellipsoid: The equatorial radius I{a} (C{meter}, conventionally),

73 an ellipsoid (L{Ellipsoid}) or a datum (L{Datum}).

74 @arg f: The ellipsoid's flattening (C{scalar}), required if B{C{a_ellipsoid})

75 is C{meter}, ignored otherwise.

76 @kwarg name: Optional C{B{name}=NN} (C{str}).

77 '''

78 _earth_datum(self, a_ellipsoid, f=f, **name) # raiser=NN

79 E = self.ellipsoid

80 with _wargs(self, *E.a_f, **name) as args:

81 _Geodesic.__init__(self, *args)

82 if name:

83 self._name, _ = _name2__(name, _or_nameof=E)

84

85 def ArcDirect(self, lat1, lon1, azi1, a12, outmask=Caps._STD):

86 '''Return the C{_Geodesic.ArcDirect} result as L{GDict}.

87 '''

88 with _wargs(self, lat1, lon1, azi1, a12, outmask) as args:

89 d = _Geodesic.ArcDirect(self, *args)

90 return GDict(d)

91

92 def ArcDirectLine(self, lat1, lon1, azi1, a12, caps=Caps._STD_LINE, **name):

93 '''Return the C{_Geodesic.ArcDirectLine} as I{wrapped} C{GeodesicLine}.

94 '''

95 return self._GenDirectLine(lat1, lon1, azi1, True, a12, caps, **name)

96

97 Area = _Geodesic.Polygon # like GeodesicExact.Area

98

99 @property_RO

100 def datum(self):

101 '''Get this geodesic's datum (C{Datum}).

102 '''

103 return self._datum

104

105 @Property

106 def debug(self):

107 '''Get the C{debug} option (C{bool}).

108 '''

109 return bool(self._debug)

110

111 @debug.setter # PYCHOK setter!

112 def debug(self, debug):

113 '''Set the C{debug} option (C{bool}) to include more

114 details in L{GDict} results.

115 '''

116 self._debug = Caps._DEBUG_ALL if debug else 0

117

118 def Direct(self, lat1, lon1, azi1, s12=0, outmask=Caps._STD):

119 '''Return the C{_Geodesic.Direct} result as L{GDict}.

120 '''

121 with _wargs(self, lat1, lon1, azi1, s12, outmask) as args:

122 d = _Geodesic.Direct(self, *args)

123 return GDict(d)

124

125 def Direct3(self, lat1, lon1, azi1, s12): # PYCHOK outmask

126 '''Return the destination lat, lon and reverse azimuth

127 in C{degrees} as L{Destination3Tuple}.

128 '''

129 d = self.Direct(lat1, lon1, azi1, s12, outmask=Caps._DIRECT3)

130 return Destination3Tuple(d.lat2, d.lon2, d.azi2)

131

132 def _DirectLine(self, ll1, azi12, s12=0, **caps_name):

133 '''(INTERNAL) Short-cut version.

134 '''

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

136

137 def DirectLine(self, lat1, lon1, azi1, s12, caps=Caps._STD_LINE, **name):

138 '''Return the C{_Geodesic.DirectLine} as I{wrapped} C{GeodesicLine}.

139 '''

140 return self._GenDirectLine(lat1, lon1, azi1, False, s12, caps, **name)

141

142 @Property_RO

143 def ellipsoid(self):

144 '''Get this geodesic's ellipsoid (C{Ellipsoid}).

145 '''

146 return self.datum.ellipsoid

147

148 @property_RO

149 def f1(self): # in .css.CassiniSoldner.reset

150 '''Get the geodesic's ellipsoid's I{1 - flattening} (C{float}).

151 '''

152 return getattr(self, _under(Geodesic.f1.name), self.ellipsoid.f1)

153

154 def _GDictDirect(self, lat, lon, azi, arcmode, s12_a12, outmask=Caps._STD):

155 '''(INTERNAL) Get C{_Geodesic._GenDirect} result as C{GDict}.

156 '''

157 with _wargs(self, lat, lon, azi, arcmode, s12_a12, outmask) as args:

158 t = _Geodesic._GenDirect(self, *args)

159 return Direct9Tuple(t).toGDict() # *t

160

161 def _GDictInverse(self, lat1, lon1, lat2, lon2, outmask=Caps._STD):

162 '''(INTERNAL) Get C{_Geodesic._GenInverse} result as L{Inverse10Tuple}.

163 '''

164 with _wargs(self, lat1, lon1, lat2, lon2, outmask) as args:

165 t = _Geodesic._GenInverse(self, *args)

166 return Inverse10Tuple(t).toGDict(lon1=lon1, lon2=lon2) # *t

167

168 def _GenDirectLine(self, lat1, lon1, azi1, arcmode, s12_a12, *caps, **name):

169 '''(INTERNAL) Invoked by C{_Geodesic.DirectLine} and C{-.ArcDirectLine},

170 returning the result as a I{wrapped} C{GeodesicLine}.

171 '''

172 with _wargs(self, lat1, lon1, azi1, arcmode, s12_a12, *caps, **name) as args:

173 t = _Geodesic._GenDirectLine(self, *args)

174 return self._Line13(t, **name)

175

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

177 '''(INTERNAL) Short-cut version, see .ellipsoidalBaseDI.intersecant2.

178 '''

179 if wrap:

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

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

182

183 def Inverse(self, lat1, lon1, lat2, lon2, outmask=Caps._STD):

184 '''Return the C{_Geodesic.Inverse} result as L{GDict}.

185 '''

186 with _wargs(self, lat1, lon1, lat2, lon2, outmask) as args:

187 d = _Geodesic.Inverse(self, *args)

188 return GDict(d)

189

190 def Inverse1(self, lat1, lon1, lat2, lon2, wrap=False):

191 '''Return the non-negative, I{angular} distance in C{degrees}.

192

193 @kwarg wrap: If C{True}, wrap or I{normalize} and unroll

194 B{C{lat2}} and BC{lon2}} (C{bool}).

195 '''

196 # see .FrechetKarney.distance, .HausdorffKarney._distance

197 # and .HeightIDWkarney._distances

198 if wrap:

199 _, lat2, lon2 = _Wrap.latlon3(lat1, lat2, lon2, True) # _Geodesic.LONG_UNROLL

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

201 # XXX _Geodesic.DISTANCE needed for 'a12'?

202 return fabs(r.a12)

203

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

205 '''Return the distance in C{meter} and the forward and reverse

206 azimuths in C{degrees} as L{Distance3Tuple}.

207 '''

208 r = self.Inverse(lat1, lon1, lat2, lon2, outmask=Caps._INVERSE3)

209 return Distance3Tuple(r.s12, wrap360(r.azi1), wrap360(r.azi2))

210

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

212 '''(INTERNAL) Short-cut version.

213 '''

214 if wrap:

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

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

217

218 def InverseLine(self, lat1, lon1, lat2, lon2, caps=Caps._STD_LINE, **name):

219 '''Return the C{_Geodesic.InverseLine} as I{wrapped} C{GeodesicLine}.

220 '''

221 with _wargs(self, lat1, lon1, lat2, lon2, caps, **name) as args:

222 t = _Geodesic.InverseLine(self, *args)

223 gl = self._Line13(t, **name)

224 p2 = gl.Position(gl.s13, outmask=Caps.AZIMUTH)

225 return _llz2line(gl, lat2=lat2, lon2=lon2, azi2=p2.azi2)

226

227 def Line(self, lat1, lon1, azi1, caps=Caps._STD_LINE, **name):

228 '''Set up a I{wrapped} C{GeodesicLine} to compute several points

229 along a single, I{wrapped} (this) geodesic.

230 '''

231 return _wrapped.GeodesicLine(self, lat1, lon1, azi1, caps=caps, **name)

232

233 def _Line13(self, t, **name):

234 '''(INTERNAL) Wrap C{_GeodesicLine}, add distance and arc length

235 to reference point 3.

236 '''

237 gl = _wrapped.GeodesicLine(self, t.lat1, t.lon1, t.azi1, caps=t.caps,

238 salp1=t.salp1, calp1=t.calp1, **name)

239 gl.a13, gl.s13 = t.a13, t.s13

240 return gl

241

242 @property_RO

243 def name(self):

244 '''Get the name (C{str}).

245 '''

246 return self._name

247

248# Polygon = _Geodesic.Polygon

249

250 # Geodesic.ArcDirect.__doc__ = _Geodesic.ArcDirect.__doc__

251 # Geodesic.Direct.__doc__ = _Geodesic.Direct.__doc__

252 # Geodesic.Inverse.__doc__ = _Geodesic.Inverse.__doc__

253 # Geodesic.InverseLine.__doc__ = _Geodesic.InverseLinr.__doc__

254 # Geodesic.Line.__doc__ = _Geodesic.Line.__doc__

255 return Geodesic

256

257 @Property_RO # MCCABE 16

258 def GeodesicLine(self):

259 '''Get the I{wrapped} C{geodesicline.GeodesicLine} class from I{Karney}'s

260 Python U{geographiclib<https://GitHub.com/geographiclib/geographiclib-python>},

261 provided the latter is installed.

262 '''

263 _GeodesicLine = self.geographiclib.GeodesicLine

264

265 class GeodesicLine(_GeodesicLine):

266 '''I{Wrapper} for I{Karney}'s Python U{geodesicline.GeodesicLine

267 <https://PyPI.org/project/geographiclib>} class.

268 '''

269 _geodesic = None

270 _name = NN

271

272 def __init__(self, geodesic, lat1, lon1, azi1, **caps_name_): # salp1=NAN, calp1=NAN

273 '''New I{wrapped} C{geodesicline.GeodesicLine} instance.

274

275 @arg geodesic: A I{wrapped} C{Geodesic} instance.

276 @arg lat1: Latitude of the first points (C{degrees}).

277 @arg lon1: Longitude of the first points (C{degrees}).

278 @arg azi1: Azimuth at the first points (compass C{degrees360}).

279 @kwarg caps_name_: Optional keyword arguments C{B{caps}=Caps.STANDARD},

280 a bit-or'ed combination of L{Caps} values specifying the

281 capabilities the C{GeodesicLine} instance should possess,

282 an optional C{B{name}=NN} plus C{salp1=NAN} and C{calp1=NAN}

283 for I{INTERNAL} use.

284 '''

285 _xinstanceof(_wrapped.Geodesic, geodesic=geodesic)

286 with _wargs(self, geodesic, lat1, lon1, azi1, **caps_name_) as args:

287 name, caps_ = _name2__(caps_name_, _or_nameof=geodesic)

288 _GeodesicLine.__init__(self, *args, **caps_) # XXX avoid updates?

289 if name:

290 self._name = name

291 self._geodesic = geodesic

292

293 @Property_RO

294 def a1(self):

295 '''Get the I{equatorial arc} (C{degrees}), the arc length between

296 the northward equatorial crossing and point C{(lat1, lon1)}.

297

298 @see: U{EquatorialArc<https://GeographicLib.SourceForge.io/

299 C++/doc/classGeographicLib_1_1GeodesicLine.html>}

300 '''

301 try:

302 return _atan2d(self._ssig1, self._csig1)

303 except AttributeError:

304 return NAN # see .geodesicx.gxline._GeodesicLineExact

305

306 equatorarc = a1

307

308 def Arc(self):

309 '''Return the angular distance to point 3 (C{degrees} or C{NAN}).

310 '''

311 return self.a13

312

313 def ArcPosition(self, a12, outmask=Caps._STD):

314 '''Return the position at C{B{a12} degrees} on this line.

315

316 @arg a12: Angular distance from this line's first point

317 (C{degrees}).

318

319 @see: Method L{Position} for further details.

320 '''

321 with _wargs(self, a12, outmask) as args:

322 d = _GeodesicLine.ArcPosition(self, *args)

323 return GDict(d)

324

325 @Property_RO

326 def azi0(self): # see .css.CassiniSoldner.forward4

327 '''Get the I{equatorial azimuth} (C{degrees}), the azimuth of the

328 geodesic line as it crosses the equator in a northward direction.

329

330 @see: U{EquatorialAzimuth<https://GeographicLib.SourceForge.io/

331 C++/doc/classGeographicLib_1_1GeodesicLine.html>}

332 '''

333 try:

334 return _atan2d(self._salp0, self._calp0)

335 except AttributeError:

336 return NAN # see .geodesicx.gxline._GeodesicLineExact

337

338 equatorazimuth = azi0

339

340 def Distance(self):

341 '''Return the distance to reference point 3 (C{meter} or C{NAN}).

342 '''

343 return self.s13

344

345 @property_RO

346 def geodesic(self):

347 '''Get the I{wrapped} geodesic (L{Geodesic}).

348 '''

349 return self._geodesic

350

351 def Intersecant2(self, lat0, lon0, radius, tol=_TOL):

352 '''Compute the intersection(s) of this geodesic line and a circle.

353

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

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

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

357 @kwarg tol: Convergence tolerance (C{scalar}).

358

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

360 a geodesic chord), each a L{GDict} from method L{Position} and

361 extended to 14 items C{lat1, lon1, azi1, lat2, lon2, azi2, a12,

362 s12, lat0, lon0, azi0, a02, s02, at} with the circle center

363 C{lat0}, C{lon0}, azimuth C{azi0} at the intersection, distance

364 C{a02} in C{degrees} and C{s02} in C{meter} along the geodesic

365 from the circle center to the intersection C{lat2, lon2} and

366 the angle C{at} between the geodesic and this line at the

367 intersection. The I{geodesic} azimuth at the intersection is

368 C{(at + azi2)}. If this line is tangential to the circle, both

369 intersections are the same L{GDict} instance.

370

371 @raise IntersectionError: The circle and this geodesic line do not

372 intersect.

373

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

375 '''

376 return _Intersecant2(self, lat0, lon0, radius, tol=tol)

377

378 def PlumbTo(self, lat0, lon0, est=None, tol=_TOL):

379 '''Compute the I{perpendicular} intersection of this geodesic line

380 with a geodesic from the given point.

381

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

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

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

385 the intersection I{along} this geodesic line (C{meter}).

386 @kwarg tol: Convergence tolerance (C(meter)).

387

388 @return: The intersection point on this geodesic line, a L{GDict}

389 from method L{Position} extended to 14 items C{lat1, lon1,

390 azi1, lat2, lon2, azi2, a12, s12, lat0, lon0, azi0, a02,

391 s02, at} with C{a02} and C{s02} the distance in C{degrees}

392 and C{meter} from the given point C{lat0, lon0} to the

393 intersection C{lat2, lon2}, azimuth C{azi0} at the given

394 point and the (perpendicular) angle C{at} between the

395 geodesic and this line at the intersection point. The

396 geodesic azimuth at the intersection is C{(at + azi2)}.

397 See method L{Position} for further details.

398

399 @see: Methods C{Intersecant2}, C{Intersection} and C{Position}.

400 '''

401 return _PlumbTo(self, lat0, lon0, est=est, tol=tol)

402

403 def Position(self, s12, outmask=Caps._STD):

404 '''Return the position at distance C{B{s12} meter} on this line.

405

406 @arg s12: Distance from this line's first point (C{meter}).

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

408 the quantities to be returned.

409

410 @return: A L{GDict} with up to 12 items C{lat1, lon1, azi1, lat2,

411 lon2, azi2, m12, a12, s12, M12, M21, S12} with C{lat1},

412 C{lon1}, C{azi1} and arc length C{a12} always included,

413 except when C{a12=NAN}.

414 '''

415 with _wargs(self, s12, outmask) as args:

416 d = _GeodesicLine.Position(self, *args)

417 return GDict(d)

418

419 # GeodesicLine.ArcPosition.__doc__ = _GeodesicLine.ArcPosition.__doc__

420 # GeodesicLine.Position.__doc__ = _GeodesicLine.Position.__doc__

421 return GeodesicLine

422

423 @Property_RO

424 def Geodesic_WGS84(self):

425 '''Get the I{wrapped} C{Geodesic(WGS84)} singleton, provided the

426 U{geographiclib<https://PyPI.org/project/geographiclib>} package

427 is installed, otherwise an C{ImportError}.

428 '''

429 return _EWGS84.geodesic

430

431_wrapped = _gWrapped() # PYCHOK singleton, .ellipsoids, .test/base.py

432

433

434 def Geodesic(a_ellipsoid=_EWGS84, f=None, **name):

435 '''Return a I{wrapped} C{geodesic.Geodesic} instance from I{Karney}'s

436 Python U{geographiclib<https://PyPI.org/project/geographiclib>},

437 provide the latter is installed, otherwise an C{ImportError}.

438

439 @arg a_ellipsoid: An ellipsoid (L{Ellipsoid}) or datum (L{Datum})

440 or the equatorial radius I{a} of the ellipsoid (C{meter}).

441 @arg f: The flattening of the ellipsoid (C{scalar}), required if

442 B{C{a_ellipsoid}}) is C{meter}, ignored otherwise.

443 @kwarg name: Optional C{B{name}=NN} (C{str}), ignored like B{C{f}}.

444 '''

445 return _wrapped.Geodesic(a_ellipsoid, f=f, **name)

446

447

448 def GeodesicLine(geodesic, lat1, lon1, azi1, caps=Caps._STD_LINE):

449 '''Return a I{wrapped} C{geodesicline.GeodesicLine} instance from I{Karney}'s

450 Python U{geographiclib<https://PyPI.org/project/geographiclib>}, provided

451 the latter is installed, otherwise an C{ImportError}.

452

453 @arg geodesic: A I{wrapped} L{Geodesic} instance.

454 @arg lat1: Latitude of the first points (C{degrees}).

455 @arg lon1: Longitude of the first points (C{degrees}).

456 @arg azi1: Azimuth at the first points (compass C{degrees360}).

457 @kwarg caps: Optional, bit-or'ed combination of L{Caps} values specifying

458 the capabilities the C{GeodesicLine} instance should possess,

459 i.e., which quantities can be returned by methods

460 C{GeodesicLine.Position} and C{GeodesicLine.ArcPosition}.

461 '''

462 return _wrapped.GeodesicLine(geodesic, lat1, lon1, azi1, caps=caps)

463

464

465 def Geodesic_WGS84():

466 '''Get the I{wrapped} L{Geodesic}C{(WGS84)} singleton, provided

467 U{geographiclib<https://PyPI.org/project/geographiclib>} is

468 installed, otherwise an C{ImportError}.

469 '''

470 return _wrapped.Geodesic_WGS84

471

472

473 class _wargs(object): # see also .formy._idllmn6, .latlonBase._toCartesian3, .vector2d._numpy

474 '''(INTERNAL) C{geographiclib} caller and exception mapper.

475 '''

476 @contextmanager # <https://www.Python.org/dev/peps/pep-0343/> Examples

477 def __call__(self, inst, *args, **kwds):

478 '''(INTERNAL) Yield C{tuple(B{args})} with any errors raised

479 as L{GeodesicError} embellished with all B{C{kwds}}.

480 '''

481 try:

482 yield args

483 except Exception as x:

484 u = _DOT_(classname(inst), callername(up=2, underOK=True))

485 raise GeodesicError(unstr(u, *args, **_name1__(kwds)), cause=x)

486

487_wargs = _wargs() # PYCHOK singleton

488

489

490 def _Intersecant2(gl, lat0, lon0, radius, tol=_TOL, form=F_D): # MCCABE in LatLonEllipsoidalBaseDI.intersecant2, .geodesicx.gxline.Intersecant2

491 # (INTERNAL) Return the intersections of a circle at C{lat0, lon0}

492 # and a geodesic line as a 2-Tuple C{(P, Q)}, each a C{GDict}.

493 r = Radius_(radius)

494 n = _dunder_nameof(_Intersecant2)[1:]

495 _P = gl.Position

496 _I = gl.geodesic.Inverse

497 _a = fabs

498

499 def _R3(s):

500 # radius, intersection, etc. at distance C{s}

501 P = _P(s)

502 d = _I(lat0, lon0, P.lat2, P.lon2)

503 return _a(d.s12), P, d

504

505 def _bisect2(s, c, Rc, r, tol, _R3):

506 _s = Fsum(c).fsumf_

507 for i in range(_TRIPS):

508 b = _s(s)

509 Rb, P, d = _R3(b)

510 if Rb > r:

511 break

512 else: # b >>> s and c >>> s

513 raise ValueError(Fmt.no_convergence(b, s))

514 __2 = _0_5 # Rb > r > Rc

515 for i in range(_TRIPS): # 47-48

516 s = (b + c) * __2

517 R, P, d = _R3(s)

518 if Rb > R > r:

519 b, Rb = s, R

520 elif Rc < R < r:

521 c, Rc = s, R

522 t = _a(b - c)

523 if t < tol: # or _a(R - r) < tol:

524 break

525 else: # t = min(t, _a(R - r))

526 raise ValueError(Fmt.no_convergence(t, tol))

527 i += C.iteration # combine iterations

528 P.set_(lat0=lat0, lon0=lon0, azi0=d.azi1, iteration=i,

529 a02=d.a12, s02=d.s12, at=d.azi2 - P.azi2, name=n)

530 return P, s

531

532 # get the perpendicular intersection of 2 geodesics,

533 # one the plumb, pseudo-rhumb line to the other

534 C = _PlumbTo(gl, lat0, lon0, tol=tol)

535 try:

536 a = _a(C.s02) # distance between centers

537 if a < r:

538 c = C.s12 # distance along pseudo-rhumb line

539 h = _copysign(r, c) # past half chord length

540 P, p = _bisect2( h, c, a, r, tol, _R3)

541 Q, q = _bisect2(-h, c, a, r, tol, _R3)

542 if _a(p - q) < max(EPS, tol):

543 Q = P

544 elif a > r:

545 raise ValueError(_too_(Fmt.distant(a)))

546 else: # tangential

547 P = Q = C

548 except Exception as x:

549 t = _LLB(C.lat2, C.lon2).toStr(form=form)

550 t = _SPACE_(x, _plumb_, _to_, Fmt.PAREN(t))

551 raise IntersectionError(t, txt=None, cause=x)

552

553 return P, Q

554

555

556 def _PlumbTo(gl, lat0, lon0, est=None, tol=_TOL):

557 # (INTERNAL) Return the I{perpendicular} intersection of

558 # a geodesic from C{(lat0, lon0)} and a geodesic (line).

559 pl = _MODS.rhumb.bases._PseudoRhumbLine(gl)

560 return pl.PlumbTo(lat0, lon0, exact=gl.geodesic,

561 est=est, tol=tol)

562

563# **) MIT License

564#

565# Copyright (C) 2016-2024 -- mrJean1 at Gmail -- All Rights Reserved.

566#

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

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

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

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

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

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

573#

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

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

576#

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

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

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

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

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

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

583# OTHER DEALINGS IN THE SOFTWARE.