Coverage for pygeodesy/geodesicw.py: 89%

221 statements  

« prev     ^ index     » next       coverage.py v7.2.2, created at 2024-05-15 16:36 -0400

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, _xkwds_pop2 

22from pygeodesy.fsums import Fsum, Fmt, unstr 

23from pygeodesy.internals import _dunder_nameof, _under 

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

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

26 _kWrapped, Inverse10Tuple 

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

28from pygeodesy.lazily import _ALL_LAZY, _ALL_MODS as _MODS 

29from pygeodesy.named import callername, classname 

30from pygeodesy.namedTuples import Destination3Tuple, Distance3Tuple 

31from 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 

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

35 

36from contextlib import contextmanager 

37# from math import fabs # from .utily 

38 

39__all__ = _ALL_LAZY.geodesicw 

40__version__ = '24.05.14' 

41 

42_plumb_ = 'plumb' 

43_TRIPS = 65 

44 

45 

46class _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=NN): # 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}), ignored if B{C{a_ellipsoid}) 

75 is not C{meter}. 

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

77 ''' 

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

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

80 _Geodesic.__init__(self, *args) 

81 if name: 

82 self._name = str(name) 

83 

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

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

86 ''' 

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

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

89 return GDict(d) 

90 

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

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

93 ''' 

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

95 

96 Area = _Geodesic.Polygon # like GeodesicExact.Area 

97 

98 @property_RO 

99 def datum(self): 

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

101 ''' 

102 return self._datum 

103 

104 @Property 

105 def debug(self): 

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

107 ''' 

108 return bool(self._debug) 

109 

110 @debug.setter # PYCHOK setter! 

111 def debug(self, debug): 

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

113 details in L{GDict} results. 

114 ''' 

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

116 

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

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

119 ''' 

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

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

122 return GDict(d) 

123 

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

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

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

127 ''' 

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

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

130 

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

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

133 ''' 

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

135 

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

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

138 ''' 

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

140 

141 @Property_RO 

142 def ellipsoid(self): 

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

144 ''' 

145 return self.datum.ellipsoid 

146 

147 @property_RO 

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

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

150 ''' 

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

152 

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

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

155 ''' 

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

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

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

159 

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

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

162 ''' 

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

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

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

166 

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

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

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

170 ''' 

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

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

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

174 

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

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

177 ''' 

178 if wrap: 

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

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

181 

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

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

184 ''' 

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

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

187 return GDict(d) 

188 

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

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

191 

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

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

194 ''' 

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

196 # and .HeightIDWkarney._distances 

197 if wrap: 

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

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

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

201 return fabs(r.a12) 

202 

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

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

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

206 ''' 

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

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

209 

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

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

212 ''' 

213 if wrap: 

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

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

216 

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

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

219 ''' 

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

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

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

223 

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

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

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

227 ''' 

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

229 

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

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

232 to reference point 3. 

233 ''' 

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

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

236 gl.a13, gl.s13 = t.a13, t.s13 

237 return gl 

238 

239 @property_RO 

240 def name(self): 

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

242 ''' 

243 return self._name 

244 

245# Polygon = _Geodesic.Polygon 

246 

247 # Geodesic.ArcDirect.__doc__ = _Geodesic.ArcDirect.__doc__ 

248 # Geodesic.Direct.__doc__ = _Geodesic.Direct.__doc__ 

249 # Geodesic.Inverse.__doc__ = _Geodesic.Inverse.__doc__ 

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

251 # Geodesic.Line.__doc__ = _Geodesic.Line.__doc__ 

252 return Geodesic 

253 

254 @Property_RO # MCCABE 16 

255 def GeodesicLine(self): 

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

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

258 provided the latter is installed. 

259 ''' 

260 _GeodesicLine = self.geographiclib.GeodesicLine 

261 

262 class GeodesicLine(_GeodesicLine): 

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

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

265 ''' 

266 _geodesic = None 

267 _name = NN 

268 

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

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

271 

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

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

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

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

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

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

278 capabilities the C{GeodesicLine} instance should possess, 

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

280 for I{INTERNAL} use. 

281 ''' 

282 _xinstanceof(_wrapped.Geodesic, geodesic=geodesic) 

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

284 name, caps_ = _xkwds_pop2(caps_name_, name=geodesic.name) 

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

286 if name: 

287 self._name = str(name) 

288 self._geodesic = geodesic 

289 

290 @Property_RO 

291 def a1(self): 

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

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

294 

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

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

297 ''' 

298 try: 

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

300 except AttributeError: 

301 return NAN # see .geodesicx.gxline._GeodesicLineExact 

302 

303 equatorarc = a1 

304 

305 def Arc(self): 

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

307 ''' 

308 return self.a13 

309 

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

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

312 

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

314 (C{degrees}). 

315 

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

317 ''' 

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

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

320 return GDict(d) 

321 

322 @Property_RO 

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

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

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

326 

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

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

329 ''' 

330 try: 

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

332 except AttributeError: 

333 return NAN # see .geodesicx.gxline._GeodesicLineExact 

334 

335 equatorazimuth = azi0 

336 

337 def Distance(self): 

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

339 ''' 

340 return self.s13 

341 

342 @property_RO 

343 def geodesic(self): 

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

345 ''' 

346 return self._geodesic 

347 

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

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

350 

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

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

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

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

355 

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

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

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

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

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

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

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

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

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

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

366 intersections are the same L{GDict} instance. 

367 

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

369 intersect. 

370 

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

372 ''' 

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

374 

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

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

377 with a geodesic from the given point. 

378 

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

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

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

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

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

384 

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

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

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

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

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

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

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

392 geodesic and this line at the intersection point. The 

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

394 See method L{Position} for further details. 

395 

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

397 ''' 

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

399 

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

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

402 

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

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

405 the quantities to be returned. 

406 

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

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

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

410 except when C{a12=NAN}. 

411 ''' 

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

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

414 return GDict(d) 

415 

416 # GeodesicLine.ArcPosition.__doc__ = _GeodesicLine.ArcPosition.__doc__ 

417 # GeodesicLine.Position.__doc__ = _GeodesicLine.Position.__doc__ 

418 return GeodesicLine 

419 

420 @Property_RO 

421 def Geodesic_WGS84(self): 

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

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

424 is installed, otherwise an C{ImportError}. 

425 ''' 

426 return _EWGS84.geodesic 

427 

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

429 

430 

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

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

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

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

435 

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

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

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

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

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

441 ''' 

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

443 

444 

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

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

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

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

449 

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

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

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

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

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

455 specifying the capabilities the C{GeodesicLine} 

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

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

458 and C{GeodesicLine.ArcPosition}. 

459 ''' 

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

461 

462 

463def Geodesic_WGS84(): 

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

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

466 installed, otherwise an C{ImportError}. 

467 ''' 

468 return _wrapped.Geodesic_WGS84 

469 

470 

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

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

473 ''' 

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

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

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

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

478 ''' 

479 try: 

480 yield args 

481 except Exception as x: 

482 n = _DOT_(classname(inst), callername(up=2, underOK=True)) 

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

484 

485_wargs = _wargs() # PYCHOK singleton 

486 

487 

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

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

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

491 r = Radius_(radius) 

492 n = _dunder_nameof(_Intersecant2)[1:] 

493 _P = gl.Position 

494 _I = gl.geodesic.Inverse 

495 _a = fabs 

496 

497 def _R3(s): 

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

499 P = _P(s) 

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

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

502 

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

504 _s = Fsum(c).fsumf_ 

505 for i in range(_TRIPS): 

506 b = _s(s) 

507 Rb, P, d = _R3(b) 

508 if Rb > r: 

509 break 

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

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

512 __2 = _0_5 # Rb > r > Rc 

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

514 s = (b + c) * __2 

515 R, P, d = _R3(s) 

516 if Rb > R > r: 

517 b, Rb = s, R 

518 elif Rc < R < r: 

519 c, Rc = s, R 

520 t = _a(b - c) 

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

522 break 

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

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

525 i += C.iteration # combine iterations 

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

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

528 return P, s 

529 

530 # get the perpendicular intersection of 2 geodesics, 

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

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

533 try: 

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

535 if a < r: 

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

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

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

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

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

541 Q = P 

542 elif a > r: 

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

544 else: # tangential 

545 P = Q = C 

546 except Exception as x: 

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

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

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

550 

551 return P, Q 

552 

553 

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

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

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

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

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

559 est=est, tol=tol) 

560 

561# **) MIT License 

562# 

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

564# 

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

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

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

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

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

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

571# 

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

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

574# 

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

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

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

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

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

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

581# OTHER DEALINGS IN THE SOFTWARE.