Secure your code as it's written. Use Snyk Code to scan source code in minutes - no build needed - and fix issues immediately.
"""Solve geodesic problems"""
GEOGRAPHICLIB_GEODESIC_ORDER = 6
nA1_ = GEOGRAPHICLIB_GEODESIC_ORDER
nC1_ = GEOGRAPHICLIB_GEODESIC_ORDER
nC1p_ = GEOGRAPHICLIB_GEODESIC_ORDER
nA2_ = GEOGRAPHICLIB_GEODESIC_ORDER
nC2_ = GEOGRAPHICLIB_GEODESIC_ORDER
nA3_ = GEOGRAPHICLIB_GEODESIC_ORDER
nA3x_ = nA3_
nC3_ = GEOGRAPHICLIB_GEODESIC_ORDER
nC3x_ = (nC3_ * (nC3_ - 1)) // 2
nC4_ = GEOGRAPHICLIB_GEODESIC_ORDER
nC4x_ = (nC4_ * (nC4_ + 1)) // 2
maxit1_ = 20
maxit2_ = maxit1_ + Math.digits + 10
tiny_ = math.sqrt(Math.minval)
tol0_ = Math.epsilon
tol1_ = 200 * tol0_
tol2_ = math.sqrt(tol0_)
tolb_ = tol0_ * tol2_
xthresh_ = 1000 * tol2_
CAP_NONE = GeodesicCapability.CAP_NONE
CAP_C1 = GeodesicCapability.CAP_C1
CAP_C1p = GeodesicCapability.CAP_C1p
CAP_C2 = GeodesicCapability.CAP_C2
CAP_C3 = GeodesicCapability.CAP_C3
CAP_C4 = GeodesicCapability.CAP_C4
CAP_ALL = GeodesicCapability.CAP_ALL
CAP_MASK = GeodesicCapability.CAP_MASK
E = Math.copysign(1, self._salp0) # East or west going?
# omg12 = omg2 - omg1
omg12 = (E * (sig12
- (math.atan2( ssig2, csig2) -
math.atan2( self._ssig1, self._csig1))
+ (math.atan2(E * somg2, comg2) -
math.atan2(E * self._somg1, self._comg1)))
if outmask & Geodesic.LONG_UNROLL
else math.atan2(somg2 * self._comg1 - comg2 * self._somg1,
comg2 * self._comg1 + somg2 * self._somg1))
lam12 = omg12 + self._A3c * (
sig12 + (Geodesic._SinCosSeries(True, ssig2, csig2, self._C3a)
- self._B31))
lon12 = math.degrees(lam12)
lon2 = (self.lon1 + lon12 if outmask & Geodesic.LONG_UNROLL else
Math.AngNormalize(Math.AngNormalize(self.lon1) +
Math.AngNormalize(lon12)))
if outmask & Geodesic.LATITUDE:
lat2 = Math.atan2d(sbet2, self._f1 * cbet2)
if outmask & Geodesic.AZIMUTH:
azi2 = Math.atan2d(salp2, calp2)
if outmask & (Geodesic.REDUCEDLENGTH | Geodesic.GEODESICSCALE):
B22 = Geodesic._SinCosSeries(True, ssig2, csig2, self._C2a)
AB2 = (1 + self._A2m1) * (B22 - self._B21)
J12 = (self._A1m1 - self._A2m1) * sig12 + (AB1 - AB2)
if outmask & Geodesic.REDUCEDLENGTH:
# Add parens around (_csig1 * ssig2) and (_ssig1 * csig2) to ensure
# accurate cancellation in the case of coincident points.
m12 = self._b * (( dn2 * (self._csig1 * ssig2) -
# 1/5 157e6 3.8e9 280e6
ssig2 = self._ssig1 * csig12 + self._csig1 * ssig12
csig2 = self._csig1 * csig12 - self._ssig1 * ssig12
B12 = Geodesic._SinCosSeries(True, ssig2, csig2, self._C1a)
serr = ((1 + self._A1m1) * (sig12 + (B12 - self._B11)) -
s12_a12 / self._b)
sig12 = sig12 - serr / math.sqrt(1 + self._k2 * Math.sq(ssig2))
ssig12 = math.sin(sig12); csig12 = math.cos(sig12)
# Update B12 below
# real omg12, lam12, lon12
# real ssig2, csig2, sbet2, cbet2, somg2, comg2, salp2, calp2
# sig2 = sig1 + sig12
ssig2 = self._ssig1 * csig12 + self._csig1 * ssig12
csig2 = self._csig1 * csig12 - self._ssig1 * ssig12
dn2 = math.sqrt(1 + self._k2 * Math.sq(ssig2))
if outmask & (
Geodesic.DISTANCE | Geodesic.REDUCEDLENGTH | Geodesic.GEODESICSCALE):
if arcmode or abs(self.f) > 0.01:
B12 = Geodesic._SinCosSeries(True, ssig2, csig2, self._C1a)
AB1 = (1 + self._A1m1) * (B12 - self._B11)
# sin(bet2) = cos(alp0) * sin(sig2)
sbet2 = self._calp0 * ssig2
# Alt: cbet2 = hypot(csig2, salp0 * ssig2)
cbet2 = math.hypot(self._salp0, self._calp0 * csig2)
if cbet2 == 0:
# I.e., salp0 = 0, csig2 = 0. Break the degeneracy in this case
cbet2 = csig2 = Geodesic.tiny_
# tan(alp0) = cos(sig2)*tan(alp2)
salp2 = self._salp0; calp2 = self._calp0 * csig2 # No need to normalize
if outmask & Geodesic.DISTANCE:
# for s and t by r^3 and r, resp.
S = p * q / 4 # S = r^3 * s
r2 = Math.sq(r)
r3 = r * r2
# The discriminant of the quadratic equation for T3. This is zero on
# the evolute curve p^(1/3)+q^(1/3) = 1
disc = S * (S + 2 * r3)
u = r
if disc >= 0:
T3 = S + r3
# Pick the sign on the sqrt to maximize abs(T3). This minimizes loss
# of precision due to cancellation. The result is unchanged because
# of the way the T is used in definition of u.
T3 += -math.sqrt(disc) if T3 < 0 else math.sqrt(disc) # T3 = (r * t)^3
# N.B. cbrt always returns the real root. cbrt(-8) = -2.
T = Math.cbrt(T3) # T = r * t
# T can be zero; but then r2 / T -> 0.
u += T + (r2 / T if T != 0 else 0)
else:
# T is complex, but the way u is defined the result is real.
ang = math.atan2(math.sqrt(-disc), -(S + r3))
# There are three possible cube roots. We choose the root which
# avoids cancellation. Note that disc < 0 implies that r < 0.
u += 2 * r * math.cos(ang / 3)
v = math.sqrt(Math.sq(u) + q) # guaranteed positive
# Avoid loss of accuracy when u < 0.
uv = q / (v - u) if u < 0 else u + v # u+v, guaranteed positive
w = (uv - q) / (2 * v) # positive?
# Rearrange expression for k to avoid loss of accuracy due to
# subtraction. Division by 0 not possible because uv > 0, w >= 0.
k = uv / (math.sqrt(uv + Math.sq(w)) + w) # guaranteed positive
else: # q == 0 && r <= 0
# for s and t by r^3 and r, resp.
S = p * q / 4 # S = r^3 * s
r2 = Math.sq(r)
r3 = r * r2
# The discrimant of the quadratic equation for T3. This is zero on
# the evolute curve p^(1/3)+q^(1/3) = 1
disc = S * (S + 2 * r3)
u = r
if (disc >= 0):
T3 = S + r3
# Pick the sign on the sqrt to maximize abs(T3). This minimizes loss
# of precision due to cancellation. The result is unchanged because
# of the way the T is used in definition of u.
T3 += -math.sqrt(disc) if T3 < 0 else math.sqrt(disc) # T3 = (r * t)^3
# N.B. cbrt always returns the real root. cbrt(-8) = -2.
T = Math.cbrt(T3) # T = r * t
# T can be zero; but then r2 / T -> 0.
u += T + (r2 / T if T != 0 else 0)
else:
# T is complex, but the way u is defined the result is real.
ang = math.atan2(math.sqrt(-disc), -(S + r3))
# There are three possible cube roots. We choose the root which
# avoids cancellation. Note that disc < 0 implies that r < 0.
u += 2 * r * math.cos(ang / 3)
v = math.sqrt(Math.sq(u) + q) # guaranteed positive
# Avoid loss of accuracy when u < 0.
uv = q / (v - u) if u < 0 else u + v # u+v, guaranteed positive
w = (uv - q) / (2 * v) # positive?
# Rearrange expression for k to avoid loss of accuracy due to
# subtraction. Division by 0 not possible because uv > 0, w >= 0.
k = uv / (math.sqrt(uv + Math.sq(w)) + w) # guaranteed positive
else: # q == 0 && r <= 0
def InverseStart(self, sbet1, cbet1, dn1, sbet2, cbet2, dn2, lam12,
# Scratch areas of the right size
C1a, C2a):
"""Private: Find a starting value for Newton's method."""
# Return a starting point for Newton's method in salp1 and calp1 (function
# value is -1). If Newton's method doesn't need to be used, return also
# salp2 and calp2 and function value is sig12.
sig12 = -1; salp2 = calp2 = dnm = Math.nan # Return values
# bet12 = bet2 - bet1 in [0, pi); bet12a = bet2 + bet1 in (-pi, 0]
sbet12 = sbet2 * cbet1 - cbet2 * sbet1
cbet12 = cbet2 * cbet1 + sbet2 * sbet1
# Volatile declaration needed to fix inverse cases
# 88.202499451857 0 -88.202499451857 179.981022032992859592
# 89.262080389218 0 -89.262080389218 179.992207982775375662
# 89.333123580033 0 -89.333123580032997687 179.99295812360148422
# which otherwise fail with g++ 4.4.4 x86 -O3
sbet12a = sbet2 * cbet1
sbet12a += cbet2 * sbet1
shortline = cbet12 >= 0 and sbet12 < 0.5 and cbet2 * lam12 < 0.5
omg12 = lam12
if shortline:
sbetm2 = Math.sq(sbet1 + sbet2)
# sin((bet1+bet2)/2)^2
# for s and t by r^3 and r, resp.
S = p * q / 4 # S = r^3 * s
r2 = Math.sq(r)
r3 = r * r2
# The discriminant of the quadratic equation for T3. This is zero on
# the evolute curve p^(1/3)+q^(1/3) = 1
disc = S * (S + 2 * r3)
u = r
if disc >= 0:
T3 = S + r3
# Pick the sign on the sqrt to maximize abs(T3). This minimizes loss
# of precision due to cancellation. The result is unchanged because
# of the way the T is used in definition of u.
T3 += -math.sqrt(disc) if T3 < 0 else math.sqrt(disc) # T3 = (r * t)^3
# N.B. cbrt always returns the real root. cbrt(-8) = -2.
T = Math.cbrt(T3) # T = r * t
# T can be zero; but then r2 / T -> 0.
u += T + (r2 / T if T != 0 else 0)
else:
# T is complex, but the way u is defined the result is real.
ang = math.atan2(math.sqrt(-disc), -(S + r3))
# There are three possible cube roots. We choose the root which
# avoids cancellation. Note that disc < 0 implies that r < 0.
u += 2 * r * math.cos(ang / 3)
v = math.sqrt(Math.sq(u) + q) # guaranteed positive
# Avoid loss of accuracy when u < 0.
uv = q / (v - u) if u < 0 else u + v # u+v, guaranteed positive
w = (uv - q) / (2 * v) # positive?
# Rearrange expression for k to avoid loss of accuracy due to
# subtraction. Division by 0 not possible because uv > 0, w >= 0.
k = uv / (math.sqrt(uv + Math.sq(w)) + w) # guaranteed positive
else: # q == 0 && r <= 0