Secure your code as it's written. Use Snyk Code to scan source code in minutes - no build needed - and fix issues immediately.
dummy, m12b, m0, dummy, dummy = self.Lengths(
self._n, math.pi + bet12a, sbet1, -cbet1, dn1, sbet2, cbet2, dn2,
cbet1, cbet2, False, C1a, C2a)
x = -1 + m12b / (cbet1 * cbet2 * m0 * math.pi)
betscale = (sbet12a / x if x < -0.01
else -self._f * Math.sq(cbet1) * math.pi)
lamscale = betscale / cbet1
y = (lam12 - math.pi) / lamscale
if y > -Geodesic.tol1_ and x > -1 - Geodesic.xthresh_:
# strip near cut
if self._f >= 0:
salp1 = min(1.0, -x); calp1 = - math.sqrt(1 - Math.sq(salp1))
else:
calp1 = max((0.0 if x > -Geodesic.tol1_ else -1.0), x)
salp1 = math.sqrt(1 - Math.sq(calp1))
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
ssig2, csig2 = Math.norm(ssig2, csig2)
# Math.norm(somg2, comg2); -- don't need to normalize!
# sig12 = sig2 - sig1, limit to [0, pi]
sig12 = math.atan2(max(0.0, csig1 * ssig2 - ssig1 * csig2),
csig1 * csig2 + ssig1 * ssig2)
# omg12 = omg2 - omg1, limit to [0, pi]
somg12 = max(0.0, comg1 * somg2 - somg1 * comg2)
comg12 = comg1 * comg2 + somg1 * somg2
# eta = omg12 - lam120
eta = math.atan2(somg12 * clam120 - comg12 * slam120,
comg12 * clam120 + somg12 * slam120)
# real B312
k2 = Math.sq(calp0) * self._ep2
eps = k2 / (2 * (1 + math.sqrt(1 + k2)) + k2)
self._C3f(eps, C3a)
B312 = (Geodesic._SinCosSeries(True, ssig2, csig2, C3a) -
Geodesic._SinCosSeries(True, ssig1, csig1, C3a))
domg12 = -self.f * self._A3f(eps) * salp0 * (sig12 + B312)
lam12 = eta + domg12
if diffp:
if calp2 == 0:
dlam12 = - 2 * self._f1 * dn1 / sbet1
else:
dummy, dlam12, dummy, dummy, dummy = self._Lengths(
eps, sig12, ssig1, csig1, dn1, ssig2, csig2, dn2, cbet1, cbet2,
Geodesic.REDUCEDLENGTH, C1a, C2a)
dlam12 *= self._f1 / (calp2 * cbet2)
else:
def _Astroid(x, y):
"""Private: solve astroid equation."""
# Solve k^4+2*k^3-(x^2+y^2-1)*k^2-2*y^2*k-y^2 = 0 for positive root k.
# This solution is adapted from Geocentric::Reverse.
p = Math.sq(x)
q = Math.sq(y)
r = (p + q - 1) / 6
if not(q == 0 and r <= 0):
# Avoid possible division by zero when r = 0 by multiplying equations
# 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:
elif (abs(self._n) >= 0.1 or # Skip astroid calc if too eccentric
csig12 >= 0 or
ssig12 >= 6 * abs(self._n) * math.pi * Math.sq(cbet1)):
# Nothing to do, zeroth order spherical approximation is OK
pass
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.
# real y, lamscale, betscale
# Volatile declaration needed to fix inverse case
# 56.320923501171 0 -56.320923501171 179.664747671772880215
# which otherwise fails with g++ 4.4.4 x86 -O3
# volatile real x
if self._f >= 0: # In fact f == 0 does not get here
# x = dlong, y = dlat
k2 = Math.sq(sbet1) * self._ep2
eps = k2 / (2 * (1 + math.sqrt(1 + k2)) + k2)
lamscale = self._f * cbet1 * self.A3f(eps) * math.pi
betscale = lamscale * cbet1
x = (lam12 - math.pi) / lamscale
y = sbet12a / betscale
else: # _f < 0
# x = dlat, y = dlong
cbet12a = cbet2 * cbet1 - sbet2 * sbet1
bet12a = math.atan2(sbet12a, cbet12a)
# real m12b, m0, dummy
# In the case of lon12 = 180, this repeats a calculation made in
# Inverse.
dummy, m12b, m0, dummy, dummy = self.Lengths(
self._n, math.pi + bet12a, sbet1, -cbet1, dn1, sbet2, cbet2, dn2,
cbet1, cbet2, False, C1a, C2a)
x = -1 + m12b / (cbet1 * cbet2 * m0 * math.pi)
salp12 = Geodesic.tiny_ * self._calp1
calp12 = -1
else:
# tan(alp) = tan(alp0) * sec(sig)
# tan(alp2-alp1) = (tan(alp2) -tan(alp1)) / (tan(alp2)*tan(alp1)+1)
# = calp0 * salp0 * (csig1-csig2) / (salp0^2 + calp0^2 * csig1*csig2)
# If csig12 > 0, write
# csig1 - csig2 = ssig12 * (csig1 * ssig12 / (1 + csig12) + ssig1)
# else
# csig1 - csig2 = csig1 * (1 - csig12) + ssig12 * ssig1
# No need to normalize
salp12 = self._calp0 * self._salp0 * (
self._csig1 * (1 - csig12) + ssig12 * self._ssig1 if csig12 <= 0
else ssig12 * (self._csig1 * ssig12 / (1 + csig12) + self._ssig1))
calp12 = (Math.sq(self._salp0) +
Math.sq(self._calp0) * self._csig1 * csig2)
S12 = self._c2 * math.atan2(salp12, calp12) + self._A4 * (B42 - self._B41)
a12 = s12_a12 if arcmode else sig12 / Math.degree
return a12, lat2, lon2, azi2, s12, m12, M12, M21, S12
def _A1m1f(eps):
"""Private: return A1-1."""
coeff = [
1, 4, 64, 0, 256,
]
m = Geodesic.nA1_//2
t = Math.polyval(m, coeff, 0, Math.sq(eps)) / coeff[m + 1]
return (t + eps) / (1 - eps)
_A1m1f = staticmethod(_A1m1f)
def _A1m1f(eps):
"""Private: return A1-1."""
coeff = [
1, 4, 64, 0, 256,
]
m = Geodesic.nA1_//2
t = Math.polyval(m, coeff, 0, Math.sq(eps)) / coeff[m + 1]
return (t + eps) / (1 - eps)
_A1m1f = staticmethod(_A1m1f)
somg12 = math.sin(omg12); comg12 = math.cos(omg12)
else:
somg12 = slam12; comg12 = clam12
salp1 = cbet2 * somg12
calp1 = (
sbet12 + cbet2 * sbet1 * Math.sq(somg12) / (1 + comg12) if comg12 >= 0
else sbet12a - cbet2 * sbet1 * Math.sq(somg12) / (1 - comg12))
ssig12 = math.hypot(salp1, calp1)
csig12 = sbet1 * sbet2 + cbet1 * cbet2 * comg12
if shortline and ssig12 < self._etol2:
# really short lines
salp2 = cbet1 * somg12
calp2 = sbet12 - cbet1 * sbet2 * (Math.sq(somg12) / (1 + comg12)
if comg12 >= 0 else 1 - comg12)
salp2, calp2 = Math.norm(salp2, calp2)
# Set return value
sig12 = math.atan2(ssig12, csig12)
elif (abs(self._n) >= 0.1 or # Skip astroid calc if too eccentric
csig12 >= 0 or
ssig12 >= 6 * abs(self._n) * math.pi * Math.sq(cbet1)):
# Nothing to do, zeroth order spherical approximation is OK
pass
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.
# real y, lamscale, betscale
# Volatile declaration needed to fix inverse case
# 56.320923501171 0 -56.320923501171 179.664747671772880215
# which otherwise fails with g++ 4.4.4 x86 -O3
# |bet1| - |bet2|. Alternatively (cbet1 >= -sbet1), abs(sbet2) + sbet1 is
# a better measure. This logic is used in assigning calp2 in Lambda12.
# 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)
if cbet1 < -sbet1:
if cbet2 == cbet1:
sbet2 = sbet1 if sbet2 < 0 else -sbet1
else:
if abs(sbet2) == -sbet1:
cbet2 = cbet1
dn1 = math.sqrt(1 + self._ep2 * Math.sq(sbet1))
dn2 = math.sqrt(1 + self._ep2 * Math.sq(sbet2))
# real a12, sig12, calp1, salp1, calp2, salp2
# index zero elements of these arrays are unused
C1a = list(range(Geodesic.nC1_ + 1))
C2a = list(range(Geodesic.nC2_ + 1))
C3a = list(range(Geodesic.nC3_))
meridian = lat1 == -90 or slam12 == 0
if meridian:
# Endpoints are on a single full meridian, so the geodesic might lie on
# a meridian.
calp1 = clam12; salp1 = slam12 # Head to the target longitude
calp2 = 1.0; salp2 = 0.0 # At the target we're heading north
def A1m1f(eps):
"""Private: return A1-1."""
eps2 = Math.sq(eps)
t = eps2*(eps2*(eps2+4)+64)/256
return (t + eps) / (1 - eps)
A1m1f = staticmethod(A1m1f)