1 2# -*- coding: utf-8 -*- 3 4u'''A Python version of I{Karney}'s C++ class U{GeodesicExact 5<https://GeographicLib.SourceForge.io/html/classGeographicLib_1_1GeodesicExact.html>}. 6 7Classes L{GeodesicExact} and L{GeodesicLineExact} follow the naming, 8methods and return values from I{Karney}s' Python classes C{Geodesic} 9and C{GeodesicLine}, respectively. 10 11Copyright (C) Charles Karney (2012-2021) <Charles@Karney.com> 12and licensed under the MIT/X11 License. For more information, 13see U{GeographicLib<https://GeographicLib.SourceForge.io>}. 14''' 15# make sure int/int division yields float quotient 16from __future__ import division 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# - beta = latitude on auxiliary sphere 25# - omega = longitude on auxiliary sphere 26# - lambda = longitude 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 copysign0, _xinstanceof, _xor, unsign0 40from pygeodesy.datums import _ellipsoidal_datum 41from pygeodesy.ellipsoids import Ellipsoid2, Ellipsoids 42from pygeodesy.fmath import cbrt, fsum_, hypot, norm2 43from pygeodesy.geodesicx.gxbases import _ALL_DOCS, Caps, _coSeries, \ 44 _GeodesicBase, _polynomial, \ 45 _sincos12, _TINY, _xnC4 46from pygeodesy.geodesicx.gxline import _GeodesicLineExact, pairs, \ 47 _update_glXs 48from pygeodesy.interns import EPS, MANT_DIG, NAN, NN, PI, PI_2, \ 49 _COMMASPACE_, _convergence_, _EPSqrt, \ 50 _no_, _0_0, _0_01, _0_5, _1_0, _2_0, \ 51 _3_0, _4_0, _6_0, _8_0, _16_0, _90_0, \ 52 _180_0, _1000_0 53from pygeodesy.interns import _0_001, _0_1 # PYCHOK used! 54from pygeodesy.karney import _around, GDict, GeodesicError, \ 55 _diff182, _fix90, _norm180 56from pygeodesy.namedTuples import Destination3Tuple, Distance3Tuple 57from pygeodesy.props import Property, Property_RO 58# from pygeodesy.streprs import pairs # from .geodesicx.gxline 59from pygeodesy.utily import atan2d, sincos2, sincos2d, unroll180, wrap360 60 61from math import atan2, cos, degrees, radians, sqrt 62 63__all__ = () 64__version__ = '21.07.15' 65 66_MAXIT1 = 20 67_MAXIT2 = 10 + _MAXIT1 + MANT_DIG # MANT_DIG == C++ digits 68_SQRT2_2 = -0.7071 # negative sqrt(2) / 2 69_1_75 = 1.75 70 71# increased multiplier in defn of _TOL1 from 100 to 200 to fix Inverse 72# case 52.784459512564 0 -52.784459512563990912 179.634407464943777557 73# which otherwise failed for Visual Studio 10 (Release and Debug) 74_TOL0 = EPS 75_TOL1 = _TOL0 * 200 76_TOL2 = _EPSqrt # == sqrt(_TOL0) 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 85class _PDict(GDict): 86 '''(INTERNAL) Parameters passed around in C{._GDictInverse} and 87 optionally returned when C{GeodesicExact.debug} is C{True}. 88 ''' 89 def setsigs(self, ssig1, csig1, ssig2, csig2): 90 '''Update the C{sig1} and C{sig2} parameters. 91 ''' 92 self.set_(ssig1=ssig1, csig1=csig1, sncndn1=(ssig1, csig1, self.dn1), # PYCHOK dn1 93 ssig2=ssig2, csig2=csig2, sncndn2=(ssig2, csig2, self.dn2)) # PYCHOK dn2 94 95 def toGDict(self): 96 '''Return as C{GDict} without attrs C{sncndn1} and C{sncndn2}. 97 ''' 98 def _rest(sncndn1=None, sncndn2=None, **rest): # PYCHOK unused 99 return GDict(rest) 100 101 return _rest(**self) 102 103 104class GeodesicExact(_GeodesicBase): 105 '''A pure Python version of I{Karney}'s C++ class U{GeodesicExact 106 <https://GeographicLib.SourceForge.io/html/classGeographicLib_1_1GeodesicExact.html>}, 107 modeled after I{Karney}'s Python class U{Geodesic<https://GeographicLib.SourceForge.io/ 108 html/python/code.html#module-geographiclib.geodesic>}. 109 ''' 110 _E = Ellipsoids.WGS84 111 _nC4 = 30 # default C4Order 112 113 def __init__(self, a_ellipsoid=Ellipsoids.WGS84, f=None, name=NN, C4Order=None): 114 '''New L{GeodesicExact} instance. 115 116 @arg a_ellipsoid: An ellipsoid (L{Ellipsoid}) or datum (L{datum}) 117 the equatorial radius of the ellipsoid (C{meter}), 118 see B{C{f}}. 119 @arg f: The flattening of the ellipsoid (C{scalar}) if B{C{a_ellipsoid}} 120 is specified as C{meter}. 121 @kwarg name: Optional name (C{str}). 122 @kwarg C4Order: Optional series expansion order (C{int}), see property 123 L{C4Order}, default C{30}. 124 125 @raise GeodesicError: Invalid B{C{C4Order}}. 126 ''' 127 if a_ellipsoid in (GeodesicExact._E, None): 128 pass # ignore f, default WGS84 129 elif f is None: 130 self._E = _ellipsoidal_datum(a_ellipsoid, name=name, raiser=True).ellipsoid 131 else: 132 self._E = Ellipsoid2(a_ellipsoid, f, name=name) 133 134 if name: 135 self.name = name 136 if C4Order: # XXX private copy, always? 137 self.C4Order = C4Order 138 139 @Property_RO 140 def a(self): 141 '''Get the I{equatorial} radius, semi-axis (C{meter}). 142 ''' 143 return self.ellipsoid.a 144 145 def ArcDirect(self, lat1, lon1, azi1, a12, outmask=Caps.STANDARD): 146 '''Solve the I{Direct} geodesic problem in terms of (spherical) arc length. 147 148 @arg lat1: Latitude of the first point (C{degrees}). 149 @arg lon1: Longitude of the first point (C{degrees}). 150 @arg azi1: Azimuth at the first point (compass C{degrees}). 151 @arg a12: Arc length between the points (C{degrees}), can be negative. 152 @kwarg outmask: Bit-or'ed combination of L{Caps} values specifying 153 the quantities to be returned. 154 155 @return: A C{dict} with up to 12 items C{lat1, lon1, azi1, lat2, 156 lon2, azi2, m12, a12, s12, M12, M21, S12} with C{lat1}, 157 C{lon1}, C{azi1} and arc length C{a12} always included. 158 159 @see: C++ U{GeodesicExact.ArcDirect 160 <https://GeographicLib.SourceForge.io/html/classGeographicLib_1_1GeodesicExact.html>} 161 and Python U{Geodesic.ArcDirect<https://GeographicLib.SourceForge.io/html/python/code.html>}. 162 ''' 163 return self._GDictDirect(lat1, lon1, azi1, True, a12, outmask) 164 165 def ArcDirectLine(self, lat1, lon1, azi1, a12, caps): 166 '''Define an L{GeodesicLineExact} in terms of the I{direct} geodesic problem and as arc length. 167 168 @arg lat1: Latitude of the first point (C{degrees}). 169 @arg lon1: Longitude of the first point (C{degrees}). 170 @arg azi1: Azimuth at the first point (compass C{degrees}). 171 @arg a12: Arc length between the points (C{degrees}), can be negative. 172 @kwarg caps: Bit-or'ed combination of L{Caps} values specifying 173 the capabilities the L{GeodesicLineExact} instance 174 should possess, i.e., which quantities can be 175 returned by calls to L{GeodesicLineExact.Position} 176 and L{GeodesicLineExact.ArcPosition}. 177 178 @return: A L{GeodesicLineExact} instance. 179 180 @note: The third point of the L{GeodesicLineExact} is set to correspond 181 to the second point of the I{Inverse} geodesic problem. 182 183 @note: Latitude B{C{lat1}} should in the range C{[-90, +90]}. 184 185 @see: C++ U{GeodesicExact.ArcDirectLine 186 <https://GeographicLib.SourceForge.io/html/classGeographicLib_1_1GeodesicExact.html>} and 187 Python U{Geodesic.ArcDirectLine<https://GeographicLib.SourceForge.io/html/python/code.html>}. 188 ''' 189 return self._GenDirectLine(lat1, lon1, azi1, True, a12, caps) 190 191 def Area(self, polyline=False, name=NN): 192 '''Set up an L{GeodesicAreaExact} to compute area and 193 perimeter of a polygon. 194 195 @kwarg polyline: If C{True} perimeter only, otherwise 196 area and perimeter (C{bool}). 197 @kwarg name: Optional name (C{str}). 198 199 @return: A L{GeodesicAreaExact} instance. 200 201 @note: The B{C{debug}} setting is passed as C{verbose} 202 to the returned L{GeodesicAreaExact} instance. 203 ''' 204 from pygeodesy.geodesicx.gxarea import GeodesicAreaExact 205 gaX = GeodesicAreaExact(self, polyline=polyline, 206 name=name or self.name) 207 if self.debug: 208 gaX.verbose = True 209 return gaX 210 211 Polygon = Area # for C{geographiclib} compatibility 212 213 @Property_RO 214 def b(self): 215 '''Get the ellipsoid's I{polar} radius, semi-axis (C{meter}). 216 ''' 217 return self.ellipsoid.b 218 219 @Property_RO 220 def c2x(self): 221 '''Get the ellipsoid's I{authalic} earth radius I{squared} (C{meter**2}). 222 ''' 223 # The Geodesic class substitutes atanh(sqrt(e2)) for asinh(sqrt(ep2)) 224 # in the definition of _c2. The latter is more accurate for very 225 # oblate ellipsoids (which the Geodesic class does not handle). Of 226 # course, the area calculation in GeodesicExact is still based on a 227 # series and only holds for moderately oblate (or prolate) ellipsoids. 228 return self.ellipsoid.c2x 229 230 c2 = c2x # in this particular case 231 232 def C4f(self, eps): 233 '''Evaluate the C{C4x} coefficients for B{C{eps}}. 234 235 @arg eps: Polynomial factor (C{float}). 236 237 @return: C{C4Order}-Tuple of C{C4x(B{eps})} coefficients. 238 ''' 239 def _c4(nC4, C4x, _p_eps, eps): 240 i, e = 0, _1_0 241 for r in range(nC4, 0, -1): 242 j = i + r 243 yield _p_eps(eps, C4x, i, j) * e 244 e *= eps 245 i = j 246 # assert i == (nC4 * (nC4 + 1)) // 2 247 248 return tuple(_c4(self._nC4, self._C4x, 249 _polynomial, eps)) 250 251 def _C4f_k2(self, k2): # in ._GDictInverse and gxline._GeodesicExactLine._C4a 252 '''(INTERNAL) Compute C{eps} from B{C{k2}} and invoke C{C4f}. 253 ''' 254 return self.C4f(k2 / fsum_(_2_0, sqrt(k2 + _1_0) * _2_0, k2)) 255 256 @Property 257 def C4Order(self): 258 '''Get the series expansion order (C{int}, 24, 27 or 30). 259 ''' 260 return self._nC4 261 262 @C4Order.setter # PYCHOK .setter 263 def C4Order(self, order): 264 '''Set the series expansion order. 265 266 @arg order: New order (C{int}, 24, 27 or 30). 267 268 @raise GeodesicError: Invalid B{C{order}}. 269 ''' 270 _xnC4(C4Order=order) 271 if self._nC4 != order: 272 GeodesicExact._C4x._update(self) 273 _update_glXs(self) 274 self._nC4 = order 275 276 @Property_RO 277 def _C4x(self): 278 '''Get this ellipsoid's C{C4} coefficients, I{cached} tuple. 279 280 @see: Property L{nC4}. 281 ''' 282 # see C4coeff() in GeographicLib.src.GeodesicExactC4.cpp 283 def _C4(nC4, coeffs, _p_n, n): 284 i = 0 285 for r in range(nC4 + 1, 1, -1): 286 for j in range(1, r): 287 j = i + j # (j - i - 1) order of polynomial 288 yield _p_n(n, coeffs, i, j) / coeffs[j] 289 i = j + 1 290 # assert i == len(cs) == (nC4 * (nC4 + 1) * (nC4 + 5)) // 6 291 292 return tuple(_C4(self._nC4, self._coeffs(self._nC4), 293 _polynomial, self.n)) # 3rd flattening 294 295 def _coeffs(self, nC4): 296 '''(INTERNAL) Get the series coefficients. 297 ''' 298 if nC4 == 30: 299 from pygeodesy.geodesicx._C4_30 import _coeffs_30 as _coeffs 300 elif nC4 == 27: 301 from pygeodesy.geodesicx._C4_27 import _coeffs_27 as _coeffs 302 elif nC4 == 24: 303 from pygeodesy.geodesicx._C4_24 import _coeffs_24 as _coeffs 304 else: # PYCHOK no cover 305 _coeffs = () 306 n = _xnC4(nC4=nC4) 307 if len(_coeffs) != n: # double check 308 raise GeodesicError(_coeffs=len(_coeffs), _xnC4=n, nC4=nC4) 309 return _coeffs 310 311 def Direct(self, lat1, lon1, azi1, s12, outmask=Caps.STANDARD): 312 '''Solve the I{Direct} geodesic problem 313 314 @arg lat1: Latitude of the first point (C{degrees}). 315 @arg lon1: Longitude of the first point (C{degrees}). 316 @arg azi1: Azimuth at the first point (compass C{degrees}). 317 @arg s12: Distance between the points (C{meter}), can be negative. 318 @kwarg outmask: Bit-or'ed combination of L{Caps} values specifying 319 the quantities to be returned. 320 321 @return: A C{dict} with up to 12 items C{lat1, lon1, azi1, lat2, 322 lon2, azi2, m12, a12, s12, M12, M21, S12} with C{lat1}, 323 C{lon1}, C{azi1} and distance C{s12} always included. 324 325 @see: C++ U{GeodesicExact.Direct 326 <https://GeographicLib.SourceForge.io/html/classGeographicLib_1_1GeodesicExact.html>} 327 and Python U{Geodesic.Direct<https://GeographicLib.SourceForge.io/html/python/code.html>}. 328 ''' 329 return self._GDictDirect(lat1, lon1, azi1, False, s12, outmask) 330 331 def Direct3(self, lat1, lon1, azi1, s12): # PYCHOK outmask 332 '''Return the destination lat, lon and reverse azimuth 333 (final bearing) in C{degrees}. 334 335 @return: L{Destination3Tuple}C{(lat, lon, final)}. 336 ''' 337 r = self._GDictDirect(lat1, lon1, azi1, False, s12, Caps._AZIMUTH_LATITUDE_LONGITUDE) 338 return Destination3Tuple(r.lat2, r.lon2, r.azi2) 339 340 def DirectLine(self, lat1, lon1, azi1, s12, caps=Caps.STANDARD): 341 '''Define a {GeodesicLineExact} in terms of the I{direct} geodesic problem and as distance. 342 343 @arg lat1: Latitude of the first point (C{degrees}). 344 @arg lon1: Longitude of the first point (C{degrees}). 345 @arg azi1: Azimuth at the first point (compass C{degrees}). 346 @arg s12: Distance between the points (C{meter}), can be negative. 347 @kwarg caps: Bit-or'ed combination of L{Caps} values specifying 348 the capabilities the L{GeodesicLineExact} instance 349 should possess, i.e., which quantities can be 350 returned by calls to L{GeodesicLineExact.Position}. 351 352 @return: A L{GeodesicLineExact} instance. 353 354 @note: The third point of the L{GeodesicLineExact} is set to correspond 355 to the second point of the I{Inverse} geodesic problem. 356 357 @note: Latitude B{C{lat1}} should in the range C{[-90, +90]}. 358 359 @see: C++ U{GeodesicExact.DirectLine 360 <https://GeographicLib.SourceForge.io/html/classGeographicLib_1_1GeodesicExact.html>} and 361 Python U{Geodesic.DirectLine<https://GeographicLib.SourceForge.io/html/python/code.html>}. 362 ''' 363 return self._GenDirectLine(lat1, lon1, azi1, False, s12, caps) 364 365 def _dn(self, sbet, cbet): # in gxline._GeodesicLineExact.__init__ 366 '''(INTERNAL) Helper. 367 ''' 368 return sqrt(_1_0 + self.ep2 * sbet**2) if self.f >= 0 else ( 369 sqrt(_1_0 - self.e2 * cbet**2) / self.f1) 370 371 @Property_RO 372 def e2(self): 373 '''Get the ellipsoid's I{(1st) eccentricity squared} (C{float}), M{f * (2 - f)}. 374 ''' 375 return self.ellipsoid.e2 376 377 @Property_RO 378 def _e2a2(self): 379 '''(INTERNAL) Cache M{E.e2 * E.a2}. 380 ''' 381 return self.e2 * self.ellipsoid.a2 382 383 @Property_RO 384 def _e2_f1(self): 385 '''(INTERNAL) Cache M{E.e2 * E.f1}. 386 ''' 387 return self.e2 / self.f1 388 389 @Property_RO 390 def _eF(self): 391 '''(INTERNAL) Get the elliptic function, aka C{.E}. 392 ''' 393 from pygeodesy.elliptic import Elliptic 394 return Elliptic(k2=-self.ep2) 395 396 def _eF_reset_e2_f1_cH(self, x, y): 397 '''(INTERNAL) Reset elliptic function and return M{e2 / f1 * cH * ...}. 398 ''' 399 self._eF_reset_k2(x) 400 return self._e2_f1 * y * self._eF.cH 401 402 def _eF_reset_k2(self, x): 403 '''(INTERNAL) Reset elliptic function and return C{k2}. 404 ''' 405 ep2 = self.ep2 406 k2 = x**2 * ep2 407 self._eF.reset(k2=-k2, alpha2=-ep2, 408 kp2=k2 + _1_0, alphap2=ep2 + _1_0) # defaults 409 _update_glXs(self) # zap cached/memoized _GeodesicLineExact attrs 410 return k2 411 412 @Property_RO 413 def ellipsoid(self): 414 '''Get the ellipsoid (C{Ellipsoid}). 415 ''' 416 return self._E 417 418 @Property_RO 419 def ep2(self): # aka .e22 420 '''Get the ellipsoid's I{2nd eccentricity squared} (C{float}), M{e2 / (1 - e2)}. 421 ''' 422 return self.ellipsoid.e22 # == self.e2 / self.f1**2 423 424 @Property_RO 425 def _eTOL2(self): 426 '''(INTERNAL) The si12 threshold for "really short". 427 ''' 428 # Using the auxiliary sphere solution with dnm computed at 429 # (bet1 + bet2) / 2, the relative error in the azimuth 430 # consistency check is sig12^2 * abs(f) * min(1, 1-f/2) / 2. 431 # (Error measured for 1/100 < b/a < 100 and abs(f) >= 1/1000. 432 # For a given f and sig12, the max error occurs for lines 433 # near the pole. If the old rule for computing dnm = (dn1 434 # + dn2)/2 is used, then the error increases by a factor of 435 # 2.) Setting this equal to epsilon gives sig12 = etol2. 436 # Here 0.1 is a safety factor (error decreased by 100) and 437 # max(0.001, abs(f)) stops etol2 getting too large in the 438 # nearly spherical case. 439 f = self.f 440 return _0_1 * _TOL2 / sqrt(min(_1_0, _1_0 - f * _0_5) * 441 max(_0_001, abs(f)) * _0_5) 442 443 @Property_RO 444 def f(self): 445 '''Get the ellipsoid's I{flattening} (C{float}), M{(a - b) / a}, C{0} for spherical, negative for prolate. 446 ''' 447 return self.ellipsoid.f 448 449 flattening = f 450 451 @Property_RO 452 def f1(self): # aka .b_a 453 '''Get the ellipsoid's ratio I{polar} over I{equatorial} radius (C{float}), M{b / a} == M{1 - f}. 454 ''' 455 return self.ellipsoid.b_a # _1_0 - self.f 456 457 @Property_RO 458 def _f_180(self): 459 '''(INTERNAL) Cached/memoized. 460 ''' 461 return self.f * _180_0 462 463 def _GDictDirect(self, lat1, lon1, azi1, arcmode, s12_a12, outmask=Caps.STANDARD): 464 '''(INTERNAL) As C{_GenDirect}, but returning a C{dict}. 465 466 @return: A C{dict} ... 467 ''' 468 if not arcmode: # supply DISTANCE_IN 469 outmask |= Caps.DISTANCE_IN 470 glX = self._Line(lat1, lon1, azi1, outmask) 471 return glX._GDictPosition(arcmode, s12_a12, outmask) 472 473 def _GDictInverse(self, lat1, lon1, lat2, lon2, outmask=Caps.STANDARD): # MCCABE 32, 41 vars 474 '''(INTERNAL) As C{_GenInverse}, but returning a C{dict}. 475 476 @return: A C{dict} ... 477 ''' 478 if self.debug: # PYCHOK no cover 479 outmask |= Caps._DEBUG_INVERSE & self._debug 480 outmask &= Caps._OUTMASK # incl. _SALPs_CALPs and _DEBUG_ 481 # compute longitude difference carefully (with _diff182): 482 # result is in [-180, +180] but -180 is only for west-going 483 # geodesics, +180 is for east-going and meridional geodesics 484 lon12, lon12s = _diff182(lon1, lon2) 485 # see C{result} from geographiclib.geodesic.Inverse 486 if (outmask & Caps.LONG_UNROLL): # == (lon1 + lon12) + lon12s 487 r = GDict(lon1=lon1, lon2=fsum_(lon1, lon12, lon12s)) 488 else: 489 r = GDict(lon1=_norm180(lon1), lon2=_norm180(lon2)) 490 # make longitude difference positive and if very close to 491 # being on the same half-meridian, then make it so. 492 if lon12 < 0: 493 lon_, lon12 = True, -_around(lon12) 494 lon12s = _around((_180_0 - lon12) + lon12s) 495 else: 496 lon_, lon12 = False, _around(lon12) 497 lon12s = _around((_180_0 - lon12) - lon12s) 498 lam12 = radians(lon12) 499 if lon12 > _90_0: 500 slam12, clam12 = sincos2d(lon12s) 501 clam12 = -clam12 502 else: 503 slam12, clam12 = sincos2d(lon12) 504 505 # If really close to the equator, treat as on equator. 506 lat1 = _around(_fix90(lat1)) 507 lat2 = _around(_fix90(lat2)) 508 r.set_(lat1=lat1, lat2=lat2) 509 # Swap points so that point with higher (abs) latitude is 510 # point 1. If one latitude is a NAN, then it becomes lat1. 511 swap_ = abs(lat1) < abs(lat2) 512 if swap_: 513 lat1, lat2 = lat2, lat1 514 lon_ = not lon_ 515 if lat1 < 0: 516 lat_ = False 517 else: # make lat1 <= 0 518 lat_ = True 519 lat1, lat2 = -lat1, -lat2 520 # Now 0 <= lon12 <= 180, -90 <= lat1 <= 0 and lat1 <= lat2 <= -lat1 521 # and lat_, lon_, swap_ register the transformation to bring the 522 # coordinates to this canonical form, where False means no change 523 # made. We make these transformations so that there are few cases 524 # to check, e.g., on verifying quadrants in atan2. In addition, 525 # this enforces some symmetries in the results returned. 526 527 # Initialize for the meridian. No longitude calculation is 528 # done in this case to let the parameter default to 0. 529 sbet1, cbet1 = self._sincos2f1(lat1) 530 sbet2, cbet2 = self._sincos2f1(lat2) 531 # If cbet1 < -sbet1, then cbet2 - cbet1 is a sensitive measure 532 # of the |bet1| - |bet2|. Alternatively (cbet1 >= -sbet1), 533 # abs(sbet2) + sbet1 is a better measure. This logic is used 534 # in assigning calp2 in _Lambda10. Sometimes these quantities 535 # vanish and in that case we force bet2 = +/- bet1 exactly. An 536 # example where is is necessary is the inverse problem 537 # 48.522876735459 0 -48.52287673545898293 179.599720456223079643 538 # which failed with Visual Studio 10 (Release and Debug) 539 if cbet1 < -sbet1: 540 if cbet2 == cbet1: 541 sbet2 = sbet1 if sbet2 < 0 else -sbet1 542 elif abs(sbet2) == -sbet1: 543 cbet2 = cbet1 544 545 p = _PDict(sbet1=sbet1, cbet1=cbet1, dn1=self._dn(sbet1, cbet1), 546 sbet2=sbet2, cbet2=cbet2, dn2=self._dn(sbet2, cbet2)) 547 548 _meridian = True # i.e. meridian = False 549 if lat1 == -90 or slam12 == 0: 550 # Endpoints are on a single full meridian, 551 # so the geodesic might lie on a meridian. 552 salp1, calp1 = slam12, clam12 # head to target lon 553 salp2, calp2 = _0_0, _1_0 # then head north 554 # tan(bet) = tan(sig) * cos(alp) 555 p.setsigs(sbet1, calp1 * cbet1, sbet2, calp2 * cbet2) 556 # sig12 = sig2 - sig1 557 sig12 = atan2(*_sincos12(sbet1, p.csig1, sbet2, p.csig2, noneg=True)) # PYCHOK csig* 558 s12x, m12x, _, \ 559 M12, M21 = self._Lengths5(sig12, outmask | Caps.REDUCEDLENGTH, p) 560 # Add the check for sig12 since zero length geodesics 561 # might yield m12 < 0. Test case was 562 # echo 20.001 0 20.001 0 | GeodSolve -i 563 # In fact, we will have sig12 > PI/2 for meridional 564 # geodesic which is not a shortest path. 565 if m12x >= 0 or sig12 < _1_0: 566 # Need at least 2 to handle 90 0 90 180 567 # Prevent negative s12 or m12 from geographiclib 1.52 568 if sig12 < _TINY3 or (sig12 < _TOL0 and (s12x < 0 or m12x < 0)): 569 sig12 = m12x = s12x = a12 = _0_0 570 else: 571 m12x *= self.b 572 s12x *= self.b 573 a12 = degrees(sig12) 574 _meridian = False 575 C = 1 576 # elif m12 < 0: # prolate and too close to anti-podal 577 # _meridian = True 578 579 if _meridian: 580 if sbet1 == 0 and ( # sbet2 == 0 and 581 self.f <= 0 or lon12s >= self._f_180): 582 # Geodesic runs along equator 583 calp1 = calp2 = _0_0 584 salp1 = salp2 = _1_0 585 sig12 = lam12 / self.f1 # == omg12 586 somg12, comg12 = sincos2(sig12) 587 m12x = self.b * somg12 588 s12x = self.a * lam12 589 if (outmask & Caps.GEODESICSCALE): 590 r.set_(M12=comg12, M21=comg12) 591 a12 = lon12 / self.f1 592 C = 2 593 else: 594 # Now point1 and point2 belong within a hemisphere bounded by a 595 # meridian and geodesic is neither meridional or equatorial. 596 p.set_(slam12=slam12, clam12=clam12) 597 # Figure a starting point for Newton's method 598 sig12, salp1, calp1, \ 599 salp2, calp2, dnm = self._InverseStart6(lam12, p) 600 if sig12 is None: # use Newton's method 601 # pre-compute the constant _Lambda10 term, once 602 p.set_(bet12=None if cbet2 == cbet1 and abs(sbet2) == -sbet1 else ( 603 ((cbet1 + cbet2) * (cbet2 - cbet1)) if cbet1 < -sbet1 else 604 ((sbet1 + sbet2) * (sbet1 - sbet2)))) 605 sig12, salp1, calp1, salp2, calp2, domg12, \ 606 s12x, m12x, M12, M21 = self._Newton10(salp1, calp1, outmask, p) 607 if (outmask & Caps.AREA): 608 # omg12 = lam12 - domg12 609 s, c = sincos2(domg12) 610 somg12, comg12 = _sincos12(s, c, slam12, clam12) 611 612 if (outmask & Caps.GEODESICSCALE): 613 if swap_: 614 M12, M21 = M21, M12 615 r.set_(M12=unsign0(M12), 616 M21=unsign0(M21)) 617 C = 3 # Newton 618 else: 619 C = 4 # Short lines, _Inverse6 set salp2, calp2, dnm 620 s, c = sincos2(sig12 / dnm) 621 m12x = s * dnm**2 622 s12x = sig12 * dnm 623 if (outmask & Caps.GEODESICSCALE): 624 r.set_(M12=c, M21=c) 625 if (outmask & Caps.AREA): 626 somg12, comg12 = sincos2(lam12 / (self.f1 * dnm)) 627 a12 = degrees(sig12) 628 m12x *= self.b 629 s12x *= self.b 630 631 else: # _meridian is False 632 somg12 = comg12 = NAN 633 634 r.set_(a12=a12) # in [0, 180] 635 if outmask == Caps._ANGLE_ONLY: 636 return r # for .Inverse1 637 638 if (outmask & Caps.DISTANCE): 639 r.set_(s12=unsign0(s12x)) 640 641 if (outmask & Caps.REDUCEDLENGTH): 642 r.set_(m12=unsign0(m12x)) 643 644 if (outmask & Caps.AREA): 645 S12 = self._InverseArea(_meridian, salp1, calp1, 646 salp2, calp2, 647 somg12, comg12, p) 648 if _xor(swap_, lat_, lon_): 649 S12 = -S12 650 r.set_(S12=unsign0(S12)) 651 652 if swap_: 653 salp1, salp2 = salp2, salp1 654 calp1, calp2 = calp2, calp1 655 if _xor(swap_, lon_): 656 salp1, salp2 = -salp1, -salp2 657 if _xor(swap_, lat_): 658 calp1, calp2 = -calp1, -calp2 659 660 if (outmask & Caps.AZIMUTH): 661 r.set_(azi1=atan2d(salp1, calp1), 662 azi2=atan2d(salp2, calp2, reverse=outmask & Caps.REVERSE2)) 663 if (outmask & Caps._SALPs_CALPs): 664 r.set_(salp1=salp1, calp1=calp1, 665 salp2=salp2, calp2=calp2) 666 667 if (outmask & Caps._DEBUG_INVERSE): # PYCHOK no cover 668 eF = self._eF 669 p.set_(C=C, a=self.a, f=self.f, f1=self.f1, 670 e=self.ellipsoid.e, e2=self.e2, ep2=self.ep2, 671 c2x=self.c2x, c2=self.ellipsoid.c2, 672 eFcD=eF.cD, eFcE=eF.cE, eFcH=eF.cH, 673 eFk2=eF.k2, eFa2=eF.alpha2) 674 p.update(r) # r overrides 675 r = p.toGDict() 676 return r 677 678 def _GenDirect(self, lat1, lon1, azi1, arcmode, s12_a12, outmask): 679 '''(INTERNAL) The general I{Inverse} geodesic calculation. 680 681 @return: L{Direct9Tuple}C{(a12, lat2, lon2, azi2, 682 s12, m12, M12, M21, S12)}. 683 ''' 684 r = self._GDictDirect(lat1, lon1, azi1, arcmode, s12_a12, outmask) 685 return r.toDirect9Tuple() 686 687 def _GenDirectLine(self, lat1, lon1, azi1, arcmode, s12_a12, caps): 688 '''(INTERNAL) Helper for C{ArcDirectLine} and C{DirectLine}. 689 690 @return: A L{GeodesicLineExact} instance. 691 ''' 692 azi1 = _norm180(azi1) 693 # guard against underflow in salp0. Also -0 is converted to +0. 694 s, c = sincos2d(_around(azi1)) 695 696 if not arcmode: # supply DISTANCE_IN 697 caps |= Caps.DISTANCE_IN 698 return _GeodesicLineExact(self, lat1, lon1, azi1, caps, 699 self._debug, s, c)._GenSet(arcmode, s12_a12) 700 701 def _GenInverse(self, lat1, lon1, lat2, lon2, outmask): 702 '''(INTERNAL) The general I{Inverse} geodesic calculation. 703 704 @return: L{Inverse10Tuple}C{(a12, s12, salp1, calp1, salp2, calp2, 705 m12, M12, M21, S12)}. 706 ''' 707 r = self._GDictInverse(lat1, lon1, lat2, lon2, outmask | Caps._SALPs_CALPs) 708 return r.toInverse10Tuple() 709 710 def Inverse(self, lat1, lon1, lat2, lon2, outmask=Caps.STANDARD): 711 '''Perform the I{Inverse} geodesic calculation. 712 713 @arg lat1: Latitude of the first point (C{degrees}). 714 @arg lon1: Longitude of the first point (C{degrees}). 715 @arg lat2: Latitude of the second point (C{degrees}). 716 @arg lon2: Longitude of the second point (C{degrees}). 717 @kwarg outmask: Bit-or'ed combination of L{Caps} values specifying 718 the quantities to be returned. 719 720 @return: A C{dict} with up to 12 items C{lat1, lon1, azi1, lat2, 721 lon2, azi2, m12, a12, s12, M12, M21, S12} with C{lat1}, 722 C{lon1}, C{azi1} and distance C{s12} always included. 723 724 @note: The third point of the L{GeodesicLineExact} is set to correspond 725 to the second point of the I{Inverse} geodesic problem. 726 727 @note: Both B{C{lat1}} and B{C{lat2}} should in the range C{[-90, +90]}. 728 729 @see: C++ U{GeodesicExact.InverseLine 730 <https://GeographicLib.SourceForge.io/html/classGeographicLib_1_1GeodesicExact.html>} and 731 Python U{Geodesic.InverseLine<https://GeographicLib.SourceForge.io/html/python/code.html>}. 732 ''' 733 return self._GDictInverse(lat1, lon1, lat2, lon2, outmask) 734 735 def Inverse1(self, lat1, lon1, lat2, lon2, wrap=False): 736 '''Return the non-negative, I{angular} distance in C{degrees}. 737 ''' 738 # see .FrechetKarney.distance, .HausdorffKarney._distance 739 # and .HeightIDWkarney._distances 740 _, lon2 = unroll180(lon1, lon2, wrap=wrap) # self.LONG_UNROLL 741 return abs(self._GDictInverse(lat1, lon1, lat2, lon2, Caps._ANGLE_ONLY).a12) 742 743 def Inverse3(self, lat1, lon1, lat2, lon2): # PYCHOK outmask 744 '''Return the distance in C{meter} and the forward and 745 reverse azimuths (initial and final bearing) in C{degrees}. 746 747 @return: L{Distance3Tuple}C{(distance, initial, final)}. 748 ''' 749 r = self._GDictInverse(lat1, lon1, lat2, lon2, Caps._AZIMUTH_DISTANCE) 750 return Distance3Tuple(r.s12, wrap360(r.azi1), wrap360(r.azi2)) 751 752 def InverseLine(self, lat1, lon1, lat2, lon2, caps=Caps.STANDARD): 753 '''Define a {GeodesicLineExact} in terms of the I{Inverse} geodesic problem. 754 755 @arg lat1: Latitude of the first point (C{degrees}). 756 @arg lon1: Longitude of the first point (C{degrees}). 757 @arg lat2: Latitude of the second point (C{degrees}). 758 @arg lon2: Longitude of the second point (C{degrees}). 759 @kwarg caps: Bit-or'ed combination of L{Caps} values specifying 760 the capabilities the L{GeodesicLineExact} instance 761 should possess, i.e., which quantities can be 762 returned by calls to L{GeodesicLineExact.Position} 763 and L{GeodesicLineExact.ArcPosition}. 764 765 @return: A L{GeodesicLineExact} instance. 766 767 @note: The third point of the L{GeodesicLineExact} is set to correspond 768 to the second point of the I{Inverse} geodesic problem. 769 770 @note: Both B{C{lat1}} and B{C{lat2}} should in the range C{[-90, +90]}. 771 772 @see: C++ U{GeodesicExact.InverseLine 773 <https://GeographicLib.SourceForge.io/html/classGeographicLib_1_1GeodesicExact.html>} and 774 Python U{Geodesic.InverseLine<https://GeographicLib.SourceForge.io/html/python/code.html>}. 775 ''' 776 r = self._GDictInverse(lat1, lon1, lat2, lon2, Caps._SALPs_CALPs) # No need for AZIMUTH 777 azi1 = atan2d(r.salp1, r.calp1) 778 if (caps & (Caps._OUTMASK & Caps.DISTANCE_IN)): 779 caps |= Caps.DISTANCE # ensure that a12 is a distance 780 return _GeodesicLineExact(self, lat1, lon1, azi1, caps, 781 self._debug, r.salp1, r.calp1)._GenSet(True, r.a12) 782 783 def _InverseArea(self, _meridian, salp1, calp1, # PYCHOK 9 args 784 salp2, calp2, 785 somg12, comg12, p): 786 '''(INTERNAL) Split off from C{_GDictInverse} to reduce complexity/length. 787 788 @return: Area C{S12}. 789 ''' 790 # from _Lambda10: sin(alp1) * cos(bet1) = sin(alp0) 791 salp0 = salp1 * p.cbet1 792 calp0 = hypot(calp1, salp1 * p.sbet1) # calp0 > 0 793 if calp0 and salp0: 794 # from _Lambda10: tan(bet) = tan(sig) * cos(alp) 795 k2 = calp0**2 * self.ep2 796 C4a = self._C4f_k2(k2) 797 B41 = _coSeries(C4a, *norm2(p.sbet1, calp1 * p.cbet1)) 798 B42 = _coSeries(C4a, *norm2(p.sbet2, calp2 * p.cbet2)) 799 # multiplier = a^2 * e^2 * cos(alpha0) * sin(alpha0) 800 A4 = self._e2a2 * calp0 * salp0 801 S12 = A4 * (B42 - B41) 802 p.set_(A4=A4, B41=B41, B42=B42, k2=k2) 803 else: # avoid problems with indeterminate sig1, sig2 on equator 804 S12 = _0_0 805 806 if (_meridian and # omg12 < 3/4 * PI 807 comg12 > _SQRT2_2 and # lon diff not too big 808 (p.sbet2 - p.sbet1) < _1_75): # lat diff not too big 809 # use tan(Gamma/2) = tan(omg12/2) * 810 # (tan(bet1/2) + tan(bet2/2)) / 811 # (tan(bet1/2) * tan(bet2/2) + 1)) 812 # with tan(x/2) = sin(x) / (1 + cos(x)) 813 cbet1 = _1_0 + p.cbet1 814 cbet2 = _1_0 + p.cbet2 815 salp12 = (p.sbet1 * cbet2 + cbet1 * p.sbet2) * somg12 816 calp12 = (p.sbet1 * p.sbet2 + cbet1 * cbet2) * (comg12 + _1_0) 817 alp12 = _2_0 * atan2(salp12, calp12) 818 else: 819 # alp12 = alp2 - alp1, used in atan2, no need to normalize 820 salp12, calp12 = _sincos12(salp1, calp1, salp2, calp2) 821 # The right thing appears to happen if alp1 = +/-180 and 822 # alp2 = 0, viz salp12 = -0 and alp12 = -180. However, 823 # this depends on the sign being attached to 0 correctly. 824 # Following ensures the correct behavior. 825 if salp12 == 0 and calp12 < 0: 826 alp12 = copysign0(PI, calp1) 827 else: 828 alp12 = atan2(salp12, calp12) 829 830 p.set_(alp12=alp12) 831 return S12 + self.c2x * alp12 832 833 def _InverseStart6(self, lam12, p): 834 '''(INTERNAL) Return a starting point for Newton's method in 835 C{salp1} and C{calp1} indicated by C{sig12=None}. If 836 Newton's method doesn't need to be used, return also 837 C{salp2}, C{calp2}, C{dnm} and C{sig12} non-C{None}. 838 839 @return: 6-Tuple C{(sig12, salp1, calp1, salp2, calp2, dnm)} 840 and updated C{p.setsigs} if C{sig12==None}. 841 ''' 842 sig12 = None # use Newton 843 salp1 = calp1 = salp2 = calp2 = dnm = NAN 844 845 # bet12 = bet2 - bet1 in [0, PI) 846 sbet12, cbet12 = _sincos12(p.sbet1, p.cbet1, p.sbet2, p.cbet2) 847 shortline = cbet12 >= 0 and sbet12 < _0_5 and (p.cbet2 * lam12) < _0_5 848 if shortline: 849 # sin((bet1 + bet2)/2)^2 = (sbet1 + sbet2)^2 / ( 850 # (cbet1 + cbet2)^2 + (sbet1 + sbet2)^2) 851 t = (p.sbet1 + p.sbet2)**2 852 s = t / ((p.cbet1 + p.cbet2)**2 + t) 853 dnm = sqrt(_1_0 + self.ep2 * s) 854 somg12, comg12 = sincos2(lam12 / (self.f1 * dnm)) 855 else: 856 somg12, comg12 = p.slam12, p.clam12 857 858 # bet12a = bet2 + bet1 in (-PI, 0], note -sbet1 859 sbet12a, cbet12a = _sincos12(-p.sbet1, p.cbet1, p.sbet2, p.cbet2) 860 861 s = somg12**2 / (abs(comg12) + _1_0) 862 t = p.cbet2 * p.sbet1 * s 863 salp1 = p.cbet2 * somg12 864 calp1 = (sbet12a - t) if comg12 < 0 else (sbet12 + t) 865 866 ssig12 = hypot(salp1, calp1) 867 csig12 = p.sbet1 * p.sbet2 + p.cbet1 * p.cbet2 * comg12 868 869 if shortline and ssig12 < self._eTOL2: # really short lines 870 if comg12 < 0: 871 s = _1_0 - comg12 872 salp2, calp2 = norm2(somg12 * p.cbet1, 873 sbet12 - p.cbet1 * p.sbet2 * s) 874 sig12 = atan2(ssig12, csig12) # return value 875 876 elif (self._n_0_1 or # Skip astroid calc if too eccentric 877 csig12 >= 0 or ssig12 >= (self._n_6PI * p.cbet1**2)): 878 pass # nothing to do, 0th order spherical approximation OK 879 880 else: 881 # Scale lam12 and bet2 to x, y coordinate system where antipodal 882 # point is at origin and singular point is at y = 0, x = -1 883 lam12x = atan2(-p.slam12, -p.clam12) # lam12 - PI 884 if self.f < 0: # x = dlat, y = dlon 885 # ssig1=sbet1, csig1=-cbet1, ssig2=sbet2, csig2=cbet2 886 p.setsigs(p.sbet1, -p.cbet1, p.sbet2, p.cbet2) 887 # if lon12 = 180, this repeats a calculation made in Inverse 888 _, m12b, m0, _, _ = self._Lengths5(atan2(sbet12a, cbet12a) + PI, 889 Caps.REDUCEDLENGTH, p) 890 x = m12b / (p.cbet1 * p.cbet2 * m0 * PI) - _1_0 891 sca = ((sbet12a / x) if x < -_0_01 else (-self.f * p.cbet1**2 * PI)) / p.cbet1 892 y = lam12x / sca 893 else: # _f >= 0, in fact f == 0 does not get here 894 sca = self._eF_reset_e2_f1_cH(p.sbet1, p.cbet1 * _2_0) 895 x = lam12x / sca # dlon 896 y = sbet12a / (sca * p.cbet1) # dlat 897 898 if y > -_TOL1 and x > -_THR1: # strip near cut 899 if self.f < 0: 900 calp1 = max((_0_0 if x > -_TOL1 else -_1_0), x) 901 salp1 = sqrt(_1_0 - calp1**2) 902 else: 903 salp1 = min(_1_0, -x) 904 calp1 = -sqrt(_1_0 - salp1**2) 905 else: 906 # Estimate alp1, by solving the astroid problem. 907 # 908 # Could estimate alpha1 = theta + PI/2, directly, i.e., 909 # calp1 = y/k; salp1 = -x/(1+k); for _f >= 0 910 # calp1 = x/(1+k); salp1 = -y/k; for _f < 0 (need to check) 911 # 912 # However, it's better to estimate omg12 from astroid and use 913 # spherical formula to compute alp1. This reduces the mean 914 # number of Newton iterations for astroid cases from 2.24 915 # (min 0, max 6) to 2.12 (min 0 max 5). The changes in the 916 # number of iterations are as follows: 917 # 918 # change percent 919 # 1 5 920 # 0 78 921 # -1 16 922 # -2 0.6 923 # -3 0.04 924 # -4 0.002 925 # 926 # The histogram of iterations is (m = number of iterations 927 # estimating alp1 directly, n = number of iterations 928 # estimating via omg12, total number of trials = 148605): 929 # 930 # iter m n 931 # 0 148 186 932 # 1 13046 13845 933 # 2 93315 102225 934 # 3 36189 32341 935 # 4 5396 7 936 # 5 455 1 937 # 6 56 0 938 # 939 # omg12 is near PI, estimate work with omg12a = PI - omg12 940 k = _Astroid(x, y) 941 k1 = k + _1_0 942 sca *= (y * k1 / k) if self.f < 0 else (x * k / k1) 943 s, c = sincos2(-sca) 944 # update spherical estimate of alp1 using omg12 instead of lam12 945 salp1 = p.cbet2 * s 946 calp1 = sbet12a - p.cbet2 * p.sbet1 * s**2 / (_1_0 + c) 947 948 # sanity check on starting guess. Backwards check allows NaN through. 949 salp1, calp1 = norm2(salp1, calp1) if salp1 > 0 else (_1_0, _0_0) 950 951 return sig12, salp1, calp1, salp2, calp2, dnm 952 953 def _Lambda5(self, salp1, calp1, diffp, p): 954 '''(INTERNAL) Helper. 955 956 @return: 5-Tuple C{(lam12, sig12, salp2, calp2, dlam12)} 957 and updated C{p.setsigs} and C{p.domg12}. 958 ''' 959 eF = self._eF 960 f1 = self.f1 961 962 if p.sbet1 == 0 and calp1 == 0: 963 # Break degeneracy of equatorial line 964 calp1 = -_TINY 965 966 # sin(alp1) * cos(bet1) = sin(alp0) 967 salp0 = salp1 * p.cbet1 968 calp0 = hypot(calp1, salp1 * p.sbet1) # calp0 > 0 969 # tan(bet1) = tan(sig1) * cos(alp1) 970 # tan(omg1) = sin(alp0) * tan(sig1) 971 # = sin(bet1) * tan(alp1) 972 somg1 = salp0 * p.sbet1 973 comg1 = calp1 * p.cbet1 974 ssig1, csig1 = norm2(p.sbet1, comg1) 975 # Without normalization we have schi1 = somg1 976 cchi1 = f1 * p.dn1 * comg1 977 978 # Enforce symmetries in the case abs(bet2) = -bet1. 979 # Need to be careful about this case, since this can 980 # yield singularities in the Newton iteration. 981 # sin(alp2) * cos(bet2) = sin(alp0) 982 salp2 = (salp0 / p.cbet2) if p.cbet2 != p.cbet1 else salp1 983 # calp2 = sqrt(1 - sq(salp2)) 984 # = sqrt(sq(calp0) - sq(sbet2)) / cbet2 985 # and subst for calp0 and rearrange to give (choose 986 # positive sqrt to give alp2 in [0, PI/2]). 987 calp2 = abs(calp1) if p.bet12 is None else ( 988 sqrt(p.bet12 + (calp1 * p.cbet1)**2) / p.cbet2) 989 # tan(bet2) = tan(sig2) * cos(alp2) 990 # tan(omg2) = sin(alp0) * tan(sig2). 991 somg2 = salp0 * p.sbet2 992 comg2 = calp2 * p.cbet2 993 ssig2, csig2 = norm2(p.sbet2, comg2) 994 # without normalization we have schi2 = somg2 995 cchi2 = f1 * p.dn2 * comg2 996 997 # omg12 = omg2 - omg1, limit to [0, PI] 998 somg12, comg12 = _sincos12(somg1, comg1, somg2, comg2, noneg=True) 999 # chi12 = chi2 - chi1, limit to [0, PI] 1000 schi12, cchi12 = _sincos12(somg1, cchi1, somg2, cchi2, noneg=True) 1001 1002 p.setsigs(ssig1, csig1, ssig2, csig2) 1003 # sig12 = sig2 - sig1, limit to [0, PI] 1004 sig12 = atan2(*_sincos12(ssig1, csig1, ssig2, csig2, noneg=True)) 1005 1006 eta12 = self._eF_reset_e2_f1_cH(calp0, salp0) / PI_2 # then ... 1007 eta12 *= fsum_(eF.deltaH(*p.sncndn2), 1008 -eF.deltaH(*p.sncndn1), sig12) 1009 # eta = chi12 - lam12 1010 lam12 = atan2(*_sincos12(p.slam12, p.clam12, schi12, cchi12)) - eta12 1011 # domg12 = chi12 - omg12 - deta12 1012 domg12 = atan2(*_sincos12( somg12, comg12, schi12, cchi12)) - eta12 1013 1014 if not diffp: 1015 dlam12 = NAN # dv > 0 in ._Newton9 1016 elif calp2: 1017 _, dlam12, _, _, _ = self._Lengths5(sig12, Caps.REDUCEDLENGTH, p) 1018 dlam12 *= f1 / (calp2 * p.cbet2) 1019 else: 1020 dlam12 = -f1 * _2_0 * p.dn1 / p.sbet1 1021 1022 # p.set_(deta12=-eta12, lam12=lam12) 1023 return lam12, sig12, salp2, calp2, domg12, dlam12 1024 1025 def _Lengths5(self, sig12, outmask, p): 1026 '''(INTERNAL) Return M{m12b = (reduced length) / self.b} and 1027 calculate M{s12b = distance / self.b} and M{m0}, the 1028 coefficient of secular term in expression for reduced length. 1029 1030 @return: 5-Tuple C{(s12b, m12b, m0, M12, M21)}. 1031 ''' 1032 s12b = m12b = m0 = M12 = M21 = NAN 1033 1034 eF = self._eF 1035 1036 # outmask &= Caps._OUTMASK 1037 if (outmask & Caps.DISTANCE): 1038 # Missing a factor of self.b 1039 s12b = eF.cE / PI_2 * fsum_(eF.deltaE(*p.sncndn2), 1040 -eF.deltaE(*p.sncndn1), sig12) 1041 1042 if (outmask & Caps._REDUCEDLENGTH_GEODESICSCALE): 1043 m0x = -eF.k2 * eF.cD / PI_2 1044 J12 = -m0x * fsum_(eF.deltaD(*p.sncndn2), 1045 -eF.deltaD(*p.sncndn1), sig12) 1046 c12 = p.csig1 * p.csig2 1047 if (outmask & Caps.REDUCEDLENGTH): 1048 m0 = m0x 1049 # Missing a factor of self.b. Add parens around 1050 # (csig1 * ssig2) and (ssig1 * csig2) to ensure 1051 # accurate cancellation for coincident points. 1052 m12b = fsum_(p.dn2 * (p.csig1 * p.ssig2), 1053 -p.dn1 * (p.ssig1 * p.csig2), J12 * c12) 1054 1055 if (outmask & Caps.GEODESICSCALE): 1056 M12 = M21 = p.ssig1 * p.ssig2 + c12 1057 t = (p.cbet1 - p.cbet2) * self.ep2 * \ 1058 (p.cbet1 + p.cbet2) / (p.dn1 + p.dn2) 1059 M12 += (p.ssig2 * t + p.csig2 * J12) * p.ssig1 / p.dn1 1060 M21 -= (p.ssig1 * t + p.csig1 * J12) * p.ssig2 / p.dn2 1061 1062 return s12b, m12b, m0, M12, M21 1063 1064 def Line(self, lat1, lon1, azi1, caps=Caps.ALL): 1065 '''Set up an L{GeodesicLineExact} to compute several points 1066 on a single geodesic. 1067 1068 @arg lat1: Latitude of the first point (C{degrees}). 1069 @arg lon1: Longitude of the first point (C{degrees}). 1070 @arg azi1: Azimuth at the first point (compass C{degrees}). 1071 @kwarg caps: Bit-or'ed combination of L{Caps} values specifying 1072 the capabilities the L{GeodesicLineExact} instance 1073 should possess, i.e., which quantities can be 1074 returnedby calls to L{GeodesicLineExact.Position} 1075 and L{GeodesicLineExact.ArcPosition}. 1076 1077 @return: A L{GeodesicLineExact} instance. 1078 1079 @note: If the point is at a pole, the azimuth is defined by keeping 1080 B{C{lon1}} fixed, writing C{B{lat1} = ±(90 − ε)}, and taking 1081 the limit C{ε → 0+}. 1082 1083 @see: C++ U{GeodesicExact.Line 1084 <https://GeographicLib.SourceForge.io/html/classGeographicLib_1_1GeodesicExact.html>} 1085 and Python U{Geodesic.Line<https://GeographicLib.SourceForge.io/html/python/code.html>}. 1086 ''' 1087 return _GeodesicLineExact(self, lat1, lon1, azi1, caps, self._debug) 1088 1089 def _Line(self, lat1, lon1, azi1, caps): # for .azimuthal._GnomonicBase.reverse 1090 '''(INTERNAL) Get a temporary L{GeodesicLineExact} instance. 1091 ''' 1092 glX = _GeodesicLineExact(self, lat1, lon1, azi1, caps, self._debug) 1093 _update_glXs(glX) # remove glX from .gxline._glXs list 1094 return glX 1095 1096 @Property_RO 1097 def n(self): 1098 '''Get the ellipsoid's I{3rd flattening} (C{float}), M{f / (2 - f) == (a - b) / (a + b)}. 1099 ''' 1100 return self.ellipsoid.n 1101 1102 @Property_RO 1103 def _n_0_1(self): 1104 '''(INTERNAL) Cached once. 1105 ''' 1106 return abs(self.n) > _0_1 1107 1108 @Property_RO 1109 def _n_6PI(self): 1110 '''(INTERNAL) Cached once. 1111 ''' 1112 return abs(self.n) * _6_0 * PI 1113 1114 def _Newton10(self, salp1, calp1, outmask, p): 1115 '''(INTERNAL) Split off from C{_GDictInverse} to reduce complexity/length. 1116 1117 @return: 10-Tuple C{(sig12, salp1, calp1, salp2, calp2, 1118 domg12, s12x, m12x, M12, M21)}. 1119 ''' 1120 # This is a straightforward solution of f(alp1) = lambda12(alp1) - 1121 # lam12 = 0 with one wrinkle. f(alp) has exactly one root in the 1122 # interval (0, PI) and its derivative is positive at the root. 1123 # Thus f(alp) is positive for alp > alp1 and negative for alp < alp1. 1124 # During the course of the iteration, a range (alp1a, alp1b) is 1125 # maintained which brackets the root and with each evaluation of 1126 # f(alp) the range is shrunk, if possible. Newton's method is 1127 # restarted whenever the derivative of f is negative (because the 1128 # new value of alp1 is then further from the solution) or if the 1129 # new estimate of alp1 lies outside (0,PI); in this case, the new 1130 # starting guess is taken to be (alp1a + alp1b) / 2. 1131 salp1a = salp1b = _TINY 1132 calp1a, calp1b = _1_0, -_1_0 1133 tripb, TOLv = False, _TOL0 1134 for it in range(_MAXIT2): 1135 # 1/4 meridian = 10e6 meter and random input, 1136 # estimated max error in nm (nano meter, by 1137 # checking Inverse problem by Direct). 1138 # 1139 # max iterations 1140 # log2(b/a) error mean sd 1141 # -7 387 5.33 3.68 1142 # -6 345 5.19 3.43 1143 # -5 269 5.00 3.05 1144 # -4 210 4.76 2.44 1145 # -3 115 4.55 1.87 1146 # -2 69 4.35 1.38 1147 # -1 36 4.05 1.03 1148 # 0 15 0.01 0.13 1149 # 1 25 5.10 1.53 1150 # 2 96 5.61 2.09 1151 # 3 318 6.02 2.74 1152 # 4 985 6.24 3.22 1153 # 5 2352 6.32 3.44 1154 # 6 6008 6.30 3.45 1155 # 7 19024 6.19 3.30 1156 v, sig12, salp2, calp2, \ 1157 domg12, dv = self._Lambda5(salp1, calp1, it < _MAXIT1, p) 1158 1159 # 2 * _TOL0 is approximately 1 ulp [0, PI] 1160 # reversed test to allow escape with NaNs 1161 if tripb or abs(v) < TOLv: 1162 break 1163 # update bracketing values 1164 if v > 0 and (it > _MAXIT1 or (calp1 / salp1) > (calp1b / salp1b)): 1165 salp1b, calp1b = salp1, calp1 1166 elif v < 0 and (it > _MAXIT1 or (calp1 / salp1) < (calp1a / salp1a)): 1167 salp1a, calp1a = salp1, calp1 1168 1169 if it < _MAXIT1 and dv > 0: 1170 dalp1 = -v / dv 1171 if abs(dalp1) < PI: 1172 s, c = sincos2(dalp1) 1173 # nalp1 = alp1 + dalp1 1174 s, c = _sincos12(-s, c, salp1, calp1) 1175 if s > 0: 1176 salp1, calp1 = norm2(s, c) 1177 # in some regimes we don't get quadratic convergence 1178 # because slope -> 0. So use convergence conditions 1179 # based on epsilon instead of sqrt(epsilon) 1180 TOLv = _TOL0 if abs(v) > _TOL016 else _TOL08 1181 continue 1182 1183 # Either dv was not positive or updated value was outside 1184 # legal range. Use the midpoint of the bracket as the next 1185 # estimate. This mechanism is not needed for the WGS84 1186 # ellipsoid, but it does catch problems with more eccentric 1187 # ellipsoids. Its efficacy is such for the WGS84 test set 1188 # with the starting guess set to alp1 = 90deg: the WGS84 1189 # test set: mean = 5.21, sd = 3.93, max = 24 WGS84 and 1190 # random input: mean = 4.74, sd = 0.99 1191 salp1, calp1 = norm2((salp1a + salp1b) * _0_5, 1192 (calp1a + calp1b) * _0_5) 1193 tripb = fsum_(calp1a, -calp1, abs(salp1a - salp1)) < _TOLb or \ 1194 fsum_(calp1b, -calp1, abs(salp1b - salp1)) < _TOLb 1195 TOLv = _TOL0 1196 1197 else: 1198 raise GeodesicError(_no_(_convergence_), _MAXIT2, txt=repr(self)) # self.toRepr() 1199 1200 s12x, m12x, _, M12, M21 = self._Lengths5(sig12, outmask, p) 1201 1202 p.set_(iter=it, trip=tripb) 1203 return (sig12, salp1, calp1, salp2, calp2, domg12, 1204 s12x, m12x, M12, M21) 1205 1206 def _sincos2f1(self, lat): 1207 '''(INTERNAL) Helper. 1208 ''' 1209 sbet, cbet = sincos2d(lat) 1210 # ensure cbet1 = +epsilon at poles; doing the fix on beta means 1211 # that sig12 will be <= 2*tiny for two points at the same pole 1212 sbet, cbet = norm2(sbet * self.f1, cbet) 1213 return sbet, max(_TINY, cbet) 1214 1215 def toStr(self, prec=6, sep=_COMMASPACE_, **unused): # PYCHOK signature 1216 '''Return this C{GeodesicExact} as string. 1217 1218 @kwarg prec: The C{float} precision, number of decimal digits (0..9). 1219 Trailing zero decimals are stripped for B{C{prec}} values 1220 of 1 and above, but kept for negative B{C{prec}} values. 1221 @kwarg sep: Optional separator to join (C{str}). 1222 1223 @return: Tuple items (C{str}). 1224 ''' 1225 d = dict(ellipsoid=self.ellipsoid, C4Order=self.C4Order) 1226 return sep.join(pairs(d, prec=prec)) 1227 1228 1229class GeodesicLineExact(_GeodesicLineExact): 1230 '''A pure Python version of I{Karney}'s C++ class U{GeodesicLineExact 1231 <https://GeographicLib.SourceForge.io/html/classGeographicLib_1_1GeodesicLineExact.html>}, 1232 modeled after I{Karney}'s Python class U{GeodesicLine<https://GeographicLib.SourceForge.io/ 1233 html/python/code.html#module-geographiclib.geodesicline>}. 1234 ''' 1235 1236 def __init__(self, geodesic, lat1, lon1, azi1, caps=Caps.STANDARD, name=NN): 1237 '''New L{GeodesicLineExact} instance, allowing points to be found along 1238 a geodesic starting at C{(B{lat1}, B{lon1})} with azimuth B{C{azi1}}. 1239 1240 @arg geodesic: The geodesic to use (L{GeodesicExact}). 1241 @arg lat1: Latitude of the first point (C{degrees}). 1242 @arg lon1: Longitude of the first point (C{degrees}). 1243 @arg azi1: Azimuth at the first points (compass C{degrees}). 1244 @kwarg caps: Bit-or'ed combination of L{Caps} values specifying 1245 the capabilities the L{GeodesicLineExact} instance 1246 should possess, i.e., which quantities can be 1247 returned by calls to L{GeodesicLineExact.Position} 1248 and L{GeodesicLineExact.ArcPosition}. 1249 @kwarg name: Optional name (C{str}). 1250 1251 @raise TypeError: Invalid B{C{geodesic}}. 1252 ''' 1253 _xinstanceof(GeodesicExact, geodesic=geodesic) 1254 1255 _GeodesicLineExact.__init__(self, geodesic, lat1, lon1, azi1, caps, 0, name=name) 1256 1257 1258def _Astroid(x, y): 1259 '''(INTERNAL) Solve M{k^4 + 2 * k^3 - (x^2 + y^2 - 1 ) * k^2 - 1260 (2 * k + 1) * y^2 = 0} for positive root k. 1261 ''' 1262 p = x**2 1263 q = y**2 1264 r = q + p - _1_0 1265 if q or r > 0: 1266 r = r / _6_0 # /= chokes PyChecker 1267 # avoid possible division by zero when r = 0 1268 # by multiplying s and t by r^3 and r, resp. 1269 S = p * q / _4_0 # S = r^3 * s 1270 r3 = r**3 1271 T3 = r3 + S 1272 # discriminant of the quadratic equation for T3 is 1273 # zero on the evolute curve p^(1/3)+q^(1/3) = 1 1274 d = S * (S + _2_0 * r3) 1275 if d < 0: 1276 # T is complex, but u is defined for a real result 1277 a = atan2(sqrt(-d), -T3) 1278 # There are three possible cube roots, choose the one which 1279 # avoids cancellation. Note d < 0 implies that r < 0. 1280 u = r * (_1_0 + _2_0 * cos(a / _3_0)) 1281 else: 1282 # pick the sign on the sqrt to maximize abs(T3) to 1283 # minimize loss of precision due to cancellation. 1284 if d: 1285 T3 += copysign0(sqrt(d), T3) # T3 = (r * t)^3 1286 # cbrt always returns the real root, cbrt(-8) = -2 1287 u = cbrt(T3) # T = r * t 1288 if u: # T can be zero; but then r2 / T -> 0 1289 u += r**2 / u 1290 u += r 1291 1292 v = sqrt(u**2 + q) 1293 # avoid loss of accuracy when u < 0 1294 uv = (q / (v - u)) if u < 0 else (u + v) 1295 w = (uv - q) / (_2_0 * v) # positive? 1296 # rearrange expression for k to avoid loss of accuracy due to 1297 # subtraction, division by 0 impossible because uv > 0, w >= 0 1298 k = uv / (sqrt(w**2 + uv) + w) # guaranteed positive 1299 1300 else: # q == 0 && r <= 0 1301 # y = 0 with |x| <= 1. Handle this case directly, for 1302 # y small, positive root is k = abs(y) / sqrt(1 - x^2) 1303 k = _0_0 1304 return k 1305 1306 1307__all__ += _ALL_DOCS(GeodesicExact, GeodesicLineExact) 1308 1309# **) MIT License 1310# 1311# Copyright (C) 2016-2021 -- mrJean1 at Gmail -- All Rights Reserved. 1312# 1313# Permission is hereby granted, free of charge, to any person obtaining a 1314# copy of this software and associated documentation files (the "Software"), 1315# to deal in the Software without restriction, including without limitation 1316# the rights to use, copy, modify, merge, publish, distribute, sublicense, 1317# and/or sell copies of the Software, and to permit persons to whom the 1318# Software is furnished to do so, subject to the following conditions: 1319# 1320# The above copyright notice and this permission notice shall be included 1321# in all copies or substantial portions of the Software. 1322# 1323# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS 1324# OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, 1325# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL 1326# THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR 1327# OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, 1328# ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR 1329# OTHER DEALINGS IN THE SOFTWARE. 1330