Secure your code as it's written. Use Snyk Code to scan source code in minutes - no build needed - and fix issues immediately.
#
# longsign, swapp, latsign register the transformation to bring the
# coordinates to this canonical form. In all cases, 1 means no change was
# 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.
# real phi, sbet1, cbet1, sbet2, cbet2, s12x, m12x
sbet1, cbet1 = Math.sincosd(lat1); sbet1 *= self._f1
# Ensure cbet1 = +epsilon at poles
sbet1, cbet1 = Math.norm(sbet1, cbet1); cbet1 = max(Geodesic.tiny_, cbet1)
sbet2, cbet2 = Math.sincosd(lat2); sbet2 *= self._f1
# Ensure cbet2 = +epsilon at poles
sbet2, cbet2 = Math.norm(sbet2, cbet2); cbet2 = max(Geodesic.tiny_, cbet2)
# 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 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
if v > 0 and (numit > Geodesic.maxit1_ or
calp1/salp1 > calp1b/salp1b):
salp1b = salp1; calp1b = calp1
elif v < 0 and (numit > Geodesic.maxit1_ or
calp1/salp1 < calp1a/salp1a):
salp1a = salp1; calp1a = calp1
numit += 1
if numit < Geodesic.maxit1_ and dv > 0:
dalp1 = -v/dv
sdalp1 = math.sin(dalp1); cdalp1 = math.cos(dalp1)
nsalp1 = salp1 * cdalp1 + calp1 * sdalp1
if nsalp1 > 0 and abs(dalp1) < math.pi:
calp1 = calp1 * cdalp1 - salp1 * sdalp1
salp1 = nsalp1
salp1, calp1 = Math.norm(salp1, calp1)
# In some regimes we don't get quadratic convergence because
# slope -> 0. So use convergence conditions based on epsilon
# instead of sqrt(epsilon).
tripn = abs(v) <= 16 * Geodesic.tol0_
continue
# 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 = 90deg:
# the WGS84 test set: mean = 5.21, sd = 3.93, max = 24
# WGS84 and random input: mean = 4.74, sd = 0.99
salp1 = (salp1a + salp1b)/2
calp1 = (calp1a + calp1b)/2
salp1, calp1 = Math.norm(salp1, calp1)
"""Private: Solve hybrid problem"""
if sbet1 == 0 and calp1 == 0:
# Break degeneracy of equatorial line. This case has already been
# handled.
calp1 = -Geodesic.tiny_
# sin(alp1) * cos(bet1) = sin(alp0)
salp0 = salp1 * cbet1
calp0 = math.hypot(calp1, salp1 * sbet1) # calp0 > 0
# real somg1, comg1, somg2, comg2, lam12
# tan(bet1) = tan(sig1) * cos(alp1)
# tan(omg1) = sin(alp0) * tan(sig1) = tan(omg1)=tan(alp1)*sin(bet1)
ssig1 = sbet1; somg1 = salp0 * sbet1
csig1 = comg1 = calp1 * cbet1
ssig1, csig1 = Math.norm(ssig1, csig1)
# Math.norm(somg1, comg1); -- don't need to normalize!
# 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)
salp2 = salp0 / cbet2 if cbet2 != cbet1 else salp1
# 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]).
calp2 = (math.sqrt(Math.sq(calp1 * cbet1) +
((cbet2 - cbet1) * (cbet1 + cbet2) if cbet1 < -sbet1
else (sbet1 - sbet2) * (sbet1 + sbet2))) / cbet2
if cbet2 != cbet1 or abs(sbet2) != -sbet1 else abs(calp1))
# tan(bet2) = tan(sig2) * cos(alp2)
#
# 0 <= lon12 <= 180
# -90 <= lat1 <= 0
# lat1 <= lat2 <= -lat1
#
# longsign, swapp, latsign register the transformation to bring the
# coordinates to this canonical form. In all cases, 1 means no change was
# 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.
# real phi, sbet1, cbet1, sbet2, cbet2, s12x, m12x
sbet1, cbet1 = Math.sincosd(lat1); sbet1 *= self._f1
# Ensure cbet1 = +epsilon at poles
sbet1, cbet1 = Math.norm(sbet1, cbet1); cbet1 = max(Geodesic.tiny_, cbet1)
sbet2, cbet2 = Math.sincosd(lat2); sbet2 *= self._f1
# Ensure cbet2 = +epsilon at poles
sbet2, cbet2 = Math.norm(sbet2, cbet2); cbet2 = max(Geodesic.tiny_, cbet2)
# 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 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
"""the longitude of the first point in degrees (readonly)"""
if Math.isnan(salp1) or Math.isnan(calp1):
self.azi1 = Math.AngNormalize(azi1)
self.salp1, self.calp1 = Math.sincosd(Math.AngRound(azi1))
else:
self.azi1 = azi1
"""the azimuth at the first point in degrees (readonly)"""
self.salp1 = salp1
"""the sine of the azimuth at the first point (readonly)"""
self.calp1 = calp1
"""the cosine of the azimuth at the first point (readonly)"""
# real cbet1, sbet1
sbet1, cbet1 = Math.sincosd(Math.AngRound(lat1)); sbet1 *= self._f1
# Ensure cbet1 = +epsilon at poles
sbet1, cbet1 = Math.norm(sbet1, cbet1); cbet1 = max(Geodesic.tiny_, cbet1)
self._dn1 = math.sqrt(1 + geod._ep2 * Math.sq(sbet1))
# Evaluate alp0 from sin(alp1) * cos(bet1) = sin(alp0),
self._salp0 = self.salp1 * cbet1 # alp0 in [0, pi/2 - |bet1|]
# Alt: calp0 = hypot(sbet1, calp1 * cbet1). The following
# is slightly better (consider the case salp1 = 0).
self._calp0 = math.hypot(self.calp1, self.salp1 * sbet1)
# Evaluate sig with tan(bet1) = tan(sig1) * cos(alp1).
# sig = 0 is nearest northward crossing of equator.
# With bet1 = 0, alp1 = pi/2, we have sig1 = 0 (equatorial line).
# With bet1 = pi/2, alp1 = -pi, sig1 = pi/2
# With bet1 = -pi/2, alp1 = 0 , sig1 = -pi/2
# Evaluate omg1 with tan(omg1) = sin(alp0) * tan(sig1).
# With alp0 in (0, pi/2], quadrants for sig and omg coincide.
# No atan2(0,0) ambiguity at poles since cbet1 = +epsilon.
# With alp0 = 0, omg1 = 0 for alp1 = 0, omg1 = pi for alp1 = pi.
# 3 36189 32341
# 4 5396 7
# 5 455 1
# 6 56 0
#
# Because omg12 is near pi, estimate work with omg12a = pi - omg12
k = Geodesic._Astroid(x, y)
omg12a = lamscale * ( -x * k/(1 + k) if self.f >= 0
else -y * (1 + k)/k )
somg12 = math.sin(omg12a); comg12 = -math.cos(omg12a)
# Update spherical estimate of alp1 using omg12 instead of lam12
salp1 = cbet2 * somg12
calp1 = sbet12a - cbet2 * sbet1 * Math.sq(somg12) / (1 - comg12)
# Sanity check on starting guess. Backwards check allows NaN through.
if not (salp1 <= 0):
salp1, calp1 = Math.norm(salp1, calp1)
else:
salp1 = 1; calp1 = 0
return sig12, salp1, calp1, salp2, calp2, dnm
m12 = 0.0 + m12x # Convert -0 to 0
if outmask & Geodesic.AREA:
# From Lambda12: sin(alp1) * cos(bet1) = sin(alp0)
salp0 = salp1 * cbet1
calp0 = math.hypot(calp1, salp1 * sbet1) # calp0 > 0
# real alp12
if calp0 != 0 and salp0 != 0:
# From Lambda12: tan(bet) = tan(sig) * cos(alp)
ssig1 = sbet1; csig1 = calp1 * cbet1
ssig2 = sbet2; csig2 = calp2 * cbet2
k2 = Math.sq(calp0) * self._ep2
eps = k2 / (2 * (1 + math.sqrt(1 + k2)) + k2)
# Multiplier = a^2 * e^2 * cos(alpha0) * sin(alpha0).
A4 = Math.sq(self.a) * calp0 * salp0 * self._e2
ssig1, csig1 = Math.norm(ssig1, csig1)
ssig2, csig2 = Math.norm(ssig2, csig2)
C4a = list(range(Geodesic.nC4_))
self._C4f(eps, C4a)
B41 = Geodesic._SinCosSeries(False, ssig1, csig1, C4a)
B42 = Geodesic._SinCosSeries(False, ssig2, csig2, C4a)
S12 = A4 * (B42 - B41)
else:
# Avoid problems with indeterminate sig1, sig2 on equator
S12 = 0.0
if not meridian and somg12 > 1:
somg12 = math.sin(omg12); comg12 = math.cos(omg12)
if (not meridian and
# omg12 < 3/4 * pi
comg12 > -0.7071 and # Long difference not too big
# 3 36189 32341
# 4 5396 7
# 5 455 1
# 6 56 0
#
# Because omg12 is near pi, estimate work with omg12a = pi - omg12
k = Geodesic._Astroid(x, y)
omg12a = lamscale * ( -x * k/(1 + k) if self.f >= 0
else -y * (1 + k)/k )
somg12 = math.sin(omg12a); comg12 = -math.cos(omg12a)
# Update spherical estimate of alp1 using omg12 instead of lam12
salp1 = cbet2 * somg12
calp1 = sbet12a - cbet2 * sbet1 * Math.sq(somg12) / (1 - comg12)
# Sanity check on starting guess. Backwards check allows NaN through.
if not (salp1 <= 0):
salp1, calp1 = Math.norm(salp1, calp1)
else:
salp1 = 1; calp1 = 0
return sig12, salp1, calp1, salp2, calp2, dnm
# iteration.
# sin(alp2) * cos(bet2) = sin(alp0)
salp2 = salp0 / cbet2 if cbet2 != cbet1 else salp1
# 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]).
calp2 = (math.sqrt(Math.sq(calp1 * cbet1) +
((cbet2 - cbet1) * (cbet1 + cbet2) if cbet1 < -sbet1
else (sbet1 - sbet2) * (sbet1 + sbet2))) / cbet2
if cbet2 != cbet1 or abs(sbet2) != -sbet1 else abs(calp1))
# tan(bet2) = tan(sig2) * cos(alp2)
# tan(omg2) = sin(alp0) * tan(sig2).
ssig2 = sbet2; somg2 = salp0 * sbet2
csig2 = comg2 = calp2 * cbet2
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)