Coverage for pygeodesy/geodesicx/gx.py: 93%

522 statements  

« prev     ^ index     » next       coverage.py v7.2.2, created at 2024-06-10 14:08 -0400

1 

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

3 

4u'''A pure Python version of I{Karney}'s C++ class U{GeodesicExact 

5<https://GeographicLib.SourceForge.io/C++/doc/classGeographicLib_1_1GeodesicExact.html>}. 

6 

7Class L{GeodesicExact} follows the naming, methods and return values 

8of class C{Geodesic} from I{Karney}'s Python U{geographiclib 

9<https://GitHub.com/geographiclib/geographiclib-python>}. 

10 

11Copyright (C) U{Charles Karney<mailto:Karney@Alum.MIT.edu>} (2008-2023) 

12and licensed under the MIT/X11 License. For more information, see the 

13U{GeographicLib<https://GeographicLib.SourceForge.io>} documentation. 

14''' 

15# make sure int/int division yields float quotient 

16from __future__ import division as _; del _ # PYCHOK semicolon 

17 

18# A copy of comments from Karney's C{GeodesicExact.cpp}: 

19# 

20# This is a reformulation of the geodesic problem. The 

21# notation is as follows: 

22# - at a general point (no suffix or 1 or 2 as suffix) 

23# - phi = latitude 

24# - lambda = longitude 

25# - beta = latitude on auxiliary sphere 

26# - omega = longitude on auxiliary sphere 

27# - alpha = azimuth of great circle 

28# - sigma = arc length along great circle 

29# - s = distance 

30# - tau = scaled distance (= sigma at multiples of PI/2) 

31# - at northwards equator crossing 

32# - beta = phi = 0 

33# - omega = lambda = 0 

34# - alpha = alpha0 

35# - sigma = s = 0 

36# - a 12 suffix means a difference, e.g., s12 = s2 - s1. 

37# - s and c prefixes mean sin and cos 

38 

39from pygeodesy.basics import _copysign, _xinstanceof, _xor, unsigned0 

40from pygeodesy.constants import EPS, EPS0, EPS02, MANT_DIG, NAN, PI, _EPSqrt, \ 

41 _SQRT2_2, isnan, _0_0, _0_001, _0_01, _0_1, _0_5, \ 

42 _1_0, _N_1_0, _1_75, _2_0, _N_2_0, _2__PI, _3_0, \ 

43 _4_0, _6_0, _8_0, _16_0, _90_0, _180_0, _1000_0 

44from pygeodesy.datums import _earth_datum, _WGS84, _EWGS84 

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

46from pygeodesy.errors import GeodesicError, _xkwds_pop2 

47from pygeodesy.fmath import hypot as _hypot 

48from pygeodesy.fsums import fsumf_, fsum1f_ 

49from pygeodesy.geodesicx.gxbases import _cosSeries, _GeodesicBase, \ 

50 _sincos12, _sin1cos2, _xnC4 

51from pygeodesy.geodesicx.gxline import _GeodesicLineExact, _TINY, _update_glXs 

52from pygeodesy.interns import NN, _COMMASPACE_, _DOT_, _UNDER_ 

53from pygeodesy.karney import GDict, _around, _atan2d, Caps, _cbrt, _diff182, \ 

54 _fix90, _K_2_0, _norm2, _norm180, _polynomial, \ 

55 _signBit, _sincos2, _sincos2d, _sincos2de, _unsigned2 

56from pygeodesy.lazily import _ALL_DOCS, _ALL_MODS as _MODS 

57from pygeodesy.namedTuples import Destination3Tuple, Distance3Tuple 

58from pygeodesy.props import deprecated_Property, Property, Property_RO, property_RO 

59from pygeodesy.streprs import Fmt, pairs 

60from pygeodesy.utily import atan2d as _atan2d_reverse, _unrollon, _Wrap, wrap360 

61 

62from math import atan2, copysign, cos, degrees, fabs, radians, sqrt 

63 

64__all__ = () 

65__version__ = '24.05.31' 

66 

67_MAXIT1 = 20 

68_MAXIT2 = 10 + _MAXIT1 + MANT_DIG # MANT_DIG == C++ digits 

69 

70# increased multiplier in defn of _TOL1 from 100 to 200 to fix Inverse 

71# case 52.784459512564 0 -52.784459512563990912 179.634407464943777557 

72# which otherwise failed for Visual Studio 10 (Release and Debug) 

73_TOL0 = EPS 

74_TOL1 = _TOL0 * -200 # negative 

75_TOL2 = _EPSqrt # == sqrt(_TOL0) 

76_TOL3 = _TOL2 * _0_1 

77_TOLb = _TOL2 * _TOL0 # Check on bisection interval 

78_THR1 = _TOL2 * _1000_0 + _1_0 

79 

80_TINY3 = _TINY * _3_0 

81_TOL08 = _TOL0 * _8_0 

82_TOL016 = _TOL0 * _16_0 

83 

84 

85def _atan12(*sincos12, **sineg0): 

86 '''(INTERNAL) Return C{ang12} in C{radians}. 

87 ''' 

88 return atan2(*_sincos12(*sincos12, **sineg0)) 

89 

90 

91def _eTOL2(f): 

92 # Using the auxiliary sphere solution with dnm computed at 

93 # (bet1 + bet2) / 2, the relative error in the azimuth 

94 # consistency check is sig12^2 * abs(f) * min(1, 1-f/2) / 2. 

95 # (Error measured for 1/100 < b/a < 100 and abs(f) >= 1/1000. 

96 

97 # For a given f and sig12, the max error occurs for lines 

98 # near the pole. If the old rule for computing dnm = (dn1 

99 # + dn2)/2 is used, then the error increases by a factor of 

100 # 2.) Setting this equal to epsilon gives sig12 = etol2. 

101 

102 # Here 0.1 is a safety factor (error decreased by 100) and 

103 # max(0.001, abs(f)) stops etol2 getting too large in the 

104 # nearly spherical case. 

105 t = min(_1_0, _1_0 - f * _0_5) * max(_0_001, fabs(f)) * _0_5 

106 return _TOL3 / (sqrt(t) if t > EPS02 else EPS0) 

107 

108 

109class _PDict(GDict): 

110 '''(INTERNAL) Parameters passed around in C{._GDictInverse} and 

111 optionally returned when C{GeodesicExact.debug} is C{True}. 

112 ''' 

113 def set_sigs(self, ssig1, csig1, ssig2, csig2): 

114 '''Update the C{sig1} and C{sig2} parameters. 

115 ''' 

116 self.set_(ssig1=ssig1, csig1=csig1, sncndn1=(ssig1, csig1, self.dn1), # PYCHOK dn1 

117 ssig2=ssig2, csig2=csig2, sncndn2=(ssig2, csig2, self.dn2)) # PYCHOK dn2 

118 

119 def toGDict(self): # PYCHOK no cover 

120 '''Return as C{GDict} without attrs C{sncndn1} and C{sncndn2}. 

121 ''' 

122 def _rest(sncndn1=None, sncndn2=None, **rest): # PYCHOK sncndn* not used 

123 return GDict(rest) 

124 

125 return _rest(**self) 

126 

127 

128class GeodesicExact(_GeodesicBase): 

129 '''A pure Python version of I{Karney}'s C++ class U{GeodesicExact 

130 <https://GeographicLib.SourceForge.io/C++/doc/classGeographicLib_1_1GeodesicExact.html>}, 

131 modeled after I{Karney}'s Python class U{geodesic.Geodesic<https://GitHub.com/ 

132 geographiclib/geographiclib-python>}. 

133 ''' 

134 _datum = _WGS84 

135 _nC4 = 30 # default C4order 

136 

137 def __init__(self, a_ellipsoid=_EWGS84, f=None, C4order=None, **name_C4Order): # for backward compatibility 

138 '''New L{GeodesicExact} instance. 

139 

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

141 the equatorial radius of the ellipsoid (C{scalar}, 

142 conventionally in C{meter}), see B{C{f}}. 

143 @arg f: The flattening of the ellipsoid (C{scalar}) if B{C{a_ellipsoid}} 

144 is specified as C{scalar}. 

145 @kwarg C4order: Optional series expansion order (C{int}), see property 

146 L{C4order}, default C{30}. 

147 @kwarg name_C4Order: Optional C{B{name}=NN} (C{str}) and the DEPRECATED 

148 keyword argument C{C4Order}, use B{C{C4order}} instead. 

149 

150 @raise GeodesicError: Invalid B{C{C4order}}. 

151 ''' 

152 if name_C4Order: 

153 C4order, name = _xkwds_pop2(name_C4Order, C4Order=C4order) 

154 if name: 

155 self.name = name 

156 else: 

157 name = {} # name_C4Order 

158 

159 _earth_datum(self, a_ellipsoid, f=f, **name) 

160 if C4order: # XXX private copy, always? 

161 self.C4order = C4order 

162 

163 @Property_RO 

164 def a(self): 

165 '''Get the I{equatorial} radius, semi-axis (C{meter}). 

166 ''' 

167 return self.ellipsoid.a 

168 

169 def ArcDirect(self, lat1, lon1, azi1, a12, outmask=Caps.STANDARD): 

170 '''Solve the I{Direct} geodesic problem in terms of (spherical) arc length. 

171 

172 @arg lat1: Latitude of the first point (C{degrees}). 

173 @arg lon1: Longitude of the first point (C{degrees}). 

174 @arg azi1: Azimuth at the first point (compass C{degrees}). 

175 @arg a12: Arc length between the points (C{degrees}), can be negative. 

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

177 the quantities to be returned. 

178 

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

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

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

182 

183 @see: C++ U{GeodesicExact.ArcDirect 

184 <https://GeographicLib.SourceForge.io/C++/doc/classGeographicLib_1_1GeodesicExact.html>} 

185 and Python U{Geodesic.ArcDirect<https://GeographicLib.SourceForge.io/Python/doc/code.html>}. 

186 ''' 

187 return self._GDictDirect(lat1, lon1, azi1, True, a12, outmask) 

188 

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

190 '''Define a L{GeodesicLineExact} in terms of the I{direct} geodesic problem and as arc length. 

191 

192 @arg lat1: Latitude of the first point (C{degrees}). 

193 @arg lon1: Longitude of the first point (C{degrees}). 

194 @arg azi1: Azimuth at the first point (compass C{degrees}). 

195 @arg a12: Arc length between the points (C{degrees}), can be negative. 

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

197 the capabilities the L{GeodesicLineExact} instance 

198 should possess, i.e., which quantities can be 

199 returned by calls to L{GeodesicLineExact.Position} 

200 and L{GeodesicLineExact.ArcPosition}. 

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

202 

203 @return: A L{GeodesicLineExact} instance. 

204 

205 @note: The third point of the L{GeodesicLineExact} is set to correspond 

206 to the second point of the I{Inverse} geodesic problem. 

207 

208 @note: Latitude B{C{lat1}} should in the range C{[-90, +90]}. 

209 

210 @see: C++ U{GeodesicExact.ArcDirectLine 

211 <https://GeographicLib.SourceForge.io/C++/doc/classGeographicLib_1_1GeodesicExact.html>} and 

212 Python U{Geodesic.ArcDirectLine<https://GeographicLib.SourceForge.io/Python/doc/code.html>}. 

213 ''' 

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

215 

216 def Area(self, polyline=False, **name): 

217 '''Set up a L{GeodesicAreaExact} to compute area and 

218 perimeter of a polygon. 

219 

220 @kwarg polyline: If C{True} perimeter only, otherwise 

221 area and perimeter (C{bool}). 

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

223 

224 @return: A L{GeodesicAreaExact} instance. 

225 

226 @note: The B{C{debug}} setting is passed as C{verbose} 

227 to the returned L{GeodesicAreaExact} instance. 

228 ''' 

229 gaX = _MODS.geodesicx.GeodesicAreaExact(self, polyline=polyline, 

230 name=self._name__(name)) 

231 if self.debug: 

232 gaX.verbose = True 

233 return gaX 

234 

235 @Property_RO 

236 def b(self): 

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

238 ''' 

239 return self.ellipsoid.b 

240 

241 @Property_RO 

242 def c2x(self): 

243 '''Get the ellipsoid's I{authalic} earth radius I{squared} (C{meter} I{squared}). 

244 ''' 

245 # The Geodesic class substitutes atanh(sqrt(e2)) for asinh(sqrt(ep2)) 

246 # in the definition of _c2. The latter is more accurate for very 

247 # oblate ellipsoids (which the Geodesic class does not handle). Of 

248 # course, the area calculation in GeodesicExact is still based on a 

249 # series and only holds for moderately oblate (or prolate) ellipsoids. 

250 return self.ellipsoid.c2x 

251 

252 c2 = c2x # in this particular case 

253 

254 def C4f(self, eps): 

255 '''Evaluate the C{C4x} coefficients for B{C{eps}}. 

256 

257 @arg eps: Polynomial factor (C{float}). 

258 

259 @return: C{C4order}-Tuple of C{C4x(B{eps})} coefficients. 

260 ''' 

261 def _c4(nC4, C4x): 

262 i, x, e = 0, _1_0, eps 

263 _p = _polynomial 

264 for r in range(nC4, 0, -1): 

265 j = i + r 

266 yield _p(e, C4x, i, j) * x 

267 x *= e 

268 i = j 

269 # assert i == (nC4 * (nC4 + 1)) // 2 

270 

271 return tuple(_c4(self._nC4, self._C4x)) 

272 

273 def _C4f_k2(self, k2): # in ._GDictInverse and gxline._GeodesicLineExact._C4a 

274 '''(INTERNAL) Compute C{eps} from B{C{k2}} and invoke C{C4f}. 

275 ''' 

276 return self.C4f(k2 / fsumf_(_2_0, sqrt(k2 + _1_0) * _2_0, k2)) 

277 

278 @Property 

279 def C4order(self): 

280 '''Get the series expansion order (C{int}, 24, 27 or 30). 

281 ''' 

282 return self._nC4 

283 

284 @C4order.setter # PYCHOK .setter! 

285 def C4order(self, order): 

286 '''Set the series expansion order (C{int}, 24, 27 or 30). 

287 

288 @raise GeodesicError: Invalid B{C{order}}. 

289 ''' 

290 _xnC4(C4order=order) 

291 if self._nC4 != order: 

292 GeodesicExact._C4x._update(self) 

293 _update_glXs(self) # zap cached _GeodesicLineExact attrs _B41, _C4a 

294 self._nC4 = order 

295 

296 @deprecated_Property 

297 def C4Order(self): 

298 '''DEPRECATED, use property C{C4order}. 

299 ''' 

300 return self.C4order 

301 

302 @C4Order.setter # PYCHOK .setter! 

303 def C4Order(self, order): 

304 '''DEPRECATED, use property C{C4order}. 

305 ''' 

306 _xnC4(C4Order=order) 

307 self.C4order = order 

308 

309 @Property_RO 

310 def _C4x(self): 

311 '''Get this ellipsoid's C{C4} coefficients, I{cached} tuple. 

312 

313 @see: Property L{C4order}. 

314 ''' 

315 # see C4coeff() in GeographicLib.src.GeodesicExactC4.cpp 

316 def _C4(nC4): 

317 i, n, cs = 0, self.n, _C4coeffs(nC4) 

318 _p = _polynomial 

319 for r in range(nC4 + 1, 1, -1): # _reverange 

320 for j in range(1, r): 

321 j = j + i # (j - i - 1) order of polynomial 

322 yield _p(n, cs, i, j) / cs[j] 

323 i = j + 1 

324 # assert i == (nC4 * (nC4 + 1) * (nC4 + 5)) // 6 

325 

326 return tuple(_C4(self._nC4)) # 3rd flattening 

327 

328 @property_RO 

329 def datum(self): 

330 '''Get the datum (C{Datum}). 

331 ''' 

332 return self._datum 

333 

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

335 '''Solve the I{Direct} geodesic problem 

336 

337 @arg lat1: Latitude of the first point (C{degrees}). 

338 @arg lon1: Longitude of the first point (C{degrees}). 

339 @arg azi1: Azimuth at the first point (compass C{degrees}). 

340 @arg s12: Distance between the points (C{meter}), can be negative. 

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

342 the quantities to be returned. 

343 

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

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

346 C{lon1}, C{azi1} and distance C{s12} always included. 

347 

348 @see: C++ U{GeodesicExact.Direct 

349 <https://GeographicLib.SourceForge.io/C++/doc/classGeographicLib_1_1GeodesicExact.html>} 

350 and Python U{Geodesic.Direct<https://GeographicLib.SourceForge.io/Python/doc/code.html>}. 

351 ''' 

352 return self._GDictDirect(lat1, lon1, azi1, False, s12, outmask) 

353 

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

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

356 (final bearing) in C{degrees}. 

357 

358 @return: L{Destination3Tuple}C{(lat, lon, final)}. 

359 ''' 

360 r = self._GDictDirect(lat1, lon1, azi1, False, s12, Caps._AZIMUTH_LATITUDE_LONGITUDE) 

361 return Destination3Tuple(r.lat2, r.lon2, r.azi2) # no iteration 

362 

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

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

365 ''' 

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

367 

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

369 '''Define a L{GeodesicLineExact} in terms of the I{direct} geodesic problem and as distance. 

370 

371 @arg lat1: Latitude of the first point (C{degrees}). 

372 @arg lon1: Longitude of the first point (C{degrees}). 

373 @arg azi1: Azimuth at the first point (compass C{degrees}). 

374 @arg s12: Distance between the points (C{meter}), can be negative. 

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

376 the capabilities the L{GeodesicLineExact} instance 

377 should possess, i.e., which quantities can be 

378 returned by calls to L{GeodesicLineExact.Position}. 

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

380 

381 @return: A L{GeodesicLineExact} instance. 

382 

383 @note: The third point of the L{GeodesicLineExact} is set to correspond 

384 to the second point of the I{Inverse} geodesic problem. 

385 

386 @note: Latitude B{C{lat1}} should in the range C{[-90, +90]}. 

387 

388 @see: C++ U{GeodesicExact.DirectLine 

389 <https://GeographicLib.SourceForge.io/C++/doc/classGeographicLib_1_1GeodesicExact.html>} and 

390 Python U{Geodesic.DirectLine<https://GeographicLib.SourceForge.io/Python/doc/code.html>}. 

391 ''' 

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

393 

394 def _dn(self, sbet, cbet): # in gxline._GeodesicLineExact.__init__ 

395 '''(INTERNAL) Helper. 

396 ''' 

397 if self.f < 0: # PYCHOK no cover 

398 dn = sqrt(_1_0 - cbet**2 * self.e2) / self.f1 

399 else: 

400 dn = sqrt(_1_0 + sbet**2 * self.ep2) 

401 return dn 

402 

403 @Property_RO 

404 def e2(self): 

405 '''Get the ellipsoid's I{(1st) eccentricity squared} (C{float}), M{f * (2 - f)}. 

406 ''' 

407 return self.ellipsoid.e2 

408 

409 @Property_RO 

410 def _e2a2(self): 

411 '''(INTERNAL) Cache M{E.e2 * E.a2}. 

412 ''' 

413 return self.e2 * self.ellipsoid.a2 

414 

415 @Property_RO 

416 def _e2_f1(self): 

417 '''(INTERNAL) Cache M{E.e2 * E.f1}. 

418 ''' 

419 return self.e2 / self.f1 

420 

421 @Property_RO 

422 def _eF(self): 

423 '''(INTERNAL) Get the elliptic function, aka C{.E}. 

424 ''' 

425 return _MODS.elliptic.Elliptic(k2=-self.ep2) 

426 

427 def _eF_reset_cHe2_f1(self, x, y): 

428 '''(INTERNAL) Reset elliptic function and return M{cH * e2 / f1 * ...}. 

429 ''' 

430 self._eF_reset_k2(x) 

431 return y * self._eF.cH * self._e2_f1 

432 

433 def _eF_reset_k2(self, x): 

434 '''(INTERNAL) Reset elliptic function and return C{k2}. 

435 ''' 

436 ep2 = self.ep2 

437 k2 = x**2 * ep2 # see .gxline._GeodesicLineExact._eF 

438 self._eF.reset(k2=-k2, alpha2=-ep2) # kp2, alphap2 defaults 

439 _update_glXs(self) # zap cached/memoized _GeodesicLineExact attrs 

440 return k2 

441 

442 @Property_RO 

443 def ellipsoid(self): 

444 '''Get the ellipsoid (C{Ellipsoid}). 

445 ''' 

446 return self.datum.ellipsoid 

447 

448 @Property_RO 

449 def ep2(self): 

450 '''Get the ellipsoid's I{2nd eccentricity squared} (C{float}), M{e2 / (1 - e2)}. 

451 ''' 

452 return self.ellipsoid.e22 # == self.e2 / self.f1**2 

453 

454 e22 = ep2 # for ellipsoid compatibility 

455 

456 @Property_RO 

457 def _eTOL2(self): 

458 '''(INTERNAL) The si12 threshold for "really short". 

459 ''' 

460 return _eTOL2(self.f) 

461 

462 @Property_RO 

463 def flattening(self): 

464 '''Get the C{ellipsoid}'s I{flattening} (C{scalar}), M{(a - b) / a}, C{0} for spherical, negative for prolate. 

465 ''' 

466 return self.ellipsoid.f 

467 

468 f = flattening 

469 

470 @Property_RO 

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

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

473 ''' 

474 return self.ellipsoid.f1 

475 

476 @Property_RO 

477 def _f180(self): 

478 '''(INTERNAL) Cached/memoized. 

479 ''' 

480 return self.f * _180_0 

481 

482 def _GDictDirect(self, lat1, lon1, azi1, arcmode, s12_a12, outmask=Caps.STANDARD): 

483 '''(INTERNAL) As C{_GenDirect}, but returning a L{GDict}. 

484 

485 @return: A L{GDict} ... 

486 ''' 

487 C = outmask if arcmode else (outmask | Caps.DISTANCE_IN) 

488 glX = self.Line(lat1, lon1, azi1, C | Caps.LINE_OFF) 

489 return glX._GDictPosition(arcmode, s12_a12, outmask) 

490 

491 def _GDictInverse(self, lat1, lon1, lat2, lon2, outmask=Caps.STANDARD): # MCCABE 33, 41 vars 

492 '''(INTERNAL) As C{_GenInverse}, but returning a L{GDict}. 

493 

494 @return: A L{GDict} ... 

495 ''' 

496 Cs = Caps 

497 if self._debug: # PYCHOK no cover 

498 outmask |= Cs._DEBUG_INVERSE & self._debug 

499 outmask &= Cs._OUT_MASK # incl. _SALPs_CALPs and _DEBUG_ 

500 # compute longitude difference carefully (with _diff182): 

501 # result is in [-180, +180] but -180 is only for west-going 

502 # geodesics, +180 is for east-going and meridional geodesics 

503 lon12, lon12s = _diff182(lon1, lon2) 

504 # see C{result} from geographiclib.geodesic.Inverse 

505 if (outmask & Cs.LONG_UNROLL): # == (lon1 + lon12) + lon12s 

506 r = GDict(lon1=lon1, lon2=fsumf_(lon1, lon12, lon12s)) 

507 else: 

508 r = GDict(lon1=_norm180(lon1), lon2=_norm180(lon2)) 

509 if _K_2_0: # GeographicLib 2.0 

510 # make longitude difference positive 

511 lon12, lon_ = _unsigned2(lon12) 

512 if lon_: 

513 lon12s = -lon12s 

514 lam12 = radians(lon12) 

515 # calculate sincosd(_around(lon12 + correction)) 

516 slam12, clam12 = _sincos2de(lon12, lon12s) 

517 # supplementary longitude difference 

518 lon12s = fsumf_(_180_0, -lon12, -lon12s) 

519 else: # GeographicLib 1.52 

520 # make longitude difference positive and if very close 

521 # to being on the same half-meridian, then make it so. 

522 if lon12 < 0: # _signBit(lon12) 

523 lon_, lon12 = True, -_around(lon12) 

524 lon12s = _around(fsumf_(_180_0, -lon12, lon12s)) 

525 else: 

526 lon_, lon12 = False, _around(lon12) 

527 lon12s = _around(fsumf_(_180_0, -lon12, -lon12s)) 

528 lam12 = radians(lon12) 

529 if lon12 > _90_0: 

530 slam12, clam12 = _sincos2d(lon12s) 

531 clam12 = -clam12 

532 else: 

533 slam12, clam12 = _sincos2(lam12) 

534 # If really close to the equator, treat as on equator. 

535 lat1 = _around(_fix90(lat1)) 

536 lat2 = _around(_fix90(lat2)) 

537 r.set_(lat1=lat1, lat2=lat2) 

538 # Swap points so that point with higher (abs) latitude is 

539 # point 1. If one latitude is a NAN, then it becomes lat1. 

540 swap_ = fabs(lat1) < fabs(lat2) or isnan(lat2) 

541 if swap_: 

542 lat1, lat2 = lat2, lat1 

543 lon_ = not lon_ 

544 if _signBit(lat1): 

545 lat_ = False # note, False 

546 else: # make lat1 <= -0 

547 lat_ = True # note, True 

548 lat1, lat2 = -lat1, -lat2 

549 # Now 0 <= lon12 <= 180, -90 <= lat1 <= -0 and lat1 <= lat2 <= -lat1 

550 # and lat_, lon_, swap_ register the transformation to bring the 

551 # coordinates to this canonical form, where False means no change 

552 # made. We make these transformations so that there are few cases 

553 # to check, e.g., on verifying quadrants in atan2. In addition, 

554 # this enforces some symmetries in the results returned. 

555 

556 # Initialize for the meridian. No longitude calculation is 

557 # done in this case to let the parameter default to 0. 

558 sbet1, cbet1 = self._sinf1cos2d(lat1) 

559 sbet2, cbet2 = self._sinf1cos2d(lat2) 

560 # If cbet1 < -sbet1, then cbet2 - cbet1 is a sensitive measure 

561 # of the |bet1| - |bet2|. Alternatively (cbet1 >= -sbet1), 

562 # abs(sbet2) + sbet1 is a better measure. This logic is used 

563 # in assigning calp2 in _Lambda6. Sometimes these quantities 

564 # vanish and in that case we force bet2 = +/- bet1 exactly. An 

565 # example where is is necessary is the inverse problem 

566 # 48.522876735459 0 -48.52287673545898293 179.599720456223079643 

567 # which failed with Visual Studio 10 (Release and Debug) 

568 if cbet1 < -sbet1: 

569 if cbet2 == cbet1: 

570 sbet2 = copysign(sbet1, sbet2) 

571 elif fabs(sbet2) == -sbet1: 

572 cbet2 = cbet1 

573 

574 p = _PDict(sbet1=sbet1, cbet1=cbet1, dn1=self._dn(sbet1, cbet1), 

575 sbet2=sbet2, cbet2=cbet2, dn2=self._dn(sbet2, cbet2)) 

576 

577 _meridian = _b = True # i.e. not meridian, not b 

578 if lat1 == -90 or slam12 == 0: 

579 # Endpoints are on a single full meridian, 

580 # so the geodesic might lie on a meridian. 

581 salp1, calp1 = slam12, clam12 # head to target lon 

582 salp2, calp2 = _0_0, _1_0 # then head north 

583 # tan(bet) = tan(sig) * cos(alp) 

584 p.set_sigs(sbet1, calp1 * cbet1, sbet2, calp2 * cbet2) 

585 # sig12 = sig2 - sig1 

586 sig12 = _atan12(sbet1, p.csig1, sbet2, p.csig2, sineg0=True) # PYCHOK csig* 

587 s12x, m12x, _, \ 

588 M12, M21 = self._Length5(sig12, outmask | Cs.REDUCEDLENGTH, p) 

589 # Add the check for sig12 since zero length geodesics 

590 # might yield m12 < 0. Test case was 

591 # echo 20.001 0 20.001 0 | GeodSolve -i 

592 # In fact, we will have sig12 > PI/2 for meridional 

593 # geodesic which is not a shortest path. 

594 if m12x >= 0 or sig12 < _1_0: 

595 # Need at least 2 to handle 90 0 90 180 

596 # Prevent negative s12 or m12 from geographiclib 1.52 

597 if sig12 < _TINY3 or (sig12 < _TOL0 and (s12x < 0 or m12x < 0)): 

598 sig12 = m12x = s12x = _0_0 

599 else: 

600 _b = False # apply .b to s12x, m12x 

601 _meridian = False 

602 C = 1 

603 # else: # m12 < 0, prolate and too close to anti-podal 

604 # _meridian = True 

605 a12 = _0_0 # if _b else degrees(sig12) 

606 

607 if _meridian: 

608 _b = sbet1 == 0 and (self.f <= 0 or lon12s >= self._f180) # and sbet2 == 0 

609 if _b: # Geodesic runs along equator 

610 calp1 = calp2 = _0_0 

611 salp1 = salp2 = _1_0 

612 sig12 = lam12 / self.f1 # == omg12 

613 somg12, comg12 = _sincos2(sig12) 

614 m12x = self.b * somg12 

615 s12x = self.a * lam12 

616 M12 = M21 = comg12 

617 a12 = lon12 / self.f1 

618 C = 2 

619 else: 

620 # Now point1 and point2 belong within a hemisphere bounded by a 

621 # meridian and geodesic is neither meridional or equatorial. 

622 p.set_(slam12=slam12, clam12=clam12) 

623 # Figure a starting point for Newton's method 

624 sig12, salp1, calp1, \ 

625 salp2, calp2, dnm = self._InverseStart6(lam12, p) 

626 if sig12 is None: # use Newton's method 

627 # pre-compute the constant _Lambda6 term, once 

628 p.set_(bet12=None if cbet2 == cbet1 and fabs(sbet2) == -sbet1 else 

629 (((cbet1 + cbet2) * (cbet2 - cbet1)) if cbet1 < -sbet1 else 

630 ((sbet1 + sbet2) * (sbet1 - sbet2)))) 

631 sig12, salp1, calp1, \ 

632 salp2, calp2, domg12 = self._Newton6(salp1, calp1, p) 

633 s12x, m12x, _, M12, M21 = self._Length5(sig12, outmask, p) 

634 if (outmask & Cs.AREA): 

635 # omg12 = lam12 - domg12 

636 s, c = _sincos2(domg12) 

637 somg12, comg12 = _sincos12(s, c, slam12, clam12) 

638 C = 3 # Newton 

639 else: # from _InverseStart6: dnm, salp*, calp* 

640 C = 4 # Short lines 

641 s, c = _sincos2(sig12 / dnm) 

642 m12x = dnm**2 * s 

643 s12x = dnm * sig12 

644 M12 = M21 = c 

645 if (outmask & Cs.AREA): 

646 somg12, comg12 = _sincos2(lam12 / (self.f1 * dnm)) 

647 

648 else: # _meridian is False 

649 somg12 = comg12 = NAN 

650 

651 r.set_(a12=a12 if _b else degrees(sig12)) # in [0, 180] 

652 

653 if (outmask & Cs.DISTANCE): 

654 r.set_(s12=unsigned0(s12x if _b else (self.b * s12x))) 

655 

656 if (outmask & Cs.REDUCEDLENGTH): 

657 r.set_(m12=unsigned0(m12x if _b else (self.b * m12x))) 

658 

659 if (outmask & Cs.GEODESICSCALE): 

660 if swap_: 

661 M12, M21 = M21, M12 

662 r.set_(M12=unsigned0(M12), 

663 M21=unsigned0(M21)) 

664 

665 if (outmask & Cs.AREA): 

666 S12 = self._InverseArea(_meridian, salp1, calp1, 

667 salp2, calp2, 

668 somg12, comg12, p) 

669 if _xor(swap_, lat_, lon_): 

670 S12 = -S12 

671 r.set_(S12=unsigned0(S12)) 

672 

673 if (outmask & (Cs.AZIMUTH | Cs._SALPs_CALPs)): 

674 if swap_: 

675 salp1, salp2 = salp2, salp1 

676 calp1, calp2 = calp2, calp1 

677 if _xor(swap_, lon_): 

678 salp1, salp2 = -salp1, -salp2 

679 if _xor(swap_, lat_): 

680 calp1, calp2 = -calp1, -calp2 

681 

682 if (outmask & Cs.AZIMUTH): 

683 r.set_(azi1=_atan2d(salp1, calp1), 

684 azi2=_atan2d_reverse(salp2, calp2, reverse=outmask & Cs.REVERSE2)) 

685 if (outmask & Cs._SALPs_CALPs): 

686 r.set_(salp1=salp1, calp1=calp1, 

687 salp2=salp2, calp2=calp2) 

688 

689 if (outmask & Cs._DEBUG_INVERSE): # PYCHOK no cover 

690 E, eF = self.ellipsoid, self._eF 

691 p.set_(C=C, a=self.a, f=self.f, f1=self.f1, 

692 e=E.e, e2=self.e2, ep2=self.ep2, 

693 c2=E.c2, c2x=self.c2x, 

694 eFcD=eF.cD, eFcE=eF.cE, eFcH=eF.cH, 

695 eFk2=eF.k2, eFa2=eF.alpha2) 

696 p.update(r) # r overrides p 

697 r = p.toGDict() 

698 return self._iter2tion(r, **p) 

699 

700 def _GenDirect(self, lat1, lon1, azi1, arcmode, s12_a12, outmask=Caps.STANDARD): 

701 '''(INTERNAL) The general I{Inverse} geodesic calculation. 

702 

703 @return: L{Direct9Tuple}C{(a12, lat2, lon2, azi2, 

704 s12, m12, M12, M21, S12)}. 

705 ''' 

706 r = self._GDictDirect(lat1, lon1, azi1, arcmode, s12_a12, outmask) 

707 return r.toDirect9Tuple() 

708 

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

710 '''(INTERNAL) Helper for C{ArcDirectLine} and C{DirectLine}. 

711 

712 @return: A L{GeodesicLineExact} instance. 

713 ''' 

714 azi1 = _norm180(azi1) 

715 # guard against underflow in salp0. Also -0 is converted to +0. 

716 s, c = _sincos2d(_around(azi1)) 

717 C = caps if arcmode else (caps | Caps.DISTANCE_IN) 

718 return _GeodesicLineExact(self, lat1, lon1, azi1, C, 

719 self._debug, s, c, **name)._GenSet(arcmode, s12_a12) 

720 

721 def _GenInverse(self, lat1, lon1, lat2, lon2, outmask=Caps.STANDARD): 

722 '''(INTERNAL) The general I{Inverse} geodesic calculation. 

723 

724 @return: L{Inverse10Tuple}C{(a12, s12, salp1, calp1, salp2, calp2, 

725 m12, M12, M21, S12)}. 

726 ''' 

727 r = self._GDictInverse(lat1, lon1, lat2, lon2, outmask | Caps._SALPs_CALPs) 

728 return r.toInverse10Tuple() 

729 

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

731 '''(INTERNAL) Short-cut version, see .base.ellipsoidalDI.intersecant2. 

732 ''' 

733 if wrap: 

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

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

736 

737 def Inverse(self, lat1, lon1, lat2, lon2, outmask=Caps.STANDARD): 

738 '''Perform the I{Inverse} geodesic calculation. 

739 

740 @arg lat1: Latitude of the first point (C{degrees}). 

741 @arg lon1: Longitude of the first point (C{degrees}). 

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

743 @arg lon2: Longitude of the second point (C{degrees}). 

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

745 the quantities to be returned. 

746 

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

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

749 C{lon1}, C{azi1} and distance C{s12} always included. 

750 

751 @note: The third point of the L{GeodesicLineExact} is set to correspond 

752 to the second point of the I{Inverse} geodesic problem. 

753 

754 @note: Both B{C{lat1}} and B{C{lat2}} should in the range C{[-90, +90]}. 

755 

756 @see: C++ U{GeodesicExact.InverseLine 

757 <https://GeographicLib.SourceForge.io/C++/doc/classGeographicLib_1_1GeodesicExact.html>} and 

758 Python U{Geodesic.InverseLine<https://GeographicLib.SourceForge.io/Python/doc/code.html>}. 

759 ''' 

760 return self._GDictInverse(lat1, lon1, lat2, lon2, outmask) 

761 

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

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

764 

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

766 B{C{lat2}} and B{C{lon2}} (C{bool}). 

767 ''' 

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

769 # and .HeightIDWkarney._distances 

770 if wrap: 

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

772 return fabs(self._GDictInverse(lat1, lon1, lat2, lon2, Caps._ANGLE_ONLY).a12) 

773 

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

775 '''Return the distance in C{meter} and the forward and 

776 reverse azimuths (initial and final bearing) in C{degrees}. 

777 

778 @return: L{Distance3Tuple}C{(distance, initial, final)}. 

779 ''' 

780 r = self._GDictInverse(lat1, lon1, lat2, lon2, Caps.AZIMUTH_DISTANCE) 

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

782 iteration=r.iteration) 

783 

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

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

786 ''' 

787 if wrap: 

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

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

790 

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

792 '''Define a L{GeodesicLineExact} in terms of the I{Inverse} geodesic problem. 

793 

794 @arg lat1: Latitude of the first point (C{degrees}). 

795 @arg lon1: Longitude of the first point (C{degrees}). 

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

797 @arg lon2: Longitude of the second point (C{degrees}). 

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

799 the capabilities the L{GeodesicLineExact} instance 

800 should possess, i.e., which quantities can be 

801 returned by calls to L{GeodesicLineExact.Position} 

802 and L{GeodesicLineExact.ArcPosition}. 

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

804 

805 @return: A L{GeodesicLineExact} instance. 

806 

807 @note: The third point of the L{GeodesicLineExact} is set to correspond 

808 to the second point of the I{Inverse} geodesic problem. 

809 

810 @note: Both B{C{lat1}} and B{C{lat2}} should in the range C{[-90, +90]}. 

811 

812 @see: C++ U{GeodesicExact.InverseLine 

813 <https://GeographicLib.SourceForge.io/C++/doc/classGeographicLib_1_1GeodesicExact.html>} and 

814 Python U{Geodesic.InverseLine<https://GeographicLib.SourceForge.io/Python/doc/code.html>}. 

815 ''' 

816 Cs = Caps 

817 r = self._GDictInverse(lat1, lon1, lat2, lon2, Cs._SALPs_CALPs) # No need for AZIMUTH 

818 C = (caps | Cs.DISTANCE) if (caps & (Cs.DISTANCE_IN & Cs._OUT_MASK)) else caps 

819 azi1 = _atan2d(r.salp1, r.calp1) 

820 return _GeodesicLineExact(self, lat1, lon1, azi1, C, # ensure a12 is distance 

821 self._debug, r.salp1, r.calp1, **name)._GenSet(True, r.a12) 

822 

823 def _InverseArea(self, _meridian, salp1, calp1, # PYCHOK 9 args 

824 salp2, calp2, 

825 somg12, comg12, p): 

826 '''(INTERNAL) Split off from C{_GDictInverse} to reduce complexity/length. 

827 

828 @return: Area C{S12}. 

829 ''' 

830 # from _Lambda6: sin(alp1) * cos(bet1) = sin(alp0), calp0 > 0 

831 salp0, calp0 = _sin1cos2(salp1, calp1, p.sbet1, p.cbet1) 

832 A4 = calp0 * salp0 

833 if A4: 

834 # from _Lambda6: tan(bet) = tan(sig) * cos(alp) 

835 k2 = calp0**2 * self.ep2 

836 C4a = self._C4f_k2(k2) 

837 B41 = _cosSeries(C4a, *_norm2(p.sbet1, calp1 * p.cbet1)) 

838 B42 = _cosSeries(C4a, *_norm2(p.sbet2, calp2 * p.cbet2)) 

839 # multiplier = a^2 * e^2 * cos(alpha0) * sin(alpha0) 

840 A4 *= self._e2a2 

841 S12 = A4 * (B42 - B41) 

842 else: # avoid problems with indeterminate sig1, sig2 on equator 

843 A4 = B41 = B42 = k2 = S12 = _0_0 

844 

845 if (_meridian and # omg12 < 3/4 * PI 

846 comg12 > -_SQRT2_2 and # lon diff not too big 

847 (p.sbet2 - p.sbet1) < _1_75): # lat diff not too big 

848 # use tan(Gamma/2) = tan(omg12/2) * 

849 # (tan(bet1/2) + tan(bet2/2)) / 

850 # (tan(bet1/2) * tan(bet2/2) + 1)) 

851 # with tan(x/2) = sin(x) / (1 + cos(x)) 

852 dbet1 = p.cbet1 + _1_0 

853 dbet2 = p.cbet2 + _1_0 

854 domg12 = comg12 + _1_0 

855 salp12 = (p.sbet1 * dbet2 + dbet1 * p.sbet2) * somg12 

856 calp12 = (p.sbet1 * p.sbet2 + dbet1 * dbet2) * domg12 

857 alp12 = atan2(salp12, calp12) * _2_0 

858 else: 

859 # alp12 = alp2 - alp1, used in atan2, no need to normalize 

860 salp12, calp12 = _sincos12(salp1, calp1, salp2, calp2) 

861 # The right thing appears to happen if alp1 = +/-180 and 

862 # alp2 = 0, viz salp12 = -0 and alp12 = -180. However, 

863 # this depends on the sign being attached to 0 correctly. 

864 # Following ensures the correct behavior. 

865 if salp12 == 0 and calp12 < 0: 

866 alp12 = _copysign(PI, calp1) 

867 else: 

868 alp12 = atan2(salp12, calp12) 

869 

870 p.set_(alp12=alp12, A4=A4, B41=B41, B42=B42, k2=k2) 

871 return S12 + self.c2x * alp12 

872 

873 def _InverseStart6(self, lam12, p): 

874 '''(INTERNAL) Return a starting point for Newton's method in 

875 C{salp1} and C{calp1} indicated by C{sig12=None}. If 

876 Newton's method doesn't need to be used, return also 

877 C{salp2}, C{calp2}, C{dnm} and C{sig12} non-C{None}. 

878 

879 @return: 6-Tuple C{(sig12, salp1, calp1, salp2, calp2, dnm)} 

880 and C{p.set_sigs} updated for Newton, C{sig12=None}. 

881 ''' 

882 sig12 = None # use Newton 

883 salp1 = calp1 = salp2 = calp2 = dnm = NAN 

884 

885 # bet12 = bet2 - bet1 in [0, PI) 

886 sbet12, cbet12 = _sincos12(p.sbet1, p.cbet1, p.sbet2, p.cbet2) 

887 shortline = cbet12 >= 0 and sbet12 < _0_5 and (p.cbet2 * lam12) < _0_5 

888 if shortline: 

889 # sin((bet1 + bet2)/2)^2 = (sbet1 + sbet2)^2 / ( 

890 # (cbet1 + cbet2)^2 + (sbet1 + sbet2)^2) 

891 t = (p.sbet1 + p.sbet2)**2 

892 s = t / ((p.cbet1 + p.cbet2)**2 + t) 

893 dnm = sqrt(_1_0 + self.ep2 * s) 

894 somg12, comg12 = _sincos2(lam12 / (self.f1 * dnm)) 

895 else: 

896 somg12, comg12 = p.slam12, p.clam12 

897 

898 # bet12a = bet2 + bet1 in (-PI, 0], note -sbet1 

899 sbet12a, cbet12a = _sincos12(-p.sbet1, p.cbet1, p.sbet2, p.cbet2) 

900 

901 c = fabs(comg12) + _1_0 # == (1 - comg12) if comg12 < 0 

902 s = somg12**2 / c 

903 t = p.sbet1 * p.cbet2 * s 

904 salp1 = p.cbet2 * somg12 

905 calp1 = (sbet12a - t) if comg12 < 0 else (sbet12 + t) 

906 

907 ssig12 = _hypot(salp1, calp1) 

908 csig12 = p.sbet1 * p.sbet2 + p.cbet1 * p.cbet2 * comg12 

909 

910 if shortline and ssig12 < self._eTOL2: # really short lines 

911 t = c if comg12 < 0 else s 

912 salp2, calp2 = _norm2(somg12 * p.cbet1, 

913 sbet12 - p.cbet1 * p.sbet2 * t) 

914 sig12 = atan2(ssig12, csig12) # do not use Newton 

915 

916 elif (self._n_0_1 or # Skip astroid calc if too eccentric 

917 csig12 >= 0 or ssig12 >= (p.cbet1**2 * self._n6PI)): 

918 pass # nothing to do, 0th order spherical approximation OK 

919 

920 else: 

921 # Scale lam12 and bet2 to x, y coordinate system where antipodal 

922 # point is at origin and singular point is at y = 0, x = -1 

923 lam12x = atan2(-p.slam12, -p.clam12) # lam12 - PI 

924 f = self.f 

925 if f < 0: # PYCHOK no cover 

926 # ssig1=sbet1, csig1=-cbet1, ssig2=sbet2, csig2=cbet2 

927 p.set_sigs(p.sbet1, -p.cbet1, p.sbet2, p.cbet2) 

928 # if lon12 = 180, this repeats a calculation made in Inverse 

929 _, m12b, m0, _, _ = self._Length5(atan2(sbet12a, cbet12a) + PI, 

930 Caps.REDUCEDLENGTH, p) 

931 t = p.cbet1 * PI # x = dlat, y = dlon 

932 x = m12b / (t * p.cbet2 * m0) - _1_0 

933 sca = (sbet12a / (x * p.cbet1)) if x < -_0_01 else (-f * t) 

934 y = lam12x / sca 

935 else: # f >= 0, however f == 0 does not get here 

936 sca = self._eF_reset_cHe2_f1(p.sbet1, p.cbet1 * _2_0) 

937 x = lam12x / sca # dlon 

938 y = sbet12a / (sca * p.cbet1) # dlat 

939 

940 if y > _TOL1 and x > -_THR1: # strip near cut 

941 if f < 0: # PYCHOK no cover 

942 calp1 = max( _0_0, x) if x > _TOL1 else max(_N_1_0, x) 

943 salp1 = sqrt(_1_0 - calp1**2) 

944 else: 

945 salp1 = min( _1_0, -x) 

946 calp1 = -sqrt(_1_0 - salp1**2) 

947 else: 

948 # Estimate alp1, by solving the astroid problem. 

949 # 

950 # Could estimate alpha1 = theta + PI/2, directly, i.e., 

951 # calp1 = y/k; salp1 = -x/(1+k); for _f >= 0 

952 # calp1 = x/(1+k); salp1 = -y/k; for _f < 0 (need to check) 

953 # 

954 # However, it's better to estimate omg12 from astroid and use 

955 # spherical formula to compute alp1. This reduces the mean 

956 # number of Newton iterations for astroid cases from 2.24 

957 # (min 0, max 6) to 2.12 (min 0, max 5). The changes in the 

958 # number of iterations are as follows: 

959 # 

960 # change percent 

961 # 1 5 

962 # 0 78 

963 # -1 16 

964 # -2 0.6 

965 # -3 0.04 

966 # -4 0.002 

967 # 

968 # The histogram of iterations is (m = number of iterations 

969 # estimating alp1 directly, n = number of iterations 

970 # estimating via omg12, total number of trials = 148605): 

971 # 

972 # iter m n 

973 # 0 148 186 

974 # 1 13046 13845 

975 # 2 93315 102225 

976 # 3 36189 32341 

977 # 4 5396 7 

978 # 5 455 1 

979 # 6 56 0 

980 # 

981 # omg12 is near PI, estimate work with omg12a = PI - omg12 

982 k = _Astroid(x, y) 

983 sca *= (y * (k + _1_0) / k) if f < 0 else \ 

984 (x * k / (k + _1_0)) 

985 s, c = _sincos2(-sca) # omg12a 

986 # update spherical estimate of alp1 using omg12 instead of lam12 

987 salp1 = p.cbet2 * s 

988 calp1 = sbet12a - s * salp1 * p.sbet1 / (c + _1_0) # c = -c 

989 

990 # sanity check on starting guess. Backwards check allows NaN through. 

991 salp1, calp1 = _norm2(salp1, calp1) if salp1 > 0 else (_1_0, _0_0) 

992 

993 return sig12, salp1, calp1, salp2, calp2, dnm 

994 

995 def _Lambda6(self, salp1, calp1, diffp, p): 

996 '''(INTERNAL) Helper. 

997 

998 @return: 6-Tuple C{(lam12, sig12, salp2, calp2, domg12, 

999 dlam12} and C{p.set_sigs} updated. 

1000 ''' 

1001 eF = self._eF 

1002 f1 = self.f1 

1003 

1004 if p.sbet1 == calp1 == 0: # PYCHOK no cover 

1005 # Break degeneracy of equatorial line 

1006 calp1 = -_TINY 

1007 

1008 # sin(alp1) * cos(bet1) = sin(alp0), # calp0 > 0 

1009 salp0, calp0 = _sin1cos2(salp1, calp1, p.sbet1, p.cbet1) 

1010 # tan(bet1) = tan(sig1) * cos(alp1) 

1011 # tan(omg1) = sin(alp0) * tan(sig1) 

1012 # = sin(bet1) * tan(alp1) 

1013 somg1 = salp0 * p.sbet1 

1014 comg1 = calp1 * p.cbet1 

1015 ssig1, csig1 = _norm2(p.sbet1, comg1) 

1016 # Without normalization we have schi1 = somg1 

1017 cchi1 = f1 * p.dn1 * comg1 

1018 

1019 # Enforce symmetries in the case abs(bet2) = -bet1. 

1020 # Need to be careful about this case, since this can 

1021 # yield singularities in the Newton iteration. 

1022 # sin(alp2) * cos(bet2) = sin(alp0) 

1023 salp2 = (salp0 / p.cbet2) if p.cbet2 != p.cbet1 else salp1 

1024 # calp2 = sqrt(1 - sq(salp2)) 

1025 # = sqrt(sq(calp0) - sq(sbet2)) / cbet2 

1026 # and subst for calp0 and rearrange to give (choose 

1027 # positive sqrt to give alp2 in [0, PI/2]). 

1028 calp2 = fabs(calp1) if p.bet12 is None else ( 

1029 sqrt((calp1 * p.cbet1)**2 + p.bet12) / p.cbet2) 

1030 # tan(bet2) = tan(sig2) * cos(alp2) 

1031 # tan(omg2) = sin(alp0) * tan(sig2). 

1032 somg2 = salp0 * p.sbet2 

1033 comg2 = calp2 * p.cbet2 

1034 ssig2, csig2 = _norm2(p.sbet2, comg2) 

1035 # without normalization we have schi2 = somg2 

1036 cchi2 = f1 * p.dn2 * comg2 

1037 

1038 # omg12 = omg2 - omg1, limit to [0, PI] 

1039 somg12, comg12 = _sincos12(somg1, comg1, somg2, comg2, sineg0=True) 

1040 # chi12 = chi2 - chi1, limit to [0, PI] 

1041 schi12, cchi12 = _sincos12(somg1, cchi1, somg2, cchi2, sineg0=True) 

1042 

1043 p.set_sigs(ssig1, csig1, ssig2, csig2) 

1044 # sig12 = sig2 - sig1, limit to [0, PI] 

1045 sig12 = _atan12(ssig1, csig1, ssig2, csig2, sineg0=True) 

1046 

1047 eta12 = self._eF_reset_cHe2_f1(calp0, salp0) * _2__PI # then ... 

1048 eta12 *= fsum1f_(eF.deltaH(*p.sncndn2), 

1049 -eF.deltaH(*p.sncndn1), sig12) 

1050 # eta = chi12 - lam12 

1051 lam12 = _atan12(p.slam12, p.clam12, schi12, cchi12) - eta12 

1052 # domg12 = chi12 - omg12 - deta12 

1053 domg12 = _atan12( somg12, comg12, schi12, cchi12) - eta12 

1054 

1055 dlam12 = NAN # dv > 0 in ._Newton6 

1056 if diffp: 

1057 d = calp2 * p.cbet2 

1058 if d: 

1059 _, dlam12, _, _, _ = self._Length5(sig12, Caps.REDUCEDLENGTH, p) 

1060 dlam12 *= f1 / d 

1061 elif p.sbet1: 

1062 dlam12 = -f1 * p.dn1 * _2_0 / p.sbet1 

1063 

1064 # p.set_(deta12=-eta12, lam12=lam12) 

1065 return lam12, sig12, salp2, calp2, domg12, dlam12 

1066 

1067 def _Length5(self, sig12, outmask, p): 

1068 '''(INTERNAL) Return M{m12b = (reduced length) / self.b} and 

1069 calculate M{s12b = distance / self.b} and M{m0}, the 

1070 coefficient of secular term in expression for reduced 

1071 length and the geodesic scales C{M12} and C{M21}. 

1072 

1073 @return: 5-Tuple C{(s12b, m12b, m0, M12, M21)}. 

1074 ''' 

1075 s12b = m12b = m0 = M12 = M21 = NAN 

1076 

1077 Cs = Caps 

1078 eF = self._eF 

1079 

1080 # outmask &= Cs._OUT_MASK 

1081 if (outmask & Cs.DISTANCE): 

1082 # Missing a factor of self.b 

1083 s12b = eF.cE * _2__PI * fsum1f_(eF.deltaE(*p.sncndn2), 

1084 -eF.deltaE(*p.sncndn1), sig12) 

1085 

1086 if (outmask & Cs._REDUCEDLENGTH_GEODESICSCALE): 

1087 m0x = -eF.k2 * eF.cD * _2__PI 

1088 J12 = -m0x * fsum1f_(eF.deltaD(*p.sncndn2), 

1089 -eF.deltaD(*p.sncndn1), sig12) 

1090 if (outmask & Cs.REDUCEDLENGTH): 

1091 m0 = m0x 

1092 # Missing a factor of self.b. Add parens around 

1093 # (csig1 * ssig2) and (ssig1 * csig2) to ensure 

1094 # accurate cancellation for coincident points. 

1095 m12b = fsum1f_(p.dn2 * (p.csig1 * p.ssig2), 

1096 -p.dn1 * (p.ssig1 * p.csig2), 

1097 J12 * (p.csig1 * p.csig2)) 

1098 if (outmask & Cs.GEODESICSCALE): 

1099 M12 = M21 = p.ssig1 * p.ssig2 + \ 

1100 p.csig1 * p.csig2 

1101 t = (p.cbet1 - p.cbet2) * self.ep2 * \ 

1102 (p.cbet1 + p.cbet2) / (p.dn1 + p.dn2) 

1103 M12 += (p.ssig2 * t + p.csig2 * J12) * p.ssig1 / p.dn1 

1104 M21 -= (p.ssig1 * t + p.csig1 * J12) * p.ssig2 / p.dn2 

1105 

1106 return s12b, m12b, m0, M12, M21 

1107 

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

1109 '''Set up a L{GeodesicLineExact} to compute several points 

1110 on a single geodesic. 

1111 

1112 @arg lat1: Latitude of the first point (C{degrees}). 

1113 @arg lon1: Longitude of the first point (C{degrees}). 

1114 @arg azi1: Azimuth at the first point (compass C{degrees}). 

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

1116 the capabilities the L{GeodesicLineExact} instance 

1117 should possess, i.e., which quantities can be 

1118 returnedby calls to L{GeodesicLineExact.Position} 

1119 and L{GeodesicLineExact.ArcPosition}. 

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

1121 

1122 @return: A L{GeodesicLineExact} instance. 

1123 

1124 @note: If the point is at a pole, the azimuth is defined by keeping 

1125 B{C{lon1}} fixed, writing C{B{lat1} = ±(90 − ε)}, and taking 

1126 the limit C{ε → 0+}. 

1127 

1128 @see: C++ U{GeodesicExact.Line 

1129 <https://GeographicLib.SourceForge.io/C++/doc/classGeographicLib_1_1GeodesicExact.html>} 

1130 and Python U{Geodesic.Line<https://GeographicLib.SourceForge.io/Python/doc/code.html>}. 

1131 ''' 

1132 return _GeodesicLineExact(self, lat1, lon1, azi1, caps, self._debug, **name) 

1133 

1134 @Property_RO 

1135 def n(self): 

1136 '''Get the C{ellipsoid}'s I{3rd flattening} (C{scalar}), M{f / (2 - f) == (a - b) / (a + b)}. 

1137 ''' 

1138 return self.ellipsoid.n 

1139 

1140 @Property_RO 

1141 def _n_0_1(self): 

1142 '''(INTERNAL) Cached once. 

1143 ''' 

1144 return fabs(self.n) > _0_1 

1145 

1146 @Property_RO 

1147 def _n6PI(self): 

1148 '''(INTERNAL) Cached once. 

1149 ''' 

1150 return fabs(self.n) * _6_0 * PI 

1151 

1152 def _Newton6(self, salp1, calp1, p): 

1153 '''(INTERNAL) Split off from C{_GDictInverse} to reduce complexity/length. 

1154 

1155 @return: 6-Tuple C{(sig12, salp1, calp1, salp2, calp2, domg12)} 

1156 and C{p.iter} and C{p.trip} updated. 

1157 ''' 

1158 _abs = fabs 

1159 # This is a straightforward solution of f(alp1) = lambda12(alp1) - 

1160 # lam12 = 0 with one wrinkle. f(alp) has exactly one root in the 

1161 # interval (0, PI) and its derivative is positive at the root. 

1162 # Thus f(alp) is positive for alp > alp1 and negative for alp < alp1. 

1163 # During the course of the iteration, a range (alp1a, alp1b) is 

1164 # maintained which brackets the root and with each evaluation of 

1165 # f(alp) the range is shrunk, if possible. Newton's method is 

1166 # restarted whenever the derivative of f is negative (because the 

1167 # new value of alp1 is then further from the solution) or if the 

1168 # new estimate of alp1 lies outside (0,PI); in this case, the new 

1169 # starting guess is taken to be (alp1a + alp1b) / 2. 

1170 salp1a = salp1b = _TINY 

1171 calp1a, calp1b = _1_0, _N_1_0 

1172 MAXIT1, TOL0 = _MAXIT1, _TOL0 

1173 HALF, TOLb = _0_5, _TOLb 

1174 tripb, TOLv = False, TOL0 

1175 for i in range(_MAXIT2): 

1176 # 1/4 meridian = 10e6 meter and random input, 

1177 # estimated max error in nm (nano meter, by 

1178 # checking Inverse problem by Direct). 

1179 # 

1180 # max iterations 

1181 # log2(b/a) error mean sd 

1182 # -7 387 5.33 3.68 

1183 # -6 345 5.19 3.43 

1184 # -5 269 5.00 3.05 

1185 # -4 210 4.76 2.44 

1186 # -3 115 4.55 1.87 

1187 # -2 69 4.35 1.38 

1188 # -1 36 4.05 1.03 

1189 # 0 15 0.01 0.13 

1190 # 1 25 5.10 1.53 

1191 # 2 96 5.61 2.09 

1192 # 3 318 6.02 2.74 

1193 # 4 985 6.24 3.22 

1194 # 5 2352 6.32 3.44 

1195 # 6 6008 6.30 3.45 

1196 # 7 19024 6.19 3.30 

1197 v, sig12, salp2, calp2, \ 

1198 domg12, dv = self._Lambda6(salp1, calp1, i < MAXIT1, p) 

1199 

1200 # 2 * _TOL0 is approximately 1 ulp [0, PI] 

1201 # reversed test to allow escape with NaNs 

1202 if tripb or _abs(v) < TOLv: 

1203 break 

1204 # update bracketing values 

1205 if v > 0 and (i > MAXIT1 or (calp1 / salp1) > (calp1b / salp1b)): 

1206 salp1b, calp1b = salp1, calp1 

1207 elif v < 0 and (i > MAXIT1 or (calp1 / salp1) < (calp1a / salp1a)): 

1208 salp1a, calp1a = salp1, calp1 

1209 

1210 if i < MAXIT1 and dv > 0: 

1211 dalp1 = -v / dv 

1212 if _abs(dalp1) < PI: 

1213 s, c = _sincos2(dalp1) 

1214 # nalp1 = alp1 + dalp1 

1215 s, c = _sincos12(-s, c, salp1, calp1) 

1216 if s > 0: 

1217 salp1, calp1 = _norm2(s, c) 

1218 # in some regimes we don't get quadratic convergence 

1219 # because slope -> 0. So use convergence conditions 

1220 # based on epsilon instead of sqrt(epsilon) 

1221 TOLv = TOL0 if _abs(v) > _TOL016 else _TOL08 

1222 continue 

1223 TOLv = TOL0 

1224 # Either dv was not positive or updated value was outside 

1225 # legal range. Use the midpoint of the bracket as the next 

1226 # estimate. This mechanism is not needed for the WGS84 

1227 # ellipsoid, but it does catch problems with more eccentric 

1228 # ellipsoids. Its efficacy is such for the WGS84 test set 

1229 # with the starting guess set to alp1 = 90 deg: the WGS84 

1230 # test set: mean = 5.21, stdev = 3.93, max = 24 and WGS84 

1231 # with random input: mean = 4.74, stdev = 0.99 

1232 salp1, calp1 = _norm2((salp1a + salp1b) * HALF, 

1233 (calp1a + calp1b) * HALF) 

1234 tripb = fsum1f_(calp1a, -calp1, _abs(salp1a - salp1)) < TOLb or \ 

1235 fsum1f_(calp1b, -calp1, _abs(salp1b - salp1)) < TOLb 

1236 else: 

1237 raise GeodesicError(Fmt.no_convergence(v, TOLv), txt=repr(self)) # self.toRepr() 

1238 

1239 p.set_(iter=i, trip=tripb) # like .geodsolve._GDictInvoke: iter NOT iteration! 

1240 return sig12, salp1, calp1, salp2, calp2, domg12 

1241 

1242 Polygon = Area # for C{geographiclib} compatibility 

1243 

1244 def _sinf1cos2d(self, lat): 

1245 '''(INTERNAL) Helper, also for C{_G_GeodesicLineExact}. 

1246 ''' 

1247 sbet, cbet = _sincos2d(lat) 

1248 # ensure cbet1 = +epsilon at poles; doing the fix on beta means 

1249 # that sig12 will be <= 2*tiny for two points at the same pole 

1250 sbet, cbet = _norm2(sbet * self.f1, cbet) 

1251 return sbet, (cbet if fabs(cbet) > _TINY else _TINY) 

1252 

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

1254 '''Return this C{GeodesicExact} as string. 

1255 

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

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

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

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

1260 

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

1262 ''' 

1263 d = dict(ellipsoid=self.ellipsoid, C4order=self.C4order) 

1264 return sep.join(pairs(d, prec=prec)) 

1265 

1266 

1267class GeodesicLineExact(_GeodesicLineExact): 

1268 '''A pure Python version of I{Karney}'s C++ class U{GeodesicLineExact 

1269 <https://GeographicLib.SourceForge.io/C++/doc/classGeographicLib_1_1GeodesicLineExact.html>}, 

1270 modeled after I{Karney}'s Python class U{geodesicline.GeodesicLine<https://GitHub.com/ 

1271 geographiclib/geographiclib-python>}. 

1272 ''' 

1273 

1274 def __init__(self, geodesic, lat1, lon1, azi1, caps=Caps.STANDARD, **name): 

1275 '''New L{GeodesicLineExact} instance, allowing points to be found along 

1276 a geodesic starting at C{(B{lat1}, B{lon1})} with azimuth B{C{azi1}}. 

1277 

1278 @arg geodesic: The geodesic to use (L{GeodesicExact}). 

1279 @arg lat1: Latitude of the first point (C{degrees}). 

1280 @arg lon1: Longitude of the first point (C{degrees}). 

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

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

1283 the capabilities the L{GeodesicLineExact} instance 

1284 should possess, i.e., which quantities can be 

1285 returned by calls to L{GeodesicLineExact.Position} 

1286 and L{GeodesicLineExact.ArcPosition}. 

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

1288 

1289 @raise TypeError: Invalid B{C{geodesic}}. 

1290 ''' 

1291 _xinstanceof(GeodesicExact, geodesic=geodesic) 

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

1293 geodesic = geodesic.copy(deep=False, name=_UNDER_(NN, geodesic.name)) # NOT _under! 

1294# _update_all(geodesic) 

1295 _GeodesicLineExact.__init__(self, geodesic, lat1, lon1, azi1, caps, 0, **name) 

1296 

1297 

1298def _Astroid(x, y): 

1299 '''(INTERNAL) Solve M{k^4 + 2 * k^3 - (x^2 + y^2 - 1) 

1300 * k^2 - (2 * k + 1) * y^2 = 0} for positive root k. 

1301 ''' 

1302 p = x**2 

1303 q = y**2 

1304 r = fsumf_(_1_0, q, p, _N_2_0) 

1305 if r > 0 or q: 

1306 # avoid possible division by zero when r = 0 

1307 # by multiplying s and t by r^3 and r, resp. 

1308 S = p * q / _4_0 # S = r^3 * s 

1309 if r: 

1310 r = r / _6_0 # /= chokes PyChecker 

1311 r3 = r**3 

1312 T3 = r3 + S 

1313 # discriminant of the quadratic equation for T3 is 

1314 # zero on the evolute curve p^(1/3) + q^(1/3) = 1 

1315 d = (r3 + T3) * S 

1316 if d < 0: 

1317 # T is complex, but u is defined for a real result 

1318 a = atan2(sqrt(-d), -T3) / _3_0 

1319 # There are 3 possible cube roots, choose the one which 

1320 # avoids cancellation. Note d < 0 implies that r < 0. 

1321 u = (cos(a) * _2_0 + _1_0) * r 

1322 else: 

1323 # pick the sign on the sqrt to maximize abs(T3) to 

1324 # minimize loss of precision due to cancellation. 

1325 if d: 

1326 T3 += _copysign(sqrt(d), T3) # T3 = (r * t)^3 

1327 # _cbrt always returns the real root, _cbrt(-8) = -2 

1328 u = _cbrt(T3) # T = r * t 

1329 if u: # T can be zero; but then r2 / T -> 0 

1330 u += r**2 / u 

1331 u += r 

1332 elif S: # d == T3**2 == S**2: sqrt(d) == abs(S) == abs(T3) 

1333 u = _cbrt(S * _2_0) # == T3 + _copysign(abs(S), T3) 

1334 else: 

1335 u = _0_0 

1336 v = _hypot(u, y) # sqrt(u**2 + q) 

1337 # avoid loss of accuracy when u < 0 

1338 u = (q / (v - u)) if u < 0 else (v + u) 

1339 w = (u - q) / (v + v) # positive? 

1340 # rearrange expression for k to avoid loss of accuracy due to 

1341 # subtraction, division by 0 impossible because u > 0, w >= 0 

1342 k = u / (sqrt(w**2 + u) + w) # guaranteed positive 

1343 

1344 else: # q == 0 && r <= 0 

1345 # y = 0 with |x| <= 1. Handle this case directly, for 

1346 # y small, positive root is k = abs(y) / sqrt(1 - x^2) 

1347 k = _0_0 

1348 

1349 return k 

1350 

1351 

1352def _C4coeffs(nC4): # in .geodesicx.__main__ 

1353 '''(INTERNAL) Get the C{C4_24}, C{_27} or C{_30} series coefficients. 

1354 ''' 

1355 try: # from pygeodesy.geodesicx._C4_xx import _coeffs_xx as _coeffs 

1356 _C4_xx = _DOT_(_MODS.geodesicx.__name__, _UNDER_('_C4', nC4)) 

1357 _coeffs = _MODS.getattr(_C4_xx, _UNDER_('_coeffs', nC4)) 

1358 except (AttributeError, ImportError, TypeError) as x: 

1359 raise GeodesicError(nC4=nC4, cause=x) 

1360 n = _xnC4(nC4=nC4) 

1361 if len(_coeffs) != n: # double check 

1362 raise GeodesicError(_coeffs=len(_coeffs), _xnC4=n, nC4=nC4) 

1363 return _coeffs 

1364 

1365 

1366__all__ += _ALL_DOCS(GeodesicExact, GeodesicLineExact) 

1367 

1368# **) MIT License 

1369# 

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

1371# 

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

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

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

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

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

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

1378# 

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

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

1381# 

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

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

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

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

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

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

1388# OTHER DEALINGS IN THE SOFTWARE.