Secure your code as it's written. Use Snyk Code to scan source code in minutes - no build needed - and fix issues immediately.
s2, c2, s1, c1, _, c21 = sincos2(phi2, phi1, lam21)
E = datum.ellipsoid
if E.f and abs(c1) > EPS and abs(c2) > EPS:
r1 = atan(E.b_a * s1 / c1)
r2 = atan(E.b_a * s2 / c2)
j = (r2 + r1) / 2.0
k = (r2 - r1) / 2.0
sj, cj, sk, ck, sl_2, _ = sincos2(j, k, lam21 / 2.0)
h = fsum_(sk**2, (ck * sl_2)**2, -(sj * sl_2)**2)
if EPS < abs(h) < EPS1:
u = 1 / (1 - h)
d = 2 * atan(sqrt(h * u)) # == acos(1 - 2 * h)
sd, cd = sincos2(d)
if abs(sd) > EPS:
u = 2 * (sj * ck)**2 * u
v = 2 * (sk * cj)**2 / h
x = u + v
y = u - v
t = d / sd
s = 4 * t**2
e = 2 * cd
a = s * e
b = 2 * d
c = t - (a - e) / 2.0
s = fsum_(a * x, c * x**2, -b * y, -e * y**2, s * x * y) * E.f / 16.0
s = fsum_(t * x, -y, -s) * E.f / 4.0
return d - s * sd
limit and/or if this and the B{C{other}}
point are coincident or near-antipodal.
'''
E = self.ellipsoids(other)
c1, s1, _ = _r3(self.lat, E.f)
c2, s2, _ = _r3(other.lat, E.f)
c1c2, s1c2 = c1 * c2, s1 * c2
c1s2, s1s2 = c1 * s2, s1 * s2
dl, _ = unroll180(self.lon, other.lon, wrap=wrap)
ll = dl = radians(dl)
for self._iteration in range(self._iterations):
ll_ = ll
sll, cll = sincos2(ll)
ss = hypot(c2 * sll, c1s2 - s1c2 * cll)
if ss < EPS: # coincident or antipodal, ...
if self.isantipodeTo(other, eps=self._epsilon):
t = '%r %sto %r' % (self, _antipodal_, other)
raise VincentyError(_ambiguous_, txt=t)
# return zeros like Karney, but unlike Veness
return Distance3Tuple(0.0, 0, 0)
cs = s1s2 + c1c2 * cll
s = atan2(ss, cs)
sa = c1c2 * sll / ss
c2a = 1 - sa**2
if abs(c2a) < EPS:
c2a = 0 # equatorial line
i = radians(bearing) # initial bearing (forward azimuth)
si, ci = sincos2(i)
s12 = atan2(t1, ci) * 2
sa = c1 * si
c2a = 1 - sa**2
if c2a < EPS:
c2a = 0
A, B = 1, 0
else: # e22 == (a / b)**2 - 1
A, B = _p2(c2a * E.e22)
s = d = distance / (E.b * A)
for self._iteration in range(self._iterations):
ss, cs = sincos2(s)
c2sm = cos(s12 + s)
s_, s = s, d + _ds(B, cs, ss, c2sm)
if abs(s - s_) < self._epsilon:
break
else:
raise VincentyError(_no_convergence_, txt=repr(self)) # self.toRepr()
t = s1 * ss - c1 * cs * ci
# final bearing (reverse azimuth +/- 180)
r = degrees360(atan2(sa, -t))
if llr:
# destination latitude in [-90, 90)
a = degrees90(atan2(s1 * cs + c1 * ss * ci,
(1 - E.f) * hypot(sa, t)))
# destination longitude in [-180, 180)
def _reverse(self, x, y, name, LatLon, LatLon_kwds, _c_t, lea):
'''(INTERNAL) Azimuthal (spherical) reverse C{x, y} to C{lat, lon}.
'''
x = Scalar(x, name=_x_)
y = Scalar(y, name=_y_)
r = hypot(x, y)
c, t = _c_t(r / self.radius)
if t:
s0, c0 = self._sc0
sc, cc = sincos2(c)
k = c / sc
z = atan2d(x, y) # (x, y) for azimuth from true North
lat = degrees(asin1(s0 * cc + c0 * sc * (y / r)))
if lea or abs(c0) > EPS:
lon = atan2(x * sc, c0 * cc * r - s0 * sc * y)
else:
lon = atan2(x, (y if s0 < 0 else -y))
lon = _norm180(self.lon0 + degrees(lon))
else:
k, z = 1, 0
lat, lon = self.latlon0
t = self._toLatLon(lat, lon, LatLon, LatLon_kwds) if LatLon else \
Azimuthal7Tuple(x, y, lat, lon, z, k, self.datum)
return _xnamed(t, name or self.name)
@arg lam2: End longitude (C{radians}).
@kwarg final: Return final bearing if C{True}, initial
otherwise (C{bool}).
@kwarg wrap: Wrap and L{unrollPI} longitudes (C{bool}).
@return: Initial or final bearing (compass C{radiansPI2}) or
zero if start and end point coincide.
'''
if final:
phi1, lam1, phi2, lam2 = phi2, lam2, phi1, lam1
r = PI2 + PI
else:
r = PI2
db, _ = unrollPI(lam1, lam2, wrap=wrap)
sa1, ca1, sa2, ca2, sdb, cdb = sincos2(phi1, phi2, db)
# see
x = ca1 * sa2 - sa1 * ca2 * cdb
y = sdb * ca2
return (atan2(y, x) + r) % PI2 # wrapPI2
@return: Angular distance (C{radians}).
@see: Functions L{vincentys}, L{cosineAndoyerLambert_},
L{cosineForsytheAndoyerLambert_}, L{cosineLaw_},
L{equirectangular_}, L{euclidean_}, L{flatLocal_}/L{hubeny_},
L{flatPolar_}, L{haversine_} and L{thomas_}.
@note: Functions L{vincentys_}, L{haversine_} and L{cosineLaw_}
produce equivalent results, but L{vincentys_} is suitable
for antipodal points and slightly more expensive (M{3 cos,
3 sin, 1 hypot, 1 atan2, 6 mul, 2 add}) than L{haversine_}
(M{2 cos, 2 sin, 2 sqrt, 1 atan2, 5 mul, 1 add}) and
L{cosineLaw_} (M{3 cos, 3 sin, 1 acos, 3 mul, 1 add}).
'''
sa1, ca1, sa2, ca2, sb21, cb21 = sincos2(phi1, phi2, lam21)
c = ca2 * cb21
x = sa1 * sa2 + ca1 * c
y = ca1 * sa2 - sa1 * c
return atan2(hypot(ca2 * sb21, y), x)
@kwarg datum: Ellipsoidal datum to use (L{Datum}).
@return: Angular distance (C{radians}).
@raise TypeError: Invalid B{C{datum}}.
@see: Functions L{thomas}, L{cosineAndoyerLambert_},
L{cosineForsytheAndoyerLambert_}, L{cosineLaw_},
L{equirectangular_}, L{euclidean_}, L{flatLocal_}/L{hubeny_},
L{flatPolar_}, L{haversine_} and L{vincentys_} and U{Geodesy-PHP
}.
'''
_xinstanceof(Datum, datum=datum)
s2, c2, s1, c1, _, c21 = sincos2(phi2, phi1, lam21)
E = datum.ellipsoid
if E.f and abs(c1) > EPS and abs(c2) > EPS:
r1 = atan(E.b_a * s1 / c1)
r2 = atan(E.b_a * s2 / c2)
j = (r2 + r1) / 2.0
k = (r2 - r1) / 2.0
sj, cj, sk, ck, sl_2, _ = sincos2(j, k, lam21 / 2.0)
h = fsum_(sk**2, (ck * sl_2)**2, -(sj * sl_2)**2)
if EPS < abs(h) < EPS1:
u = 1 / (1 - h)
d = 2 * atan(sqrt(h * u)) # == acos(1 - 2 * h)
sd, cd = sincos2(d)
if abs(sd) > EPS:
u = 2 * (sj * ck)**2 * u
def _destination2(a, b, r, t):
'''(INTERNAL) Destination phi- and longitude in C{radians}.
@arg a: Latitude (C{radians}).
@arg b: Longitude (C{radians}).
@arg r: Angular distance (C{radians}).
@arg t: Bearing (compass C{radians}).
@return: 2-Tuple (phi, lam) of (C{radians}, C{radiansPI}).
'''
# see
sa, ca, sr, cr, st, ct = sincos2(a, r, t)
a = asin(ct * sr * ca + cr * sa)
d = atan2(st * sr * ca, cr - sa * sin(a))
# note, in EdWilliams.org/avform.htm W is + and E is -
return a, b + d
@arg lat: Latitude at the crossing (C{degrees}).
@kwarg wrap: Wrap and unroll longitudes (C{bool}).
@return: 2-Tuple C{(lon1, lon2)}, both in C{degrees180} or
C{None} if the great circle doesn't reach B{C{lat}}.
'''
self.others(other)
a1, b1 = self.philam
a2, b2 = other.philam
a = radians(lat)
db, b2 = unrollPI(b1, b2, wrap=wrap)
sa, ca, sa1, ca1, \
sa2, ca2, sdb, cdb = sincos2(a, a1, a2, db)
x = sa1 * ca2 * ca * sdb
y = sa1 * ca2 * ca * cdb - ca1 * sa2 * ca
z = ca1 * ca2 * sa * sdb
h = hypot(x, y)
if h < EPS or abs(z) > h:
return None # great circle doesn't reach latitude
m = atan2(-y, x) + b1 # longitude at max latitude
d = acos1(z / h) # delta longitude to intersections
return degrees180(m - d), degrees180(m + d)