Coverage for pygeodesy/geodesicw.py: 89%

209 statements  

« prev     ^ index     » next       coverage.py v7.2.2, created at 2023-12-24 11:18 -0500

1 

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

3 

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

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

6that package is installed. 

7 

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

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

10 

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

12L{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 

16from pygeodesy.basics import _copysign, _xinstanceof 

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

18from pygeodesy.datums import _earth_datum, _WGS84, _EWGS84 

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

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

21from pygeodesy.errors import IntersectionError, GeodesicError 

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

23 _to_, _too_,_under 

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

25 _kWrapped, Inverse10Tuple 

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

27from pygeodesy.lazily import _ALL_LAZY, _ALL_MODS as _MODS 

28from pygeodesy.named import callername, classname 

29from pygeodesy.namedTuples import Destination3Tuple, Distance3Tuple 

30from pygeodesy.props import Property, Property_RO, property_RO 

31from pygeodesy.streprs import Fmt, unstr 

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

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

34 

35from contextlib import contextmanager 

36# from math import fabs # from .utily 

37 

38__all__ = _ALL_LAZY.geodesicw 

39__version__ = '23.12.09' 

40 

41_plumb_ = 'plumb' 

42_TRIPS = 129 

43 

44 

45class _gWrapped(_kWrapped): 

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

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

48 ''' 

49 

50 @Property_RO # MCCABE 24 

51 def Geodesic(self): 

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

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

54 provided the latter is installed. 

55 ''' 

56 _Geodesic = self.geographiclib.Geodesic 

57 # assert Caps._STD == _Geodesic.STANDARD 

58 

59 class Geodesic(_Geodesic): 

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

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

62 ''' 

63 _datum = _WGS84 

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

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

66 

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

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

69 

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

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

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

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

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

75 ''' 

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

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

78 _Geodesic.__init__(self, *args) 

79 

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

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

82 ''' 

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

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

85 return GDict(d) 

86 

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

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

89 ''' 

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

91 

92 Area = _Geodesic.Polygon # like GeodesicExact.Area 

93 

94 @property_RO 

95 def datum(self): 

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

97 ''' 

98 return self._datum 

99 

100 @Property 

101 def debug(self): 

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

103 ''' 

104 return bool(self._debug) 

105 

106 @debug.setter # PYCHOK setter! 

107 def debug(self, debug): 

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

109 details in L{GDict} results. 

110 ''' 

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

112 

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

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

115 ''' 

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

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

118 return GDict(d) 

119 

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

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

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

123 ''' 

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

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

126 

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

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

129 ''' 

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

131 

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

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

134 ''' 

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

136 

137 @Property_RO 

138 def ellipsoid(self): 

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

140 ''' 

141 return self.datum.ellipsoid 

142 

143 @property_RO 

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

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

146 ''' 

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

148 

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

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

151 ''' 

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

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

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

155 

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

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

158 ''' 

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

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

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

162 

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

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

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

166 ''' 

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

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

169 return self._Line13(t) 

170 

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

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

173 ''' 

174 if wrap: 

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

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

177 

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

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

180 ''' 

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

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

183 return GDict(d) 

184 

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

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

187 

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

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

190 ''' 

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

192 # and .HeightIDWkarney._distances 

193 if wrap: 

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

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

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

197 return fabs(r.a12) 

198 

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

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

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

202 ''' 

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

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

205 

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

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

208 ''' 

209 if wrap: 

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

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

212 

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

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

215 ''' 

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

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

218 return self._Line13(t) 

219 

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

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

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

223 ''' 

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

225 

226 def _Line13(self, t): 

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

228 to reference point 3. 

229 ''' 

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

231 salp1=t.salp1, calp1=t.calp1) 

232 gl.a13, gl.s13 = t.a13, t.s13 

233 return gl 

234 

235# Polygon = _Geodesic.Polygon 

236 

237 # Geodesic.ArcDirect.__doc__ = _Geodesic.ArcDirect.__doc__ 

238 # Geodesic.Direct.__doc__ = _Geodesic.Direct.__doc__ 

239 # Geodesic.Inverse.__doc__ = _Geodesic.Inverse.__doc__ 

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

241 # Geodesic.Line.__doc__ = _Geodesic.Line.__doc__ 

242 return Geodesic 

243 

244 @Property_RO # MCCABE 16 

245 def GeodesicLine(self): 

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

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

248 provided the latter is installed. 

249 ''' 

250 _GeodesicLine = self.geographiclib.GeodesicLine 

251 

252 class GeodesicLine(_GeodesicLine): 

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

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

255 ''' 

256 _geodesic = None 

257 

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

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

260 

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

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

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

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

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

266 specifying the capabilities the C{GeodesicLine} 

267 instance should possess (plus optional keyword 

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

269 ''' 

270 _xinstanceof(_wrapped.Geodesic, geodesic=geodesic) 

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

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

273 self._geodesic = geodesic 

274 

275 @Property_RO 

276 def a1(self): 

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

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

279 

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

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

282 ''' 

283 try: 

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

285 except AttributeError: 

286 return NAN # see .geodesicx.gxline._GeodesicLineExact 

287 

288 equatorarc = a1 

289 

290 def Arc(self): 

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

292 ''' 

293 return self.a13 

294 

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

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

297 

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

299 (C{degrees}). 

300 

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

302 ''' 

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

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

305 return GDict(d) 

306 

307 @Property_RO 

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

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

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

311 

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

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

314 ''' 

315 try: 

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

317 except AttributeError: 

318 return NAN # see .geodesicx.gxline._GeodesicLineExact 

319 

320 equatorazimuth = azi0 

321 

322 def Distance(self): 

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

324 ''' 

325 return self.s13 

326 

327 @property_RO 

328 def geodesic(self): 

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

330 ''' 

331 return self._geodesic 

332 

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

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

335 

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

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

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

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

340 

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

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

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

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

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

346 the geodesic from the circle center to the intersection 

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

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

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

350 tangential to the circle, both intersections are the same 

351 L{GDict} instance. 

352 

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

354 intersect. 

355 

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

357 ''' 

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

359 

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

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

362 with a geodesic from the given point. 

363 

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

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

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

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

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

369 

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

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

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

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

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

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

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

377 geodesic and this line at the intersection point. The 

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

379 See method L{Position} for further details. 

380 

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

382 ''' 

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

384 

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

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

387 

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

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

390 the quantities to be returned. 

391 

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

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

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

395 except when C{a12=NAN}. 

396 ''' 

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

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

399 return GDict(d) 

400 

401 # GeodesicLine.ArcPosition.__doc__ = _GeodesicLine.ArcPosition.__doc__ 

402 # GeodesicLine.Position.__doc__ = _GeodesicLine.Position.__doc__ 

403 return GeodesicLine 

404 

405 @Property_RO 

406 def Geodesic_WGS84(self): 

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

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

409 is installed, otherwise an C{ImportError}. 

410 ''' 

411 return _EWGS84.geodesic 

412 

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

414 

415 

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

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

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

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

420 

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

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

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

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

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

426 ''' 

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

428 

429 

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

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

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

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

434 

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

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

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

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

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

440 specifying the capabilities the C{GeodesicLine} 

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

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

443 and C{GeodesicLine.ArcPosition}. 

444 ''' 

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

446 

447 

448def Geodesic_WGS84(): 

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

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

451 installed, otherwise an C{ImportError}. 

452 ''' 

453 return _wrapped.Geodesic_WGS84 

454 

455 

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

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

458 ''' 

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

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

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

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

463 ''' 

464 try: 

465 yield args 

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

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

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

469 

470_wargs = _wargs() # PYCHOK singleton 

471 

472 

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

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

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

476 r = Radius_(radius) 

477 n = _dunder_nameof(_Intersecant2)[1:] 

478 _P = gl.Position 

479 _I = gl.geodesic.Inverse 

480 _a = fabs 

481 

482 def _R3(s): 

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

484 P = _P(s) 

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

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

487 

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

489 b = c 

490 while True: 

491 b += s 

492 Rb, P, d = _R3(b) 

493 if Rb > r: 

494 break 

495 # assert Rb > r > Rc 

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

497 s = (b + c) * _0_5 

498 R, P, d = _R3(s) 

499 if Rb > R > r: 

500 b, Rb = s, R 

501 elif Rc < R < r: 

502 c, Rc = s, R 

503 t = _a(b - c) 

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

505 break 

506 else: 

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

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

509 i += C.iteration # combine iterations 

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 

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

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

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

541 pl = _MODS.rhumb.bases._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-2024 -- 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.