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

Hot-keys on this page
r m x p toggle line displays
j k next/prev highlighted chunk
0 (zero) top of page
1 (one) first highlighted chunk
# -*- coding: utf-8 -*-
<https://GeographicLib.SourceForge.io/C++/doc/classGeographicLib_1_1GeodesicExact.html>}.
Classes L{GeodesicExact} and L{GeodesicLineExact} follow the naming, methods and return values from I{Karney}s' Python classes C{Geodesic} and C{GeodesicLine}, respectively.
Copyright (C) U{Charles Karney<mailto:Charles@Karney.com>} (2008-2022) and licensed under the MIT/X11 License. For more information, see the U{GeographicLib<https://GeographicLib.SourceForge.io>} documentation. ''' # make sure int/int division yields float quotient
# A copy of comments from Karney's C{GeodesicExact.cpp}: # # This is a reformulation of the geodesic problem. The # notation is as follows: # - at a general point (no suffix or 1 or 2 as suffix) # - phi = latitude # - beta = latitude on auxiliary sphere # - omega = longitude on auxiliary sphere # - lambda = longitude # - alpha = azimuth of great circle # - sigma = arc length along great circle # - s = distance # - tau = scaled distance (= sigma at multiples of PI/2) # - at northwards equator crossing # - beta = phi = 0 # - omega = lambda = 0 # - alpha = alpha0 # - sigma = s = 0 # - a 12 suffix means a difference, e.g., s12 = s2 - s1. # - s and c prefixes mean sin and cos
_SQRT2_2, isnan, _0_0, _0_001, _0_01, _0_1, _0_5, \ _1_0, _N_1_0, _1_75, _2_0, _N_2_0, _2__PI, _3_0, \ _4_0, _6_0, _8_0, _16_0, _90_0, _180_0, _1000_0 # from pygeodesy.datums import _a_ellipsoid # from .karney _sincos12, _sin1cos2, _TINY, _xnC4 _a_ellipsoid, _EWGS84, GDict, GeodesicError, _fix90, \ _hypot, _K_2_0, _norm2, _norm180, _polynomial, \ _signBit, _sincos2, _sincos2d, _sincos2de
# increased multiplier in defn of _TOL1 from 100 to 200 to fix Inverse # case 52.784459512564 0 -52.784459512563990912 179.634407464943777557 # which otherwise failed for Visual Studio 10 (Release and Debug)
'''(INTERNAL) Return C{ang12} in C{radians}. ''' _sincos12(*sincos12)
# Using the auxiliary sphere solution with dnm computed at # (bet1 + bet2) / 2, the relative error in the azimuth # consistency check is sig12^2 * abs(f) * min(1, 1-f/2) / 2. # (Error measured for 1/100 < b/a < 100 and abs(f) >= 1/1000.
# For a given f and sig12, the max error occurs for lines # near the pole. If the old rule for computing dnm = (dn1 # + dn2)/2 is used, then the error increases by a factor of # 2.) Setting this equal to epsilon gives sig12 = etol2.
# Here 0.1 is a safety factor (error decreased by 100) and # max(0.001, abs(f)) stops etol2 getting too large in the # nearly spherical case.
'''(INTERNAL) Parameters passed around in C{._GDictInverse} and optionally returned when C{GeodesicExact.debug} is C{True}. ''' '''Update the C{sig1} and C{sig2} parameters. ''' ssig2=ssig2, csig2=csig2, sncndn2=(ssig2, csig2, self.dn2)) # PYCHOK dn2
def toGDict(self): # PYCHOK no cover '''Return as C{GDict} without attrs C{sncndn1} and C{sncndn2}. ''' def _rest(sncndn1=None, sncndn2=None, **rest): # PYCHOK sncndn* not used return GDict(rest)
return _rest(**self)
'''A pure Python version of I{Karney}'s C++ class U{GeodesicExact <https://GeographicLib.SourceForge.io/C++/doc/classGeographicLib_1_1GeodesicExact.html>}, modeled after I{Karney}'s Python class U{Geodesic<https://GeographicLib.SourceForge.io/ html/python/code.html#module-geographiclib.geodesic>}. '''
C4Order=None): # for backward compatibility '''New L{GeodesicExact} instance.
@arg a_ellipsoid: An ellipsoid (L{Ellipsoid}) or datum (L{Datum}) or the equatorial radius of the ellipsoid (C{scalar}, conventionally in C{meter}), see B{C{f}}. @arg f: The flattening of the ellipsoid (C{scalar}) if B{C{a_ellipsoid}} is specified as C{scalar}. @kwarg name: Optional name (C{str}). @kwarg C4order: Optional series expansion order (C{int}), see property L{C4order}, default C{30}. @kwarg C4Order: DEPRECATED, use keyword argument B{C{C4order}}.
@raise GeodesicError: Invalid B{C{C4order}}. '''
self.C4order = C4order self.C4Order = C4Order
'''Get the I{equatorial} radius, semi-axis (C{meter}). '''
'''Solve the I{Direct} geodesic problem in terms of (spherical) arc length.
@arg lat1: Latitude of the first point (C{degrees}). @arg lon1: Longitude of the first point (C{degrees}). @arg azi1: Azimuth at the first point (compass C{degrees}). @arg a12: Arc length between the points (C{degrees}), can be negative. @kwarg outmask: Bit-or'ed combination of L{Caps} values specifying the quantities to be returned.
@return: A L{GDict} with up to 12 items C{lat1, lon1, azi1, lat2, lon2, azi2, m12, a12, s12, M12, M21, S12} with C{lat1}, C{lon1}, C{azi1} and arc length C{a12} always included.
@see: C++ U{GeodesicExact.ArcDirect <https://GeographicLib.SourceForge.io/C++/doc/classGeographicLib_1_1GeodesicExact.html>} and Python U{Geodesic.ArcDirect<https://GeographicLib.SourceForge.io/C++/doc/python/code.html>}. '''
'''Define a L{GeodesicLineExact} in terms of the I{direct} geodesic problem and as arc length.
@arg lat1: Latitude of the first point (C{degrees}). @arg lon1: Longitude of the first point (C{degrees}). @arg azi1: Azimuth at the first point (compass C{degrees}). @arg a12: Arc length between the points (C{degrees}), can be negative. @kwarg caps: Bit-or'ed combination of L{Caps} values specifying the capabilities the L{GeodesicLineExact} instance should possess, i.e., which quantities can be returned by calls to L{GeodesicLineExact.Position} and L{GeodesicLineExact.ArcPosition}.
@return: A L{GeodesicLineExact} instance.
@note: The third point of the L{GeodesicLineExact} is set to correspond to the second point of the I{Inverse} geodesic problem.
@note: Latitude B{C{lat1}} should in the range C{[-90, +90]}.
@see: C++ U{GeodesicExact.ArcDirectLine <https://GeographicLib.SourceForge.io/C++/doc/classGeographicLib_1_1GeodesicExact.html>} and Python U{Geodesic.ArcDirectLine<https://GeographicLib.SourceForge.io/C++/doc/python/code.html>}. '''
'''Set up a L{GeodesicAreaExact} to compute area and perimeter of a polygon.
@kwarg polyline: If C{True} perimeter only, otherwise area and perimeter (C{bool}). @kwarg name: Optional name (C{str}).
@return: A L{GeodesicAreaExact} instance.
@note: The B{C{debug}} setting is passed as C{verbose} to the returned L{GeodesicAreaExact} instance. ''' name=name or self.name) gaX.verbose = True
'''Get the ellipsoid's I{polar} radius, semi-axis (C{meter}). '''
'''Get the ellipsoid's I{authalic} earth radius I{squared} (C{meter} I{squared}). ''' # The Geodesic class substitutes atanh(sqrt(e2)) for asinh(sqrt(ep2)) # in the definition of _c2. The latter is more accurate for very # oblate ellipsoids (which the Geodesic class does not handle). Of # course, the area calculation in GeodesicExact is still based on a # series and only holds for moderately oblate (or prolate) ellipsoids.
'''Evaluate the C{C4x} coefficients for B{C{eps}}.
@arg eps: Polynomial factor (C{float}).
@return: C{C4order}-Tuple of C{C4x(B{eps})} coefficients. ''' # assert i == (nC4 * (nC4 + 1)) // 2
'''(INTERNAL) Compute C{eps} from B{C{k2}} and invoke C{C4f}. '''
'''Get the series expansion order (C{int}, 24, 27 or 30). '''
'''Set the series expansion order (C{int}, 24, 27 or 30).
@raise GeodesicError: Invalid B{C{order}}. '''
'''DEPRECATED, use property C{C4order}. ''' return self.C4order
'''DEPRECATED, use property C{C4order}. '''
'''(INTERNAL) Get the C{C4_24}, C{_27} or C{_30} series coefficients. ''' except (AttributeError, ImportError, TypeError) as x: raise GeodesicError(nC4=nC4, txt=str(x)) raise GeodesicError(_coeffs=len(_coeffs), _xnC4=n, nC4=nC4)
'''Get this ellipsoid's C{C4} coefficients, I{cached} tuple.
@see: Property L{C4order}. ''' # see C4coeff() in GeographicLib.src.GeodesicExactC4.cpp # assert i == (nC4 * (nC4 + 1) * (nC4 + 5)) // 6
'''Solve the I{Direct} geodesic problem
@arg lat1: Latitude of the first point (C{degrees}). @arg lon1: Longitude of the first point (C{degrees}). @arg azi1: Azimuth at the first point (compass C{degrees}). @arg s12: Distance between the points (C{meter}), can be negative. @kwarg outmask: Bit-or'ed combination of L{Caps} values specifying the quantities to be returned.
@return: A L{GDict} with up to 12 items C{lat1, lon1, azi1, lat2, lon2, azi2, m12, a12, s12, M12, M21, S12} with C{lat1}, C{lon1}, C{azi1} and distance C{s12} always included.
@see: C++ U{GeodesicExact.Direct <https://GeographicLib.SourceForge.io/C++/doc/classGeographicLib_1_1GeodesicExact.html>} and Python U{Geodesic.Direct<https://GeographicLib.SourceForge.io/C++/doc/python/code.html>}. '''
'''Return the destination lat, lon and reverse azimuth (final bearing) in C{degrees}.
@return: L{Destination3Tuple}C{(lat, lon, final)}. '''
'''Define a L{GeodesicLineExact} in terms of the I{direct} geodesic problem and as distance.
@arg lat1: Latitude of the first point (C{degrees}). @arg lon1: Longitude of the first point (C{degrees}). @arg azi1: Azimuth at the first point (compass C{degrees}). @arg s12: Distance between the points (C{meter}), can be negative. @kwarg caps: Bit-or'ed combination of L{Caps} values specifying the capabilities the L{GeodesicLineExact} instance should possess, i.e., which quantities can be returned by calls to L{GeodesicLineExact.Position}.
@return: A L{GeodesicLineExact} instance.
@note: The third point of the L{GeodesicLineExact} is set to correspond to the second point of the I{Inverse} geodesic problem.
@note: Latitude B{C{lat1}} should in the range C{[-90, +90]}.
@see: C++ U{GeodesicExact.DirectLine <https://GeographicLib.SourceForge.io/C++/doc/classGeographicLib_1_1GeodesicExact.html>} and Python U{Geodesic.DirectLine<https://GeographicLib.SourceForge.io/C++/doc/python/code.html>}. '''
'''(INTERNAL) Helper. ''' if self.f < 0: # PYCHOK no cover dn = sqrt(_1_0 - cbet**2 * self.e2) / self.f1 else:
'''Get the ellipsoid's I{(1st) eccentricity squared} (C{float}), M{f * (2 - f)}. '''
'''(INTERNAL) Cache M{E.e2 * E.a2}. '''
'''(INTERNAL) Cache M{E.e2 * E.f1}. '''
'''(INTERNAL) Get the elliptic function, aka C{.E}. '''
'''(INTERNAL) Reset elliptic function and return M{cH * e2 / f1 * ...}. '''
'''(INTERNAL) Reset elliptic function and return C{k2}. '''
'''Get the ellipsoid (C{Ellipsoid}). '''
'''Get the ellipsoid's I{2nd eccentricity squared} (C{float}), M{e2 / (1 - e2)}. '''
'''(INTERNAL) The si12 threshold for "really short". '''
'''Get the ellipsoid's I{flattening} (C{float}), M{(a - b) / a}, C{0} for spherical, negative for prolate. '''
'''Get the ellipsoid's I{1 - flattening} (C{float}). '''
'''(INTERNAL) Cached/memoized. '''
'''(INTERNAL) As C{_GenDirect}, but returning a L{GDict}.
@return: A L{GDict} ... '''
'''(INTERNAL) As C{_GenInverse}, but returning a L{GDict}.
@return: A L{GDict} ... ''' if self._debug: # PYCHOK no cover outmask |= Caps._DEBUG_INVERSE & self._debug # compute longitude difference carefully (with _diff182): # result is in [-180, +180] but -180 is only for west-going # geodesics, +180 is for east-going and meridional geodesics # see C{result} from geographiclib.geodesic.Inverse else: # make longitude difference positive # calculate sincosd(_around(lon12 + correction)) # supplementary longitude difference else: # GeographicLib 1.52 # make longitude difference positive and if very close # to being on the same half-meridian, then make it so. if lon12 < 0: # _signBit(lon12) lon_, lon12 = True, -_around(lon12) lon12s = _around(fsum_(_180_0, -lon12, lon12s)) else: lon_, lon12 = False, _around(lon12) lon12s = _around(fsum_(_180_0, -lon12, -lon12s)) lam12 = radians(lon12) if lon12 > _90_0: slam12, clam12 = _sincos2d(lon12s) clam12 = -clam12 else: slam12, clam12 = _sincos2(lam12) # If really close to the equator, treat as on equator. # Swap points so that point with higher (abs) latitude is # point 1. If one latitude is a NAN, then it becomes lat1. else: # make lat1 <= -0 # Now 0 <= lon12 <= 180, -90 <= lat1 <= -0 and lat1 <= lat2 <= -lat1 # and lat_, lon_, swap_ register the transformation to bring the # coordinates to this canonical form, where False means no change # made. We make these transformations so that there are few cases # to check, e.g., on verifying quadrants in atan2. In addition, # this enforces some symmetries in the results returned.
# Initialize for the meridian. No longitude calculation is # done in this case to let the parameter default to 0. # If cbet1 < -sbet1, then cbet2 - cbet1 is a sensitive measure # of the |bet1| - |bet2|. Alternatively (cbet1 >= -sbet1), # abs(sbet2) + sbet1 is a better measure. This logic is used # in assigning calp2 in _Lambda6. Sometimes these quantities # vanish and in that case we force bet2 = +/- bet1 exactly. An # example where is is necessary is the inverse problem # 48.522876735459 0 -48.52287673545898293 179.599720456223079643 # which failed with Visual Studio 10 (Release and Debug)
sbet2=sbet2, cbet2=cbet2, dn2=self._dn(sbet2, cbet2))
# Endpoints are on a single full meridian, # so the geodesic might lie on a meridian. # tan(bet) = tan(sig) * cos(alp) # sig12 = sig2 - sig1 M12, M21 = self._Length5(sig12, outmask | Caps.REDUCEDLENGTH, p) # Add the check for sig12 since zero length geodesics # might yield m12 < 0. Test case was # echo 20.001 0 20.001 0 | GeodSolve -i # In fact, we will have sig12 > PI/2 for meridional # geodesic which is not a shortest path. # Need at least 2 to handle 90 0 90 180 # Prevent negative s12 or m12 from geographiclib 1.52 else: # else: # m12 < 0, prolate and too close to anti-podal # _meridian = True
else: # Now point1 and point2 belong within a hemisphere bounded by a # meridian and geodesic is neither meridional or equatorial. # Figure a starting point for Newton's method salp2, calp2, dnm = self._InverseStart6(lam12, p) # pre-compute the constant _Lambda6 term, once (((cbet1 + cbet2) * (cbet2 - cbet1)) if cbet1 < -sbet1 else ((sbet1 + sbet2) * (sbet1 - sbet2)))) salp2, calp2, domg12 = self._Newton6(salp1, calp1, p) # omg12 = lam12 - domg12 else: # from _InverseStart6: dnm, salp*, calp*
else: # _meridian is False
M21=unsigned0(M21))
salp2, calp2, somg12, comg12, p)
azi2=_atan2d_reverse(salp2, calp2, reverse=outmask & Caps.REVERSE2)) salp2=salp2, calp2=calp2)
if (outmask & Caps._DEBUG_INVERSE): # PYCHOK no cover E, eF = self.ellipsoid, self._eF p.set_(C=C, a=self.a, f=self.f, f1=self.f1, e=E.e, e2=self.e2, ep2=self.ep2, c2=E.c2, c2x=self.c2x, eFcD=eF.cD, eFcE=eF.cE, eFcH=eF.cH, eFk2=eF.k2, eFa2=eF.alpha2) p.update(r) # r overrides p r = p.toGDict()
'''(INTERNAL) The general I{Inverse} geodesic calculation.
@return: L{Direct9Tuple}C{(a12, lat2, lon2, azi2, s12, m12, M12, M21, S12)}. ''' r = self._GDictDirect(lat1, lon1, azi1, arcmode, s12_a12, outmask) return r.toDirect9Tuple()
'''(INTERNAL) Helper for C{ArcDirectLine} and C{DirectLine}.
@return: A L{GeodesicLineExact} instance. ''' # guard against underflow in salp0. Also -0 is converted to +0. self._debug, s, c)._GenSet(arcmode, s12_a12)
'''(INTERNAL) The general I{Inverse} geodesic calculation.
@return: L{Inverse10Tuple}C{(a12, s12, salp1, calp1, salp2, calp2, m12, M12, M21, S12)}. ''' r = self._GDictInverse(lat1, lon1, lat2, lon2, outmask | Caps._SALPs_CALPs) return r.toInverse10Tuple()
'''Perform the I{Inverse} geodesic calculation.
@arg lat1: Latitude of the first point (C{degrees}). @arg lon1: Longitude of the first point (C{degrees}). @arg lat2: Latitude of the second point (C{degrees}). @arg lon2: Longitude of the second point (C{degrees}). @kwarg outmask: Bit-or'ed combination of L{Caps} values specifying the quantities to be returned.
@return: A L{GDict} with up to 12 items C{lat1, lon1, azi1, lat2, lon2, azi2, m12, a12, s12, M12, M21, S12} with C{lat1}, C{lon1}, C{azi1} and distance C{s12} always included.
@note: The third point of the L{GeodesicLineExact} is set to correspond to the second point of the I{Inverse} geodesic problem.
@note: Both B{C{lat1}} and B{C{lat2}} should in the range C{[-90, +90]}.
@see: C++ U{GeodesicExact.InverseLine <https://GeographicLib.SourceForge.io/C++/doc/classGeographicLib_1_1GeodesicExact.html>} and Python U{Geodesic.InverseLine<https://GeographicLib.SourceForge.io/C++/doc/python/code.html>}. '''
'''Return the non-negative, I{angular} distance in C{degrees}. ''' # see .FrechetKarney.distance, .HausdorffKarney._distance # and .HeightIDWkarney._distances
'''Return the distance in C{meter} and the forward and reverse azimuths (initial and final bearing) in C{degrees}.
@return: L{Distance3Tuple}C{(distance, initial, final)}. ''' iteration=r.iteration)
'''Define a L{GeodesicLineExact} in terms of the I{Inverse} geodesic problem.
@arg lat1: Latitude of the first point (C{degrees}). @arg lon1: Longitude of the first point (C{degrees}). @arg lat2: Latitude of the second point (C{degrees}). @arg lon2: Longitude of the second point (C{degrees}). @kwarg caps: Bit-or'ed combination of L{Caps} values specifying the capabilities the L{GeodesicLineExact} instance should possess, i.e., which quantities can be returned by calls to L{GeodesicLineExact.Position} and L{GeodesicLineExact.ArcPosition}.
@return: A L{GeodesicLineExact} instance.
@note: The third point of the L{GeodesicLineExact} is set to correspond to the second point of the I{Inverse} geodesic problem.
@note: Both B{C{lat1}} and B{C{lat2}} should in the range C{[-90, +90]}.
@see: C++ U{GeodesicExact.InverseLine <https://GeographicLib.SourceForge.io/C++/doc/classGeographicLib_1_1GeodesicExact.html>} and Python U{Geodesic.InverseLine<https://GeographicLib.SourceForge.io/C++/doc/python/code.html>}. ''' self._debug, r.salp1, r.calp1)._GenSet(True, r.a12)
salp2, calp2, somg12, comg12, p): '''(INTERNAL) Split off from C{_GDictInverse} to reduce complexity/length.
@return: Area C{S12}. ''' # from _Lambda6: sin(alp1) * cos(bet1) = sin(alp0), calp0 > 0 # from _Lambda6: tan(bet) = tan(sig) * cos(alp) # multiplier = a^2 * e^2 * cos(alpha0) * sin(alpha0) else: # avoid problems with indeterminate sig1, sig2 on equator
comg12 > -_SQRT2_2 and # lon diff not too big (p.sbet2 - p.sbet1) < _1_75): # lat diff not too big # use tan(Gamma/2) = tan(omg12/2) * # (tan(bet1/2) + tan(bet2/2)) / # (tan(bet1/2) * tan(bet2/2) + 1)) # with tan(x/2) = sin(x) / (1 + cos(x)) else: # alp12 = alp2 - alp1, used in atan2, no need to normalize # The right thing appears to happen if alp1 = +/-180 and # alp2 = 0, viz salp12 = -0 and alp12 = -180. However, # this depends on the sign being attached to 0 correctly. # Following ensures the correct behavior. alp12 = _copysign(PI, calp1) else:
'''(INTERNAL) Return a starting point for Newton's method in C{salp1} and C{calp1} indicated by C{sig12=None}. If Newton's method doesn't need to be used, return also C{salp2}, C{calp2}, C{dnm} and C{sig12} non-C{None}.
@return: 6-Tuple C{(sig12, salp1, calp1, salp2, calp2, dnm)} and C{p.setsigs} updated for Newton, C{sig12=None}. '''
# bet12 = bet2 - bet1 in [0, PI) # sin((bet1 + bet2)/2)^2 = (sbet1 + sbet2)^2 / ( # (cbet1 + cbet2)^2 + (sbet1 + sbet2)^2) else:
# bet12a = bet2 + bet1 in (-PI, 0], note -sbet1
sbet12 - p.cbet1 * p.sbet2 * t)
csig12 >= 0 or ssig12 >= (p.cbet1**2 * self._n6PI)):
else: # Scale lam12 and bet2 to x, y coordinate system where antipodal # point is at origin and singular point is at y = 0, x = -1 if f < 0: # PYCHOK no cover # ssig1=sbet1, csig1=-cbet1, ssig2=sbet2, csig2=cbet2 p.setsigs(p.sbet1, -p.cbet1, p.sbet2, p.cbet2) # if lon12 = 180, this repeats a calculation made in Inverse _, m12b, m0, _, _ = self._Length5(atan2(sbet12a, cbet12a) + PI, Caps.REDUCEDLENGTH, p) t = p.cbet1 * PI # x = dlat, y = dlon x = m12b / (t * p.cbet2 * m0) - _1_0 sca = (sbet12a / (x * p.cbet1)) if x < -_0_01 else (-f * t) y = lam12x / sca else: # f >= 0, however f == 0 does not get here
if f < 0: # PYCHOK no cover calp1 = max(_0_0, x) if x > _TOL1 else max(_N_1_0, x) salp1 = sqrt(_1_0 - calp1**2) else: salp1 = min(_1_0, -x) calp1 = -sqrt(_1_0 - salp1**2) else: # Estimate alp1, by solving the astroid problem. # # Could estimate alpha1 = theta + PI/2, directly, i.e., # calp1 = y/k; salp1 = -x/(1+k); for _f >= 0 # calp1 = x/(1+k); salp1 = -y/k; for _f < 0 (need to check) # # However, it's better to estimate omg12 from astroid and use # spherical formula to compute alp1. This reduces the mean # number of Newton iterations for astroid cases from 2.24 # (min 0, max 6) to 2.12 (min 0, max 5). The changes in the # number of iterations are as follows: # # change percent # 1 5 # 0 78 # -1 16 # -2 0.6 # -3 0.04 # -4 0.002 # # The histogram of iterations is (m = number of iterations # estimating alp1 directly, n = number of iterations # estimating via omg12, total number of trials = 148605): # # iter m n # 0 148 186 # 1 13046 13845 # 2 93315 102225 # 3 36189 32341 # 4 5396 7 # 5 455 1 # 6 56 0 # # omg12 is near PI, estimate work with omg12a = PI - omg12 (x * k / (k + _1_0)) # update spherical estimate of alp1 using omg12 instead of lam12
# sanity check on starting guess. Backwards check allows NaN through.
'''(INTERNAL) Helper.
@return: 6-Tuple C{(lam12, sig12, salp2, calp2, domg12, dlam12} and C{p.setsigs} updated. '''
if p.sbet1 == 0 and calp1 == 0: # PYCHOK no cover # Break degeneracy of equatorial line calp1 = -_TINY
# sin(alp1) * cos(bet1) = sin(alp0), # calp0 > 0 # tan(bet1) = tan(sig1) * cos(alp1) # tan(omg1) = sin(alp0) * tan(sig1) # = sin(bet1) * tan(alp1) # Without normalization we have schi1 = somg1
# Enforce symmetries in the case abs(bet2) = -bet1. # Need to be careful about this case, since this can # yield singularities in the Newton iteration. # sin(alp2) * cos(bet2) = sin(alp0) # calp2 = sqrt(1 - sq(salp2)) # = sqrt(sq(calp0) - sq(sbet2)) / cbet2 # and subst for calp0 and rearrange to give (choose # positive sqrt to give alp2 in [0, PI/2]). sqrt((calp1 * p.cbet1)**2 + p.bet12) / p.cbet2) # tan(bet2) = tan(sig2) * cos(alp2) # tan(omg2) = sin(alp0) * tan(sig2). # without normalization we have schi2 = somg2
# omg12 = omg2 - omg1, limit to [0, PI] # chi12 = chi2 - chi1, limit to [0, PI]
# sig12 = sig2 - sig1, limit to [0, PI]
-eF.deltaH(*p.sncndn1), sig12) # eta = chi12 - lam12 # domg12 = chi12 - omg12 - deta12
elif p.sbet1: dlam12 = -f1 * p.dn1 * _2_0 / p.sbet1
# p.set_(deta12=-eta12, lam12=lam12)
'''(INTERNAL) Return M{m12b = (reduced length) / self.b} and calculate M{s12b = distance / self.b} and M{m0}, the coefficient of secular term in expression for reduced length and the geodesic scales C{M12} and C{M21}.
@return: 5-Tuple C{(s12b, m12b, m0, M12, M21)}. '''
# outmask &= Caps._OUT_MASK # Missing a factor of self.b -eF.deltaE(*p.sncndn1), sig12)
-eF.deltaD(*p.sncndn1), sig12) # Missing a factor of self.b. Add parens around # (csig1 * ssig2) and (ssig1 * csig2) to ensure # accurate cancellation for coincident points. -p.dn1 * (p.ssig1 * p.csig2), J12 * (p.csig1 * p.csig2)) p.csig1 * p.csig2 (p.cbet1 + p.cbet2) / (p.dn1 + p.dn2)
'''Set up a L{GeodesicLineExact} to compute several points on a single geodesic.
@arg lat1: Latitude of the first point (C{degrees}). @arg lon1: Longitude of the first point (C{degrees}). @arg azi1: Azimuth at the first point (compass C{degrees}). @kwarg caps: Bit-or'ed combination of L{Caps} values specifying the capabilities the L{GeodesicLineExact} instance should possess, i.e., which quantities can be returnedby calls to L{GeodesicLineExact.Position} and L{GeodesicLineExact.ArcPosition}.
@return: A L{GeodesicLineExact} instance.
@note: If the point is at a pole, the azimuth is defined by keeping B{C{lon1}} fixed, writing C{B{lat1} = ±(90 − ε)}, and taking the limit C{ε → 0+}.
@see: C++ U{GeodesicExact.Line <https://GeographicLib.SourceForge.io/C++/doc/classGeographicLib_1_1GeodesicExact.html>} and Python U{Geodesic.Line<https://GeographicLib.SourceForge.io/C++/doc/python/code.html>}. '''
'''Get the ellipsoid's I{3rd flattening} (C{float}), M{f / (2 - f) == (a - b) / (a + b)}. '''
'''(INTERNAL) Cached once. '''
'''(INTERNAL) Cached once. '''
'''(INTERNAL) Split off from C{_GDictInverse} to reduce complexity/length.
@return: 6-Tuple C{(sig12, salp1, calp1, salp2, calp2, domg12)} and C{p.iter} and C{p.trip} updated. ''' # This is a straightforward solution of f(alp1) = lambda12(alp1) - # lam12 = 0 with one wrinkle. f(alp) has exactly one root in the # interval (0, PI) and its derivative is positive at the root. # Thus f(alp) is positive for alp > alp1 and negative for alp < alp1. # During the course of the iteration, a range (alp1a, alp1b) is # maintained which brackets the root and with each evaluation of # f(alp) the range is shrunk, if possible. Newton's method is # restarted whenever the derivative of f is negative (because the # new value of alp1 is then further from the solution) or if the # new estimate of alp1 lies outside (0,PI); in this case, the new # starting guess is taken to be (alp1a + alp1b) / 2. # 1/4 meridian = 10e6 meter and random input, # estimated max error in nm (nano meter, by # checking Inverse problem by Direct). # # max iterations # log2(b/a) error mean sd # -7 387 5.33 3.68 # -6 345 5.19 3.43 # -5 269 5.00 3.05 # -4 210 4.76 2.44 # -3 115 4.55 1.87 # -2 69 4.35 1.38 # -1 36 4.05 1.03 # 0 15 0.01 0.13 # 1 25 5.10 1.53 # 2 96 5.61 2.09 # 3 318 6.02 2.74 # 4 985 6.24 3.22 # 5 2352 6.32 3.44 # 6 6008 6.30 3.45 # 7 19024 6.19 3.30 domg12, dv = self._Lambda6(salp1, calp1, i < _MAXIT1, p)
# 2 * _TOL0 is approximately 1 ulp [0, PI] # reversed test to allow escape with NaNs # update bracketing values
# nalp1 = alp1 + dalp1 # in some regimes we don't get quadratic convergence # because slope -> 0. So use convergence conditions # based on epsilon instead of sqrt(epsilon)
# Either dv was not positive or updated value was outside # legal range. Use the midpoint of the bracket as the next # estimate. This mechanism is not needed for the WGS84 # ellipsoid, but it does catch problems with more eccentric # ellipsoids. Its efficacy is such for the WGS84 test set # with the starting guess set to alp1 = 90 deg: the WGS84 # test set: mean = 5.21, sd = 3.93, max = 24 and WGS84 with # random input: mean = 4.74, sd = 0.99 salp1, calp1 = _norm2((salp1a + salp1b) * _0_5, (calp1a + calp1b) * _0_5) tripb = fsum1_(calp1a, -calp1, fabs(salp1a - salp1)) < _TOLb or \ fsum1_(calp1b, -calp1, fabs(salp1b - salp1)) < _TOLb TOLv = _TOL0
else: raise GeodesicError(Fmt.no_convergence(v, TOLv), txt=repr(self)) # self.toRepr()
'''(INTERNAL) Helper, also for C{_G_GeodesicLineExact}. ''' # ensure cbet1 = +epsilon at poles; doing the fix on beta means # that sig12 will be <= 2*tiny for two points at the same pole
'''Return this C{GeodesicExact} as string.
@kwarg prec: The C{float} precision, number of decimal digits (0..9). Trailing zero decimals are stripped for B{C{prec}} values of 1 and above, but kept for negative B{C{prec}} values. @kwarg sep: Separator to join (C{str}).
@return: Tuple items (C{str}). '''
'''A pure Python version of I{Karney}'s C++ class U{GeodesicLineExact <https://GeographicLib.SourceForge.io/C++/doc/classGeographicLib_1_1GeodesicLineExact.html>}, modeled after I{Karney}'s Python class U{GeodesicLine<https://GeographicLib.SourceForge.io/ html/python/code.html#module-geographiclib.geodesicline>}. '''
'''New L{GeodesicLineExact} instance, allowing points to be found along a geodesic starting at C{(B{lat1}, B{lon1})} with azimuth B{C{azi1}}.
@arg geodesic: The geodesic to use (L{GeodesicExact}). @arg lat1: Latitude of the first point (C{degrees}). @arg lon1: Longitude of the first point (C{degrees}). @arg azi1: Azimuth at the first points (compass C{degrees}). @kwarg caps: Bit-or'ed combination of L{Caps} values specifying the capabilities the L{GeodesicLineExact} instance should possess, i.e., which quantities can be returned by calls to L{GeodesicLineExact.Position} and L{GeodesicLineExact.ArcPosition}. @kwarg name: Optional name (C{str}).
@raise TypeError: Invalid B{C{geodesic}}. ''' geodesic = geodesic.copy(deep=False, name=NN(_UNDER_, geodesic.name)) # _update_all(geodesic)
'''(INTERNAL) Solve M{k^4 + 2 * k^3 - (x^2 + y^2 - 1 ) * k^2 - (2 * k + 1) * y^2 = 0} for positive root k. ''' # avoid possible division by zero when r = 0 # by multiplying s and t by r^3 and r, resp. # discriminant of the quadratic equation for T3 is # zero on the evolute curve p^(1/3) + q^(1/3) = 1 # T is complex, but u is defined for a real result a = atan2(sqrt(-d), -T3) / _3_0 # There are 3 possible cube roots, choose the one which # avoids cancellation. Note d < 0 implies that r < 0. u = (cos(a) * _2_0 + _1_0) * r else: # pick the sign on the sqrt to maximize abs(T3) to # minimize loss of precision due to cancellation. # _cbrt always returns the real root, _cbrt(-8) = -2 # avoid loss of accuracy when u < 0 # rearrange expression for k to avoid loss of accuracy due to # subtraction, division by 0 impossible because u > 0, w >= 0
else: # q == 0 && r <= 0 # y = 0 with |x| <= 1. Handle this case directly, for # y small, positive root is k = abs(y) / sqrt(1 - x^2) k = _0_0
# **) MIT License # # Copyright (C) 2016-2022 -- mrJean1 at Gmail -- All Rights Reserved. # # Permission is hereby granted, free of charge, to any person obtaining a # copy of this software and associated documentation files (the "Software"), # to deal in the Software without restriction, including without limitation # the rights to use, copy, modify, merge, publish, distribute, sublicense, # and/or sell copies of the Software, and to permit persons to whom the # Software is furnished to do so, subject to the following conditions: # # The above copyright notice and this permission notice shall be included # in all copies or substantial portions of the Software. # # THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS # OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, # FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL # THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR # OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, # ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR # OTHER DEALINGS IN THE SOFTWARE. |