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 # from .karney

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

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

21 from pygeodesy.errors import IntersectionError, GeodesicError

22 from pygeodesy.interns import NN, _DOT_, _dunder_nameof, _SPACE_, \

23 _to_, _too_,_under

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

25 _kWrapped, Inverse10Tuple, _earth_datum, \

26 _ALL_LAZY, _EWGS84, _MODS, _WGS84 # PYCHOK used!

27 from pygeodesy.latlonBase import LatLonBase as _LLB, F_D

28# from pygeodesy.lazily import _ALL_LAZY, _ALL_MODS as _MODS # from .karney

29 from pygeodesy.named import callername, classname

30 from pygeodesy.namedTuples import Destination3Tuple, Distance3Tuple

31 from pygeodesy.props import Property, Property_RO, property_RO

32 from pygeodesy.streprs import Fmt, unstr

33# from pygeodesy.units import Radius_ # _MODS

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__ = '23.11.11'

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

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

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

70

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

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

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

74 B{C{a_ellipsoid}) is not specified as C{scalar}.

75 @kwarg name: Optional ellipsoid name (C{str}), ignored like B{C{f}}.

76 '''

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

78 with _wargs(self, *self.ellipsoid.a_f, name=name) as args:

79 _Geodesic.__init__(self, *args)

80

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

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

83 '''

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

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

86 return GDict(d)

87

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

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

90 '''

91 return self._GenDirectLine(lat1, lon1, azi1, True, a12, caps)

92

93 Area = _Geodesic.Polygon # like GeodesicExact.Area

94

95 @property_RO

96 def datum(self):

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

98 '''

99 return self._datum

100

101 @Property

102 def debug(self):

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

104 '''

105 return bool(self._debug)

106

107 @debug.setter # PYCHOK setter!

108 def debug(self, debug):

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

110 details in L{GDict} results.

111 '''

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

113

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

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

116 '''

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

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

119 return GDict(d)

120

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

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

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

124 '''

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

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

127

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

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

130 '''

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

132

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

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

135 '''

136 return self._GenDirectLine(lat1, lon1, azi1, False, s12, caps)

137

138 @Property_RO

139 def ellipsoid(self):

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

141 '''

142 return self.datum.ellipsoid

143

144 @property_RO

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

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

147 '''

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

149

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

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

152 '''

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

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

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

156

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

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

159 '''

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

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

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

163

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

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

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

167 '''

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

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

170 return self._Line13(t)

171

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

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

174 '''

175 if wrap:

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

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

178

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

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

181 '''

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

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

184 return GDict(d)

185

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

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

188

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

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

191 '''

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

193 # and .HeightIDWkarney._distances

194 if wrap:

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

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

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

198 return fabs(r.a12)

199

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

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

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

203 '''

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

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

206

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

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

209 '''

210 if wrap:

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

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

213

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

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

216 '''

217 with _wargs(self, lat1, lon1, lat2, lon2, caps) as args:

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

219 return self._Line13(t)

220

221 def Line(self, lat1, lon1, azi1, caps=Caps._STD_LINE):

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

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

224 '''

225 return _wrapped.GeodesicLine(self, lat1, lon1, azi1, caps=caps)

226

227 def _Line13(self, t):

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

229 to reference point 3.

230 '''

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

232 salp1=t.salp1, calp1=t.calp1)

233 gl.a13, gl.s13 = t.a13, t.s13

234 return gl

235

236# Polygon = _Geodesic.Polygon

237

238 # Geodesic.ArcDirect.__doc__ = _Geodesic.ArcDirect.__doc__

239 # Geodesic.Direct.__doc__ = _Geodesic.Direct.__doc__

240 # Geodesic.Inverse.__doc__ = _Geodesic.Inverse.__doc__

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

242 # Geodesic.Line.__doc__ = _Geodesic.Line.__doc__

243 return Geodesic

244

245 @Property_RO # MCCABE 16

246 def GeodesicLine(self):

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

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

249 provided the latter is installed.

250 '''

251 _GeodesicLine = self.geographiclib.GeodesicLine

252

253 class GeodesicLine(_GeodesicLine):

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

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

256 '''

257 _geodesic = None

258

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

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

261

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

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

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

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

266 @kwarg caps_: Optional, bit-or'ed combination of L{Caps} values

267 specifying the capabilities the C{GeodesicLine}

268 instance should possess (plus optional keyword

269 arguments C{salp1=NAN} and C{calp1=NAN}).

270 '''

271 _xinstanceof(_wrapped.Geodesic, geodesic=geodesic)

272 with _wargs(self, geodesic, lat1, lon1, azi1, **caps_) as args:

273 _GeodesicLine.__init__(self, *args, **caps_)

274 self._geodesic = geodesic

275

276 @Property_RO

277 def a1(self):

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

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

280

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

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

283 '''

284 try:

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

286 except AttributeError:

287 return NAN # see .geodesicx.gxline._GeodesicLineExact

288

289 equatorarc = a1

290

291 def Arc(self):

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

293 '''

294 return self.a13

295

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

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

298

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

300 (C{degrees}).

301

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

303 '''

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

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

306 return GDict(d)

307

308 @Property_RO

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

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

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

312

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

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

315 '''

316 try:

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

318 except AttributeError:

319 return NAN # see .geodesicx.gxline._GeodesicLineExact

320

321 equatorazimuth = azi0

322

323 def Distance(self):

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

325 '''

326 return self.s13

327

328 @property_RO

329 def geodesic(self):

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

331 '''

332 return self._geodesic

333

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

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

336

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

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

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

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

341

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

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

344 extended to 14 items by C{lon0, lat0, azi0, a02, s02, at}

345 with the circle center C{lat0}, C{lon0}, azimuth C{azi0} at,

346 distance C{a02} in C{degrees} and C{s02} in C{meter} along

347 the geodesic from the circle center to the intersection

348 C{lat2}, C{lon2} and the angle C{at} between the geodesic

349 and this line at the intersection. The I{geodesic} azimuth

350 at the intersection is C{(at + azi2)}. If this line is

351 tangential to the circle, both intersections are the same

352 L{GDict} instance.

353

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

355 intersect.

356

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

358 '''

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

360

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

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

363 with a geodesic from the given point.

364

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

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

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

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

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

370

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

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

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

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

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

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

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

378 geodesic and this line at the intersection point. The

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

380 See method L{Position} for further details.

381

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

383 '''

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

385

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

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

388

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

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

391 the quantities to be returned.

392

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

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

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

396 except when C{a12=NAN}.

397 '''

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

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

400 return GDict(d)

401

402 # GeodesicLine.ArcPosition.__doc__ = _GeodesicLine.ArcPosition.__doc__

403 # GeodesicLine.Position.__doc__ = _GeodesicLine.Position.__doc__

404 return GeodesicLine

405

406 @Property_RO

407 def Geodesic_WGS84(self):

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

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

410 is installed, otherwise an C{ImportError}.

411 '''

412 return _EWGS84.geodesic

413

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

415

416

417 def Geodesic(a_ellipsoid, f=None, name=NN):

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

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

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

421

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

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

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

425 B{C{a_ellipsoid}}) is not specified as C{meter}.

426 @kwarg name: Optional ellipsoid name (C{str}), ignored like B{C{f}}.

427 '''

428 return _wrapped.Geodesic(a_ellipsoid, f=f, name=name)

429

430

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

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

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

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

435

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

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

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

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

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

441 specifying the capabilities the C{GeodesicLine}

442 instance should possess, i.e., which quantities can

443 be returned by calls to C{GeodesicLine.Position}

444 and C{GeodesicLine.ArcPosition}.

445 '''

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

447

448

449 def Geodesic_WGS84():

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

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

452 installed, otherwise an C{ImportError}.

453 '''

454 return _wrapped.Geodesic_WGS84

455

456

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

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

459 '''

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

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

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

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

464 '''

465 try:

466 yield args

467 except (AttributeError, TypeError, ValueError) as x:

468 n = _DOT_(classname(inst), callername(up=3, underOK=True))

469 raise GeodesicError(unstr(n, *args, **kwds), cause=x)

470

471_wargs = _wargs() # PYCHOK singleton

472

473

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

475 # (INTERNAL) Return a 2-Tuple C{(P, Q)} with the intersections of

476 # a circle at C{lat0, lon0} and a geodesic line, each a C{GDict}.

477 r = _MODS.units.Radius_(radius)

478 n = _dunder_nameof(_Intersecant2)[1:]

479 _P = gl.Position

480 _I = gl.geodesic.Inverse

481 _a = fabs

482

483 def _R3(s):

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

485 P = _P(s)

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

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

488

489 def _bisect2(s, c, Rc, r, tol):

490 b = c

491 while True:

492 b += s

493 Rb, P, d = _R3(b)

494 if Rb > r:

495 break

496 assert Rb > r > Rc

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

498 s = (b + c) * _0_5

499 R, P, d = _R3(s)

500 if Rb > R > r:

501 b, Rb = s, R

502 elif Rc < R < r:

503 c, Rc = s, R

504 t = _a(b - c)

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

506 break

507 else:

508 # t = min(t, _a(R - r))

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

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

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

512 return P, s

513

514 # get the perpendicular intersection of 2

515 # geodesics, one as a pseudo-rhumb line

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

517 try:

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

519 if a < r:

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

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

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

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

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

525 Q = P

526 elif a > r:

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

528 else: # tangential

529 P = Q = C

530 except Exception as x:

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

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

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

534

535 return P, Q

536

537

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

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

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

541 pl = _MODS.rhumbBase._PseudoRhumbLine(gl)

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

543 est=est, tol=tol)

544

545# **) MIT License

546#

547# Copyright (C) 2016-2023 -- mrJean1 at Gmail -- All Rights Reserved.

548#

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

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

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

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

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

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

555#

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

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

558#

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

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

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

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

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

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

565# OTHER DEALINGS IN THE SOFTWARE.