geodesic.py
Go to the documentation of this file.
1 """Define the :class:`~geographiclib.geodesic.Geodesic` class
2 
3 The ellipsoid parameters are defined by the constructor. The direct and
4 inverse geodesic problems are solved by
5 
6  * :meth:`~geographiclib.geodesic.Geodesic.Inverse` Solve the inverse
7  geodesic problem
8  * :meth:`~geographiclib.geodesic.Geodesic.Direct` Solve the direct
9  geodesic problem
10  * :meth:`~geographiclib.geodesic.Geodesic.ArcDirect` Solve the direct
11  geodesic problem in terms of spherical arc length
12 
13 :class:`~geographiclib.geodesicline.GeodesicLine` objects can be created
14 with
15 
16  * :meth:`~geographiclib.geodesic.Geodesic.Line`
17  * :meth:`~geographiclib.geodesic.Geodesic.DirectLine`
18  * :meth:`~geographiclib.geodesic.Geodesic.ArcDirectLine`
19  * :meth:`~geographiclib.geodesic.Geodesic.InverseLine`
20 
21 :class:`~geographiclib.polygonarea.PolygonArea` objects can be created
22 with
23 
24  * :meth:`~geographiclib.geodesic.Geodesic.Polygon`
25 
26 The public attributes for this class are
27 
28  * :attr:`~geographiclib.geodesic.Geodesic.a`
29  :attr:`~geographiclib.geodesic.Geodesic.f`
30 
31 *outmask* and *caps* bit masks are
32 
33  * :const:`~geographiclib.geodesic.Geodesic.EMPTY`
34  * :const:`~geographiclib.geodesic.Geodesic.LATITUDE`
35  * :const:`~geographiclib.geodesic.Geodesic.LONGITUDE`
36  * :const:`~geographiclib.geodesic.Geodesic.AZIMUTH`
37  * :const:`~geographiclib.geodesic.Geodesic.DISTANCE`
38  * :const:`~geographiclib.geodesic.Geodesic.STANDARD`
39  * :const:`~geographiclib.geodesic.Geodesic.DISTANCE_IN`
40  * :const:`~geographiclib.geodesic.Geodesic.REDUCEDLENGTH`
41  * :const:`~geographiclib.geodesic.Geodesic.GEODESICSCALE`
42  * :const:`~geographiclib.geodesic.Geodesic.AREA`
43  * :const:`~geographiclib.geodesic.Geodesic.ALL`
44  * :const:`~geographiclib.geodesic.Geodesic.LONG_UNROLL`
45 
46 :Example:
47 
48  >>> from geographiclib.geodesic import Geodesic
49  >>> # The geodesic inverse problem
50  ... Geodesic.WGS84.Inverse(-41.32, 174.81, 40.96, -5.50)
51  {'lat1': -41.32,
52  'a12': 179.6197069334283,
53  's12': 19959679.26735382,
54  'lat2': 40.96,
55  'azi2': 18.825195123248392,
56  'azi1': 161.06766998615882,
57  'lon1': 174.81,
58  'lon2': -5.5}
59 
60 """
61 # geodesic.py
62 #
63 # This is a rather literal translation of the GeographicLib::Geodesic class to
64 # python. See the documentation for the C++ class for more information at
65 #
66 # https://geographiclib.sourceforge.io/html/annotated.html
67 #
68 # The algorithms are derived in
69 #
70 # Charles F. F. Karney,
71 # Algorithms for geodesics, J. Geodesy 87, 43-55 (2013),
72 # https://doi.org/10.1007/s00190-012-0578-z
73 # Addenda: https://geographiclib.sourceforge.io/geod-addenda.html
74 #
75 # Copyright (c) Charles Karney (2011-2017) <charles@karney.com> and licensed
76 # under the MIT/X11 License. For more information, see
77 # https://geographiclib.sourceforge.io/
78 
79 
80 import math
81 from geographiclib.geomath import Math
82 from geographiclib.constants import Constants
83 from geographiclib.geodesiccapability import GeodesicCapability
84 
86  """Solve geodesic problems"""
87 
88  GEOGRAPHICLIB_GEODESIC_ORDER = 6
89  nA1_ = GEOGRAPHICLIB_GEODESIC_ORDER
90  nC1_ = GEOGRAPHICLIB_GEODESIC_ORDER
91  nC1p_ = GEOGRAPHICLIB_GEODESIC_ORDER
92  nA2_ = GEOGRAPHICLIB_GEODESIC_ORDER
93  nC2_ = GEOGRAPHICLIB_GEODESIC_ORDER
94  nA3_ = GEOGRAPHICLIB_GEODESIC_ORDER
95  nA3x_ = nA3_
96  nC3_ = GEOGRAPHICLIB_GEODESIC_ORDER
97  nC3x_ = (nC3_ * (nC3_ - 1)) // 2
98  nC4_ = GEOGRAPHICLIB_GEODESIC_ORDER
99  nC4x_ = (nC4_ * (nC4_ + 1)) // 2
100  maxit1_ = 20
101  maxit2_ = maxit1_ + Math.digits + 10
102 
103  tiny_ = math.sqrt(Math.minval)
104  tol0_ = Math.epsilon
105  tol1_ = 200 * tol0_
106  tol2_ = math.sqrt(tol0_)
107  tolb_ = tol0_ * tol2_
108  xthresh_ = 1000 * tol2_
109 
110  CAP_NONE = GeodesicCapability.CAP_NONE
111  CAP_C1 = GeodesicCapability.CAP_C1
112  CAP_C1p = GeodesicCapability.CAP_C1p
113  CAP_C2 = GeodesicCapability.CAP_C2
114  CAP_C3 = GeodesicCapability.CAP_C3
115  CAP_C4 = GeodesicCapability.CAP_C4
116  CAP_ALL = GeodesicCapability.CAP_ALL
117  CAP_MASK = GeodesicCapability.CAP_MASK
118  OUT_ALL = GeodesicCapability.OUT_ALL
119  OUT_MASK = GeodesicCapability.OUT_MASK
120 
121  def _SinCosSeries(sinp, sinx, cosx, c):
122  """Private: Evaluate a trig series using Clenshaw summation."""
123  # Evaluate
124  # y = sinp ? sum(c[i] * sin( 2*i * x), i, 1, n) :
125  # sum(c[i] * cos((2*i+1) * x), i, 0, n-1)
126  # using Clenshaw summation. N.B. c[0] is unused for sin series
127  # Approx operation count = (n + 5) mult and (2 * n + 2) add
128  k = len(c) # Point to one beyond last element
129  n = k - sinp
130  ar = 2 * (cosx - sinx) * (cosx + sinx) # 2 * cos(2 * x)
131  y1 = 0 # accumulators for sum
132  if n & 1:
133  k -= 1; y0 = c[k]
134  else:
135  y0 = 0
136  # Now n is even
137  n = n // 2
138  while n: # while n--:
139  n -= 1
140  # Unroll loop x 2, so accumulators return to their original role
141  k -= 1; y1 = ar * y0 - y1 + c[k]
142  k -= 1; y0 = ar * y1 - y0 + c[k]
143  return ( 2 * sinx * cosx * y0 if sinp # sin(2 * x) * y0
144  else cosx * (y0 - y1) ) # cos(x) * (y0 - y1)
145  _SinCosSeries = staticmethod(_SinCosSeries)
146 
147  def _Astroid(x, y):
148  """Private: solve astroid equation."""
149  # Solve k^4+2*k^3-(x^2+y^2-1)*k^2-2*y^2*k-y^2 = 0 for positive root k.
150  # This solution is adapted from Geocentric::Reverse.
151  p = Math.sq(x)
152  q = Math.sq(y)
153  r = (p + q - 1) / 6
154  if not(q == 0 and r <= 0):
155  # Avoid possible division by zero when r = 0 by multiplying equations
156  # for s and t by r^3 and r, resp.
157  S = p * q / 4 # S = r^3 * s
158  r2 = Math.sq(r)
159  r3 = r * r2
160  # The discriminant of the quadratic equation for T3. This is zero on
161  # the evolute curve p^(1/3)+q^(1/3) = 1
162  disc = S * (S + 2 * r3)
163  u = r
164  if disc >= 0:
165  T3 = S + r3
166  # Pick the sign on the sqrt to maximize abs(T3). This minimizes loss
167  # of precision due to cancellation. The result is unchanged because
168  # of the way the T is used in definition of u.
169  T3 += -math.sqrt(disc) if T3 < 0 else math.sqrt(disc) # T3 = (r * t)^3
170  # N.B. cbrt always returns the real root. cbrt(-8) = -2.
171  T = Math.cbrt(T3) # T = r * t
172  # T can be zero; but then r2 / T -> 0.
173  u += T + (r2 / T if T != 0 else 0)
174  else:
175  # T is complex, but the way u is defined the result is real.
176  ang = math.atan2(math.sqrt(-disc), -(S + r3))
177  # There are three possible cube roots. We choose the root which
178  # avoids cancellation. Note that disc < 0 implies that r < 0.
179  u += 2 * r * math.cos(ang / 3)
180  v = math.sqrt(Math.sq(u) + q) # guaranteed positive
181  # Avoid loss of accuracy when u < 0.
182  uv = q / (v - u) if u < 0 else u + v # u+v, guaranteed positive
183  w = (uv - q) / (2 * v) # positive?
184  # Rearrange expression for k to avoid loss of accuracy due to
185  # subtraction. Division by 0 not possible because uv > 0, w >= 0.
186  k = uv / (math.sqrt(uv + Math.sq(w)) + w) # guaranteed positive
187  else: # q == 0 && r <= 0
188  # y = 0 with |x| <= 1. Handle this case directly.
189  # for y small, positive root is k = abs(y)/sqrt(1-x^2)
190  k = 0
191  return k
192  _Astroid = staticmethod(_Astroid)
193 
194  def _A1m1f(eps):
195  """Private: return A1-1."""
196  coeff = [
197  1, 4, 64, 0, 256,
198  ]
199  m = Geodesic.nA1_//2
200  t = Math.polyval(m, coeff, 0, Math.sq(eps)) / coeff[m + 1]
201  return (t + eps) / (1 - eps)
202  _A1m1f = staticmethod(_A1m1f)
203 
204  def _C1f(eps, c):
205  """Private: return C1."""
206  coeff = [
207  -1, 6, -16, 32,
208  -9, 64, -128, 2048,
209  9, -16, 768,
210  3, -5, 512,
211  -7, 1280,
212  -7, 2048,
213  ]
214  eps2 = Math.sq(eps)
215  d = eps
216  o = 0
217  for l in range(1, Geodesic.nC1_ + 1): # l is index of C1p[l]
218  m = (Geodesic.nC1_ - l) // 2 # order of polynomial in eps^2
219  c[l] = d * Math.polyval(m, coeff, o, eps2) / coeff[o + m + 1]
220  o += m + 2
221  d *= eps
222  _C1f = staticmethod(_C1f)
223 
224  def _C1pf(eps, c):
225  """Private: return C1'"""
226  coeff = [
227  205, -432, 768, 1536,
228  4005, -4736, 3840, 12288,
229  -225, 116, 384,
230  -7173, 2695, 7680,
231  3467, 7680,
232  38081, 61440,
233  ]
234  eps2 = Math.sq(eps)
235  d = eps
236  o = 0
237  for l in range(1, Geodesic.nC1p_ + 1): # l is index of C1p[l]
238  m = (Geodesic.nC1p_ - l) // 2 # order of polynomial in eps^2
239  c[l] = d * Math.polyval(m, coeff, o, eps2) / coeff[o + m + 1]
240  o += m + 2
241  d *= eps
242  _C1pf = staticmethod(_C1pf)
243 
244  def _A2m1f(eps):
245  """Private: return A2-1"""
246  coeff = [
247  -11, -28, -192, 0, 256,
248  ]
249  m = Geodesic.nA2_//2
250  t = Math.polyval(m, coeff, 0, Math.sq(eps)) / coeff[m + 1]
251  return (t - eps) / (1 + eps)
252  _A2m1f = staticmethod(_A2m1f)
253 
254  def _C2f(eps, c):
255  """Private: return C2"""
256  coeff = [
257  1, 2, 16, 32,
258  35, 64, 384, 2048,
259  15, 80, 768,
260  7, 35, 512,
261  63, 1280,
262  77, 2048,
263  ]
264  eps2 = Math.sq(eps)
265  d = eps
266  o = 0
267  for l in range(1, Geodesic.nC2_ + 1): # l is index of C2[l]
268  m = (Geodesic.nC2_ - l) // 2 # order of polynomial in eps^2
269  c[l] = d * Math.polyval(m, coeff, o, eps2) / coeff[o + m + 1]
270  o += m + 2
271  d *= eps
272  _C2f = staticmethod(_C2f)
273 
274  def __init__(self, a, f):
275  """Construct a Geodesic object
276 
277  :param a: the equatorial radius of the ellipsoid in meters
278  :param f: the flattening of the ellipsoid
279 
280  An exception is thrown if *a* or the polar semi-axis *b* = *a* (1 -
281  *f*) is not a finite positive quantity.
282 
283  """
284 
285  self.a = float(a)
286  """The equatorial radius in meters (readonly)"""
287  self.f = float(f)
288  """The flattening (readonly)"""
289  self._f1 = 1 - self.f
290  self._e2 = self.f * (2 - self.f)
291  self._ep2 = self._e2 / Math.sq(self._f1) # e2 / (1 - e2)
292  self._n = self.f / ( 2 - self.f)
293  self._b = self.a * self._f1
294  # authalic radius squared
295  self._c2 = (Math.sq(self.a) + Math.sq(self._b) *
296  (1 if self._e2 == 0 else
297  (Math.atanh(math.sqrt(self._e2)) if self._e2 > 0 else
298  math.atan(math.sqrt(-self._e2))) /
299  math.sqrt(abs(self._e2))))/2
300  # The sig12 threshold for "really short". Using the auxiliary sphere
301  # solution with dnm computed at (bet1 + bet2) / 2, the relative error in
302  # the azimuth consistency check is sig12^2 * abs(f) * min(1, 1-f/2) / 2.
303  # (Error measured for 1/100 < b/a < 100 and abs(f) >= 1/1000. For a given
304  # f and sig12, the max error occurs for lines near the pole. If the old
305  # rule for computing dnm = (dn1 + dn2)/2 is used, then the error increases
306  # by a factor of 2.) Setting this equal to epsilon gives sig12 = etol2.
307  # Here 0.1 is a safety factor (error decreased by 100) and max(0.001,
308  # abs(f)) stops etol2 getting too large in the nearly spherical case.
309  self._etol2 = 0.1 * Geodesic.tol2_ / math.sqrt( max(0.001, abs(self.f)) *
310  min(1.0, 1-self.f/2) / 2 )
311  if not(Math.isfinite(self.a) and self.a > 0):
312  raise ValueError("Equatorial radius is not positive")
313  if not(Math.isfinite(self._b) and self._b > 0):
314  raise ValueError("Polar semi-axis is not positive")
315  self._A3x = list(range(Geodesic.nA3x_))
316  self._C3x = list(range(Geodesic.nC3x_))
317  self._C4x = list(range(Geodesic.nC4x_))
318  self._A3coeff()
319  self._C3coeff()
320  self._C4coeff()
321 
322  def _A3coeff(self):
323  """Private: return coefficients for A3"""
324  coeff = [
325  -3, 128,
326  -2, -3, 64,
327  -1, -3, -1, 16,
328  3, -1, -2, 8,
329  1, -1, 2,
330  1, 1,
331  ]
332  o = 0; k = 0
333  for j in range(Geodesic.nA3_ - 1, -1, -1): # coeff of eps^j
334  m = min(Geodesic.nA3_ - j - 1, j) # order of polynomial in n
335  self._A3x[k] = Math.polyval(m, coeff, o, self._n) / coeff[o + m + 1]
336  k += 1
337  o += m + 2
338 
339  def _C3coeff(self):
340  """Private: return coefficients for C3"""
341  coeff = [
342  3, 128,
343  2, 5, 128,
344  -1, 3, 3, 64,
345  -1, 0, 1, 8,
346  -1, 1, 4,
347  5, 256,
348  1, 3, 128,
349  -3, -2, 3, 64,
350  1, -3, 2, 32,
351  7, 512,
352  -10, 9, 384,
353  5, -9, 5, 192,
354  7, 512,
355  -14, 7, 512,
356  21, 2560,
357  ]
358  o = 0; k = 0
359  for l in range(1, Geodesic.nC3_): # l is index of C3[l]
360  for j in range(Geodesic.nC3_ - 1, l - 1, -1): # coeff of eps^j
361  m = min(Geodesic.nC3_ - j - 1, j) # order of polynomial in n
362  self._C3x[k] = Math.polyval(m, coeff, o, self._n) / coeff[o + m + 1]
363  k += 1
364  o += m + 2
365 
366  def _C4coeff(self):
367  """Private: return coefficients for C4"""
368  coeff = [
369  97, 15015,
370  1088, 156, 45045,
371  -224, -4784, 1573, 45045,
372  -10656, 14144, -4576, -858, 45045,
373  64, 624, -4576, 6864, -3003, 15015,
374  100, 208, 572, 3432, -12012, 30030, 45045,
375  1, 9009,
376  -2944, 468, 135135,
377  5792, 1040, -1287, 135135,
378  5952, -11648, 9152, -2574, 135135,
379  -64, -624, 4576, -6864, 3003, 135135,
380  8, 10725,
381  1856, -936, 225225,
382  -8448, 4992, -1144, 225225,
383  -1440, 4160, -4576, 1716, 225225,
384  -136, 63063,
385  1024, -208, 105105,
386  3584, -3328, 1144, 315315,
387  -128, 135135,
388  -2560, 832, 405405,
389  128, 99099,
390  ]
391  o = 0; k = 0
392  for l in range(Geodesic.nC4_): # l is index of C4[l]
393  for j in range(Geodesic.nC4_ - 1, l - 1, -1): # coeff of eps^j
394  m = Geodesic.nC4_ - j - 1 # order of polynomial in n
395  self._C4x[k] = Math.polyval(m, coeff, o, self._n) / coeff[o + m + 1]
396  k += 1
397  o += m + 2
398 
399  def _A3f(self, eps):
400  """Private: return A3"""
401  # Evaluate A3
402  return Math.polyval(Geodesic.nA3_ - 1, self._A3x, 0, eps)
403 
404  def _C3f(self, eps, c):
405  """Private: return C3"""
406  # Evaluate C3
407  # Elements c[1] thru c[nC3_ - 1] are set
408  mult = 1
409  o = 0
410  for l in range(1, Geodesic.nC3_): # l is index of C3[l]
411  m = Geodesic.nC3_ - l - 1 # order of polynomial in eps
412  mult *= eps
413  c[l] = mult * Math.polyval(m, self._C3x, o, eps)
414  o += m + 1
415 
416  def _C4f(self, eps, c):
417  """Private: return C4"""
418  # Evaluate C4 coeffs by Horner's method
419  # Elements c[0] thru c[nC4_ - 1] are set
420  mult = 1
421  o = 0
422  for l in range(Geodesic.nC4_): # l is index of C4[l]
423  m = Geodesic.nC4_ - l - 1 # order of polynomial in eps
424  c[l] = mult * Math.polyval(m, self._C4x, o, eps)
425  o += m + 1
426  mult *= eps
427 
428  # return s12b, m12b, m0, M12, M21
429  def _Lengths(self, eps, sig12,
430  ssig1, csig1, dn1, ssig2, csig2, dn2, cbet1, cbet2, outmask,
431  # Scratch areas of the right size
432  C1a, C2a):
433  """Private: return a bunch of lengths"""
434  # Return s12b, m12b, m0, M12, M21, where
435  # m12b = (reduced length)/_b; s12b = distance/_b,
436  # m0 = coefficient of secular term in expression for reduced length.
437  outmask &= Geodesic.OUT_MASK
438  # outmask & DISTANCE: set s12b
439  # outmask & REDUCEDLENGTH: set m12b & m0
440  # outmask & GEODESICSCALE: set M12 & M21
441 
442  s12b = m12b = m0 = M12 = M21 = Math.nan
443  if outmask & (Geodesic.DISTANCE | Geodesic.REDUCEDLENGTH |
444  Geodesic.GEODESICSCALE):
445  A1 = Geodesic._A1m1f(eps)
446  Geodesic._C1f(eps, C1a)
447  if outmask & (Geodesic.REDUCEDLENGTH | Geodesic.GEODESICSCALE):
448  A2 = Geodesic._A2m1f(eps)
449  Geodesic._C2f(eps, C2a)
450  m0x = A1 - A2
451  A2 = 1 + A2
452  A1 = 1 + A1
453  if outmask & Geodesic.DISTANCE:
454  B1 = (Geodesic._SinCosSeries(True, ssig2, csig2, C1a) -
455  Geodesic._SinCosSeries(True, ssig1, csig1, C1a))
456  # Missing a factor of _b
457  s12b = A1 * (sig12 + B1)
458  if outmask & (Geodesic.REDUCEDLENGTH | Geodesic.GEODESICSCALE):
459  B2 = (Geodesic._SinCosSeries(True, ssig2, csig2, C2a) -
460  Geodesic._SinCosSeries(True, ssig1, csig1, C2a))
461  J12 = m0x * sig12 + (A1 * B1 - A2 * B2)
462  elif outmask & (Geodesic.REDUCEDLENGTH | Geodesic.GEODESICSCALE):
463  # Assume here that nC1_ >= nC2_
464  for l in range(1, Geodesic.nC2_):
465  C2a[l] = A1 * C1a[l] - A2 * C2a[l]
466  J12 = m0x * sig12 + (Geodesic._SinCosSeries(True, ssig2, csig2, C2a) -
467  Geodesic._SinCosSeries(True, ssig1, csig1, C2a))
468  if outmask & Geodesic.REDUCEDLENGTH:
469  m0 = m0x
470  # Missing a factor of _b.
471  # Add parens around (csig1 * ssig2) and (ssig1 * csig2) to ensure
472  # accurate cancellation in the case of coincident points.
473  m12b = (dn2 * (csig1 * ssig2) - dn1 * (ssig1 * csig2) -
474  csig1 * csig2 * J12)
475  if outmask & Geodesic.GEODESICSCALE:
476  csig12 = csig1 * csig2 + ssig1 * ssig2
477  t = self._ep2 * (cbet1 - cbet2) * (cbet1 + cbet2) / (dn1 + dn2)
478  M12 = csig12 + (t * ssig2 - csig2 * J12) * ssig1 / dn1
479  M21 = csig12 - (t * ssig1 - csig1 * J12) * ssig2 / dn2
480  return s12b, m12b, m0, M12, M21
481 
482  # return sig12, salp1, calp1, salp2, calp2, dnm
483  def _InverseStart(self, sbet1, cbet1, dn1, sbet2, cbet2, dn2,
484  lam12, slam12, clam12,
485  # Scratch areas of the right size
486  C1a, C2a):
487  """Private: Find a starting value for Newton's method."""
488  # Return a starting point for Newton's method in salp1 and calp1 (function
489  # value is -1). If Newton's method doesn't need to be used, return also
490  # salp2 and calp2 and function value is sig12.
491  sig12 = -1; salp2 = calp2 = dnm = Math.nan # Return values
492  # bet12 = bet2 - bet1 in [0, pi); bet12a = bet2 + bet1 in (-pi, 0]
493  sbet12 = sbet2 * cbet1 - cbet2 * sbet1
494  cbet12 = cbet2 * cbet1 + sbet2 * sbet1
495  # Volatile declaration needed to fix inverse cases
496  # 88.202499451857 0 -88.202499451857 179.981022032992859592
497  # 89.262080389218 0 -89.262080389218 179.992207982775375662
498  # 89.333123580033 0 -89.333123580032997687 179.99295812360148422
499  # which otherwise fail with g++ 4.4.4 x86 -O3
500  sbet12a = sbet2 * cbet1
501  sbet12a += cbet2 * sbet1
502 
503  shortline = cbet12 >= 0 and sbet12 < 0.5 and cbet2 * lam12 < 0.5
504  if shortline:
505  sbetm2 = Math.sq(sbet1 + sbet2)
506  # sin((bet1+bet2)/2)^2
507  # = (sbet1 + sbet2)^2 / ((sbet1 + sbet2)^2 + (cbet1 + cbet2)^2)
508  sbetm2 /= sbetm2 + Math.sq(cbet1 + cbet2)
509  dnm = math.sqrt(1 + self._ep2 * sbetm2)
510  omg12 = lam12 / (self._f1 * dnm)
511  somg12 = math.sin(omg12); comg12 = math.cos(omg12)
512  else:
513  somg12 = slam12; comg12 = clam12
514 
515  salp1 = cbet2 * somg12
516  calp1 = (
517  sbet12 + cbet2 * sbet1 * Math.sq(somg12) / (1 + comg12) if comg12 >= 0
518  else sbet12a - cbet2 * sbet1 * Math.sq(somg12) / (1 - comg12))
519 
520  ssig12 = math.hypot(salp1, calp1)
521  csig12 = sbet1 * sbet2 + cbet1 * cbet2 * comg12
522 
523  if shortline and ssig12 < self._etol2:
524  # really short lines
525  salp2 = cbet1 * somg12
526  calp2 = sbet12 - cbet1 * sbet2 * (Math.sq(somg12) / (1 + comg12)
527  if comg12 >= 0 else 1 - comg12)
528  salp2, calp2 = Math.norm(salp2, calp2)
529  # Set return value
530  sig12 = math.atan2(ssig12, csig12)
531  elif (abs(self._n) >= 0.1 or # Skip astroid calc if too eccentric
532  csig12 >= 0 or
533  ssig12 >= 6 * abs(self._n) * math.pi * Math.sq(cbet1)):
534  # Nothing to do, zeroth order spherical approximation is OK
535  pass
536  else:
537  # Scale lam12 and bet2 to x, y coordinate system where antipodal point
538  # is at origin and singular point is at y = 0, x = -1.
539  # real y, lamscale, betscale
540  # Volatile declaration needed to fix inverse case
541  # 56.320923501171 0 -56.320923501171 179.664747671772880215
542  # which otherwise fails with g++ 4.4.4 x86 -O3
543  # volatile real x
544  lam12x = math.atan2(-slam12, -clam12)
545  if self.f >= 0: # In fact f == 0 does not get here
546  # x = dlong, y = dlat
547  k2 = Math.sq(sbet1) * self._ep2
548  eps = k2 / (2 * (1 + math.sqrt(1 + k2)) + k2)
549  lamscale = self.f * cbet1 * self._A3f(eps) * math.pi
550  betscale = lamscale * cbet1
551  x = lam12x / lamscale
552  y = sbet12a / betscale
553  else: # _f < 0
554  # x = dlat, y = dlong
555  cbet12a = cbet2 * cbet1 - sbet2 * sbet1
556  bet12a = math.atan2(sbet12a, cbet12a)
557  # real m12b, m0, dummy
558  # In the case of lon12 = 180, this repeats a calculation made in
559  # Inverse.
560  dummy, m12b, m0, dummy, dummy = self._Lengths(
561  self._n, math.pi + bet12a, sbet1, -cbet1, dn1, sbet2, cbet2, dn2,
562  cbet1, cbet2, Geodesic.REDUCEDLENGTH, C1a, C2a)
563  x = -1 + m12b / (cbet1 * cbet2 * m0 * math.pi)
564  betscale = (sbet12a / x if x < -0.01
565  else -self.f * Math.sq(cbet1) * math.pi)
566  lamscale = betscale / cbet1
567  y = lam12x / lamscale
568 
569  if y > -Geodesic.tol1_ and x > -1 - Geodesic.xthresh_:
570  # strip near cut
571  if self.f >= 0:
572  salp1 = min(1.0, -x); calp1 = - math.sqrt(1 - Math.sq(salp1))
573  else:
574  calp1 = max((0.0 if x > -Geodesic.tol1_ else -1.0), x)
575  salp1 = math.sqrt(1 - Math.sq(calp1))
576  else:
577  # Estimate alp1, by solving the astroid problem.
578  #
579  # Could estimate alpha1 = theta + pi/2, directly, i.e.,
580  # calp1 = y/k; salp1 = -x/(1+k); for _f >= 0
581  # calp1 = x/(1+k); salp1 = -y/k; for _f < 0 (need to check)
582  #
583  # However, it's better to estimate omg12 from astroid and use
584  # spherical formula to compute alp1. This reduces the mean number of
585  # Newton iterations for astroid cases from 2.24 (min 0, max 6) to 2.12
586  # (min 0 max 5). The changes in the number of iterations are as
587  # follows:
588  #
589  # change percent
590  # 1 5
591  # 0 78
592  # -1 16
593  # -2 0.6
594  # -3 0.04
595  # -4 0.002
596  #
597  # The histogram of iterations is (m = number of iterations estimating
598  # alp1 directly, n = number of iterations estimating via omg12, total
599  # number of trials = 148605):
600  #
601  # iter m n
602  # 0 148 186
603  # 1 13046 13845
604  # 2 93315 102225
605  # 3 36189 32341
606  # 4 5396 7
607  # 5 455 1
608  # 6 56 0
609  #
610  # Because omg12 is near pi, estimate work with omg12a = pi - omg12
611  k = Geodesic._Astroid(x, y)
612  omg12a = lamscale * ( -x * k/(1 + k) if self.f >= 0
613  else -y * (1 + k)/k )
614  somg12 = math.sin(omg12a); comg12 = -math.cos(omg12a)
615  # Update spherical estimate of alp1 using omg12 instead of lam12
616  salp1 = cbet2 * somg12
617  calp1 = sbet12a - cbet2 * sbet1 * Math.sq(somg12) / (1 - comg12)
618  # Sanity check on starting guess. Backwards check allows NaN through.
619  if not (salp1 <= 0):
620  salp1, calp1 = Math.norm(salp1, calp1)
621  else:
622  salp1 = 1; calp1 = 0
623  return sig12, salp1, calp1, salp2, calp2, dnm
624 
625  # return lam12, salp2, calp2, sig12, ssig1, csig1, ssig2, csig2, eps,
626  # domg12, dlam12
627  def _Lambda12(self, sbet1, cbet1, dn1, sbet2, cbet2, dn2, salp1, calp1,
628  slam120, clam120, diffp,
629  # Scratch areas of the right size
630  C1a, C2a, C3a):
631  """Private: Solve hybrid problem"""
632  if sbet1 == 0 and calp1 == 0:
633  # Break degeneracy of equatorial line. This case has already been
634  # handled.
635  calp1 = -Geodesic.tiny_
636 
637  # sin(alp1) * cos(bet1) = sin(alp0)
638  salp0 = salp1 * cbet1
639  calp0 = math.hypot(calp1, salp1 * sbet1) # calp0 > 0
640 
641  # real somg1, comg1, somg2, comg2, lam12
642  # tan(bet1) = tan(sig1) * cos(alp1)
643  # tan(omg1) = sin(alp0) * tan(sig1) = tan(omg1)=tan(alp1)*sin(bet1)
644  ssig1 = sbet1; somg1 = salp0 * sbet1
645  csig1 = comg1 = calp1 * cbet1
646  ssig1, csig1 = Math.norm(ssig1, csig1)
647  # Math.norm(somg1, comg1); -- don't need to normalize!
648 
649  # Enforce symmetries in the case abs(bet2) = -bet1. Need to be careful
650  # about this case, since this can yield singularities in the Newton
651  # iteration.
652  # sin(alp2) * cos(bet2) = sin(alp0)
653  salp2 = salp0 / cbet2 if cbet2 != cbet1 else salp1
654  # calp2 = sqrt(1 - sq(salp2))
655  # = sqrt(sq(calp0) - sq(sbet2)) / cbet2
656  # and subst for calp0 and rearrange to give (choose positive sqrt
657  # to give alp2 in [0, pi/2]).
658  calp2 = (math.sqrt(Math.sq(calp1 * cbet1) +
659  ((cbet2 - cbet1) * (cbet1 + cbet2) if cbet1 < -sbet1
660  else (sbet1 - sbet2) * (sbet1 + sbet2))) / cbet2
661  if cbet2 != cbet1 or abs(sbet2) != -sbet1 else abs(calp1))
662  # tan(bet2) = tan(sig2) * cos(alp2)
663  # tan(omg2) = sin(alp0) * tan(sig2).
664  ssig2 = sbet2; somg2 = salp0 * sbet2
665  csig2 = comg2 = calp2 * cbet2
666  ssig2, csig2 = Math.norm(ssig2, csig2)
667  # Math.norm(somg2, comg2); -- don't need to normalize!
668 
669  # sig12 = sig2 - sig1, limit to [0, pi]
670  sig12 = math.atan2(max(0.0, csig1 * ssig2 - ssig1 * csig2),
671  csig1 * csig2 + ssig1 * ssig2)
672 
673  # omg12 = omg2 - omg1, limit to [0, pi]
674  somg12 = max(0.0, comg1 * somg2 - somg1 * comg2)
675  comg12 = comg1 * comg2 + somg1 * somg2
676  # eta = omg12 - lam120
677  eta = math.atan2(somg12 * clam120 - comg12 * slam120,
678  comg12 * clam120 + somg12 * slam120)
679 
680  # real B312
681  k2 = Math.sq(calp0) * self._ep2
682  eps = k2 / (2 * (1 + math.sqrt(1 + k2)) + k2)
683  self._C3f(eps, C3a)
684  B312 = (Geodesic._SinCosSeries(True, ssig2, csig2, C3a) -
685  Geodesic._SinCosSeries(True, ssig1, csig1, C3a))
686  domg12 = -self.f * self._A3f(eps) * salp0 * (sig12 + B312)
687  lam12 = eta + domg12
688 
689  if diffp:
690  if calp2 == 0:
691  dlam12 = - 2 * self._f1 * dn1 / sbet1
692  else:
693  dummy, dlam12, dummy, dummy, dummy = self._Lengths(
694  eps, sig12, ssig1, csig1, dn1, ssig2, csig2, dn2, cbet1, cbet2,
695  Geodesic.REDUCEDLENGTH, C1a, C2a)
696  dlam12 *= self._f1 / (calp2 * cbet2)
697  else:
698  dlam12 = Math.nan
699 
700  return (lam12, salp2, calp2, sig12, ssig1, csig1, ssig2, csig2, eps,
701  domg12, dlam12)
702 
703  # return a12, s12, salp1, calp1, salp2, calp2, m12, M12, M21, S12
704  def _GenInverse(self, lat1, lon1, lat2, lon2, outmask):
705  """Private: General version of the inverse problem"""
706  a12 = s12 = m12 = M12 = M21 = S12 = Math.nan # return vals
707 
708  outmask &= Geodesic.OUT_MASK
709  # Compute longitude difference (AngDiff does this carefully). Result is
710  # in [-180, 180] but -180 is only for west-going geodesics. 180 is for
711  # east-going and meridional geodesics.
712  lon12, lon12s = Math.AngDiff(lon1, lon2)
713  # Make longitude difference positive.
714  lonsign = 1 if lon12 >= 0 else -1
715  # If very close to being on the same half-meridian, then make it so.
716  lon12 = lonsign * Math.AngRound(lon12)
717  lon12s = Math.AngRound((180 - lon12) - lonsign * lon12s)
718  lam12 = math.radians(lon12)
719  if lon12 > 90:
720  slam12, clam12 = Math.sincosd(lon12s); clam12 = -clam12
721  else:
722  slam12, clam12 = Math.sincosd(lon12)
723 
724  # If really close to the equator, treat as on equator.
725  lat1 = Math.AngRound(Math.LatFix(lat1))
726  lat2 = Math.AngRound(Math.LatFix(lat2))
727  # Swap points so that point with higher (abs) latitude is point 1
728  # If one latitude is a nan, then it becomes lat1.
729  swapp = -1 if abs(lat1) < abs(lat2) else 1
730  if swapp < 0:
731  lonsign *= -1
732  lat2, lat1 = lat1, lat2
733  # Make lat1 <= 0
734  latsign = 1 if lat1 < 0 else -1
735  lat1 *= latsign
736  lat2 *= latsign
737  # Now we have
738  #
739  # 0 <= lon12 <= 180
740  # -90 <= lat1 <= 0
741  # lat1 <= lat2 <= -lat1
742  #
743  # longsign, swapp, latsign register the transformation to bring the
744  # coordinates to this canonical form. In all cases, 1 means no change was
745  # made. We make these transformations so that there are few cases to
746  # check, e.g., on verifying quadrants in atan2. In addition, this
747  # enforces some symmetries in the results returned.
748 
749  # real phi, sbet1, cbet1, sbet2, cbet2, s12x, m12x
750 
751  sbet1, cbet1 = Math.sincosd(lat1); sbet1 *= self._f1
752  # Ensure cbet1 = +epsilon at poles
753  sbet1, cbet1 = Math.norm(sbet1, cbet1); cbet1 = max(Geodesic.tiny_, cbet1)
754 
755  sbet2, cbet2 = Math.sincosd(lat2); sbet2 *= self._f1
756  # Ensure cbet2 = +epsilon at poles
757  sbet2, cbet2 = Math.norm(sbet2, cbet2); cbet2 = max(Geodesic.tiny_, cbet2)
758 
759  # If cbet1 < -sbet1, then cbet2 - cbet1 is a sensitive measure of the
760  # |bet1| - |bet2|. Alternatively (cbet1 >= -sbet1), abs(sbet2) + sbet1 is
761  # a better measure. This logic is used in assigning calp2 in Lambda12.
762  # Sometimes these quantities vanish and in that case we force bet2 = +/-
763  # bet1 exactly. An example where is is necessary is the inverse problem
764  # 48.522876735459 0 -48.52287673545898293 179.599720456223079643
765  # which failed with Visual Studio 10 (Release and Debug)
766 
767  if cbet1 < -sbet1:
768  if cbet2 == cbet1:
769  sbet2 = sbet1 if sbet2 < 0 else -sbet1
770  else:
771  if abs(sbet2) == -sbet1:
772  cbet2 = cbet1
773 
774  dn1 = math.sqrt(1 + self._ep2 * Math.sq(sbet1))
775  dn2 = math.sqrt(1 + self._ep2 * Math.sq(sbet2))
776 
777  # real a12, sig12, calp1, salp1, calp2, salp2
778  # index zero elements of these arrays are unused
779  C1a = list(range(Geodesic.nC1_ + 1))
780  C2a = list(range(Geodesic.nC2_ + 1))
781  C3a = list(range(Geodesic.nC3_))
782 
783  meridian = lat1 == -90 or slam12 == 0
784 
785  if meridian:
786 
787  # Endpoints are on a single full meridian, so the geodesic might lie on
788  # a meridian.
789 
790  calp1 = clam12; salp1 = slam12 # Head to the target longitude
791  calp2 = 1.0; salp2 = 0.0 # At the target we're heading north
792 
793  # tan(bet) = tan(sig) * cos(alp)
794  ssig1 = sbet1; csig1 = calp1 * cbet1
795  ssig2 = sbet2; csig2 = calp2 * cbet2
796 
797  # sig12 = sig2 - sig1
798  sig12 = math.atan2(max(0.0, csig1 * ssig2 - ssig1 * csig2),
799  csig1 * csig2 + ssig1 * ssig2)
800 
801  s12x, m12x, dummy, M12, M21 = self._Lengths(
802  self._n, sig12, ssig1, csig1, dn1, ssig2, csig2, dn2, cbet1, cbet2,
803  outmask | Geodesic.DISTANCE | Geodesic.REDUCEDLENGTH, C1a, C2a)
804 
805  # Add the check for sig12 since zero length geodesics might yield m12 <
806  # 0. Test case was
807  #
808  # echo 20.001 0 20.001 0 | GeodSolve -i
809  #
810  # In fact, we will have sig12 > pi/2 for meridional geodesic which is
811  # not a shortest path.
812  if sig12 < 1 or m12x >= 0:
813  if sig12 < 3 * Geodesic.tiny_:
814  sig12 = m12x = s12x = 0.0
815  m12x *= self._b
816  s12x *= self._b
817  a12 = math.degrees(sig12)
818  else:
819  # m12 < 0, i.e., prolate and too close to anti-podal
820  meridian = False
821  # end if meridian:
822 
823  # somg12 > 1 marks that it needs to be calculated
824  somg12 = 2.0; comg12 = 0.0; omg12 = 0.0
825  if (not meridian and
826  sbet1 == 0 and # and sbet2 == 0
827  # Mimic the way Lambda12 works with calp1 = 0
828  (self.f <= 0 or lon12s >= self.f * 180)):
829 
830  # Geodesic runs along equator
831  calp1 = calp2 = 0.0; salp1 = salp2 = 1.0
832  s12x = self.a * lam12
833  sig12 = omg12 = lam12 / self._f1
834  m12x = self._b * math.sin(sig12)
835  if outmask & Geodesic.GEODESICSCALE:
836  M12 = M21 = math.cos(sig12)
837  a12 = lon12 / self._f1
838 
839  elif not meridian:
840 
841  # Now point1 and point2 belong within a hemisphere bounded by a
842  # meridian and geodesic is neither meridional or equatorial.
843 
844  # Figure a starting point for Newton's method
845  sig12, salp1, calp1, salp2, calp2, dnm = self._InverseStart(
846  sbet1, cbet1, dn1, sbet2, cbet2, dn2, lam12, slam12, clam12, C1a, C2a)
847 
848  if sig12 >= 0:
849  # Short lines (InverseStart sets salp2, calp2, dnm)
850  s12x = sig12 * self._b * dnm
851  m12x = (Math.sq(dnm) * self._b * math.sin(sig12 / dnm))
852  if outmask & Geodesic.GEODESICSCALE:
853  M12 = M21 = math.cos(sig12 / dnm)
854  a12 = math.degrees(sig12)
855  omg12 = lam12 / (self._f1 * dnm)
856  else:
857 
858  # Newton's method. This is a straightforward solution of f(alp1) =
859  # lambda12(alp1) - lam12 = 0 with one wrinkle. f(alp) has exactly one
860  # root in the interval (0, pi) and its derivative is positive at the
861  # root. Thus f(alp) is positive for alp > alp1 and negative for alp <
862  # alp1. During the course of the iteration, a range (alp1a, alp1b) is
863  # maintained which brackets the root and with each evaluation of f(alp)
864  # the range is shrunk if possible. Newton's method is restarted
865  # whenever the derivative of f is negative (because the new value of
866  # alp1 is then further from the solution) or if the new estimate of
867  # alp1 lies outside (0,pi); in this case, the new starting guess is
868  # taken to be (alp1a + alp1b) / 2.
869  # real ssig1, csig1, ssig2, csig2, eps
870  numit = 0
871  tripn = tripb = False
872  # Bracketing range
873  salp1a = Geodesic.tiny_; calp1a = 1.0
874  salp1b = Geodesic.tiny_; calp1b = -1.0
875 
876  while numit < Geodesic.maxit2_:
877  # the WGS84 test set: mean = 1.47, sd = 1.25, max = 16
878  # WGS84 and random input: mean = 2.85, sd = 0.60
879  (v, salp2, calp2, sig12, ssig1, csig1, ssig2, csig2,
880  eps, domg12, dv) = self._Lambda12(
881  sbet1, cbet1, dn1, sbet2, cbet2, dn2,
882  salp1, calp1, slam12, clam12, numit < Geodesic.maxit1_,
883  C1a, C2a, C3a)
884  # 2 * tol0 is approximately 1 ulp for a number in [0, pi].
885  # Reversed test to allow escape with NaNs
886  if tripb or not (abs(v) >= (8 if tripn else 1) * Geodesic.tol0_):
887  break
888  # Update bracketing values
889  if v > 0 and (numit > Geodesic.maxit1_ or
890  calp1/salp1 > calp1b/salp1b):
891  salp1b = salp1; calp1b = calp1
892  elif v < 0 and (numit > Geodesic.maxit1_ or
893  calp1/salp1 < calp1a/salp1a):
894  salp1a = salp1; calp1a = calp1
895 
896  numit += 1
897  if numit < Geodesic.maxit1_ and dv > 0:
898  dalp1 = -v/dv
899  sdalp1 = math.sin(dalp1); cdalp1 = math.cos(dalp1)
900  nsalp1 = salp1 * cdalp1 + calp1 * sdalp1
901  if nsalp1 > 0 and abs(dalp1) < math.pi:
902  calp1 = calp1 * cdalp1 - salp1 * sdalp1
903  salp1 = nsalp1
904  salp1, calp1 = Math.norm(salp1, calp1)
905  # In some regimes we don't get quadratic convergence because
906  # slope -> 0. So use convergence conditions based on epsilon
907  # instead of sqrt(epsilon).
908  tripn = abs(v) <= 16 * Geodesic.tol0_
909  continue
910  # Either dv was not positive or updated value was outside
911  # legal range. Use the midpoint of the bracket as the next
912  # estimate. This mechanism is not needed for the WGS84
913  # ellipsoid, but it does catch problems with more eccentric
914  # ellipsoids. Its efficacy is such for
915  # the WGS84 test set with the starting guess set to alp1 = 90deg:
916  # the WGS84 test set: mean = 5.21, sd = 3.93, max = 24
917  # WGS84 and random input: mean = 4.74, sd = 0.99
918  salp1 = (salp1a + salp1b)/2
919  calp1 = (calp1a + calp1b)/2
920  salp1, calp1 = Math.norm(salp1, calp1)
921  tripn = False
922  tripb = (abs(salp1a - salp1) + (calp1a - calp1) < Geodesic.tolb_ or
923  abs(salp1 - salp1b) + (calp1 - calp1b) < Geodesic.tolb_)
924 
925  lengthmask = (outmask |
926  (Geodesic.DISTANCE
927  if (outmask & (Geodesic.REDUCEDLENGTH |
928  Geodesic.GEODESICSCALE))
929  else Geodesic.EMPTY))
930  s12x, m12x, dummy, M12, M21 = self._Lengths(
931  eps, sig12, ssig1, csig1, dn1, ssig2, csig2, dn2, cbet1, cbet2,
932  lengthmask, C1a, C2a)
933 
934  m12x *= self._b
935  s12x *= self._b
936  a12 = math.degrees(sig12)
937  if outmask & Geodesic.AREA:
938  # omg12 = lam12 - domg12
939  sdomg12 = math.sin(domg12); cdomg12 = math.cos(domg12)
940  somg12 = slam12 * cdomg12 - clam12 * sdomg12
941  comg12 = clam12 * cdomg12 + slam12 * sdomg12
942 
943  # end elif not meridian
944 
945  if outmask & Geodesic.DISTANCE:
946  s12 = 0.0 + s12x # Convert -0 to 0
947 
948  if outmask & Geodesic.REDUCEDLENGTH:
949  m12 = 0.0 + m12x # Convert -0 to 0
950 
951  if outmask & Geodesic.AREA:
952  # From Lambda12: sin(alp1) * cos(bet1) = sin(alp0)
953  salp0 = salp1 * cbet1
954  calp0 = math.hypot(calp1, salp1 * sbet1) # calp0 > 0
955  # real alp12
956  if calp0 != 0 and salp0 != 0:
957  # From Lambda12: tan(bet) = tan(sig) * cos(alp)
958  ssig1 = sbet1; csig1 = calp1 * cbet1
959  ssig2 = sbet2; csig2 = calp2 * cbet2
960  k2 = Math.sq(calp0) * self._ep2
961  eps = k2 / (2 * (1 + math.sqrt(1 + k2)) + k2)
962  # Multiplier = a^2 * e^2 * cos(alpha0) * sin(alpha0).
963  A4 = Math.sq(self.a) * calp0 * salp0 * self._e2
964  ssig1, csig1 = Math.norm(ssig1, csig1)
965  ssig2, csig2 = Math.norm(ssig2, csig2)
966  C4a = list(range(Geodesic.nC4_))
967  self._C4f(eps, C4a)
968  B41 = Geodesic._SinCosSeries(False, ssig1, csig1, C4a)
969  B42 = Geodesic._SinCosSeries(False, ssig2, csig2, C4a)
970  S12 = A4 * (B42 - B41)
971  else:
972  # Avoid problems with indeterminate sig1, sig2 on equator
973  S12 = 0.0
974 
975  if not meridian and somg12 > 1:
976  somg12 = math.sin(omg12); comg12 = math.cos(omg12)
977 
978  if (not meridian and
979  # omg12 < 3/4 * pi
980  comg12 > -0.7071 and # Long difference not too big
981  sbet2 - sbet1 < 1.75): # Lat difference not too big
982  # Use tan(Gamma/2) = tan(omg12/2)
983  # * (tan(bet1/2)+tan(bet2/2))/(1+tan(bet1/2)*tan(bet2/2))
984  # with tan(x/2) = sin(x)/(1+cos(x))
985  domg12 = 1 + comg12; dbet1 = 1 + cbet1; dbet2 = 1 + cbet2
986  alp12 = 2 * math.atan2( somg12 * ( sbet1 * dbet2 + sbet2 * dbet1 ),
987  domg12 * ( sbet1 * sbet2 + dbet1 * dbet2 ) )
988  else:
989  # alp12 = alp2 - alp1, used in atan2 so no need to normalize
990  salp12 = salp2 * calp1 - calp2 * salp1
991  calp12 = calp2 * calp1 + salp2 * salp1
992  # The right thing appears to happen if alp1 = +/-180 and alp2 = 0, viz
993  # salp12 = -0 and alp12 = -180. However this depends on the sign
994  # being attached to 0 correctly. The following ensures the correct
995  # behavior.
996  if salp12 == 0 and calp12 < 0:
997  salp12 = Geodesic.tiny_ * calp1
998  calp12 = -1.0
999  alp12 = math.atan2(salp12, calp12)
1000  S12 += self._c2 * alp12
1001  S12 *= swapp * lonsign * latsign
1002  # Convert -0 to 0
1003  S12 += 0.0
1004 
1005  # Convert calp, salp to azimuth accounting for lonsign, swapp, latsign.
1006  if swapp < 0:
1007  salp2, salp1 = salp1, salp2
1008  calp2, calp1 = calp1, calp2
1009  if outmask & Geodesic.GEODESICSCALE:
1010  M21, M12 = M12, M21
1011 
1012  salp1 *= swapp * lonsign; calp1 *= swapp * latsign
1013  salp2 *= swapp * lonsign; calp2 *= swapp * latsign
1014 
1015  return a12, s12, salp1, calp1, salp2, calp2, m12, M12, M21, S12
1016 
1017  def Inverse(self, lat1, lon1, lat2, lon2,
1018  outmask = GeodesicCapability.STANDARD):
1019  """Solve the inverse geodesic problem
1020 
1021  :param lat1: latitude of the first point in degrees
1022  :param lon1: longitude of the first point in degrees
1023  :param lat2: latitude of the second point in degrees
1024  :param lon2: longitude of the second point in degrees
1025  :param outmask: the :ref:`output mask <outmask>`
1026  :return: a :ref:`dict`
1027 
1028  Compute geodesic between (*lat1*, *lon1*) and (*lat2*, *lon2*).
1029  The default value of *outmask* is STANDARD, i.e., the *lat1*,
1030  *lon1*, *azi1*, *lat2*, *lon2*, *azi2*, *s12*, *a12* entries are
1031  returned.
1032 
1033  """
1034 
1035  a12, s12, salp1,calp1, salp2,calp2, m12, M12, M21, S12 = self._GenInverse(
1036  lat1, lon1, lat2, lon2, outmask)
1037  outmask &= Geodesic.OUT_MASK
1038  if outmask & Geodesic.LONG_UNROLL:
1039  lon12, e = Math.AngDiff(lon1, lon2)
1040  lon2 = (lon1 + lon12) + e
1041  else:
1042  lon2 = Math.AngNormalize(lon2)
1043  result = {'lat1': Math.LatFix(lat1),
1044  'lon1': lon1 if outmask & Geodesic.LONG_UNROLL else
1045  Math.AngNormalize(lon1),
1046  'lat2': Math.LatFix(lat2),
1047  'lon2': lon2}
1048  result['a12'] = a12
1049  if outmask & Geodesic.DISTANCE: result['s12'] = s12
1050  if outmask & Geodesic.AZIMUTH:
1051  result['azi1'] = Math.atan2d(salp1, calp1)
1052  result['azi2'] = Math.atan2d(salp2, calp2)
1053  if outmask & Geodesic.REDUCEDLENGTH: result['m12'] = m12
1054  if outmask & Geodesic.GEODESICSCALE:
1055  result['M12'] = M12; result['M21'] = M21
1056  if outmask & Geodesic.AREA: result['S12'] = S12
1057  return result
1058 
1059  # return a12, lat2, lon2, azi2, s12, m12, M12, M21, S12
1060  def _GenDirect(self, lat1, lon1, azi1, arcmode, s12_a12, outmask):
1061  """Private: General version of direct problem"""
1062  from geographiclib.geodesicline import GeodesicLine
1063  # Automatically supply DISTANCE_IN if necessary
1064  if not arcmode: outmask |= Geodesic.DISTANCE_IN
1065  line = GeodesicLine(self, lat1, lon1, azi1, outmask)
1066  return line._GenPosition(arcmode, s12_a12, outmask)
1067 
1068  def Direct(self, lat1, lon1, azi1, s12,
1069  outmask = GeodesicCapability.STANDARD):
1070  """Solve the direct geodesic problem
1071 
1072  :param lat1: latitude of the first point in degrees
1073  :param lon1: longitude of the first point in degrees
1074  :param azi1: azimuth at the first point in degrees
1075  :param s12: the distance from the first point to the second in
1076  meters
1077  :param outmask: the :ref:`output mask <outmask>`
1078  :return: a :ref:`dict`
1079 
1080  Compute geodesic starting at (*lat1*, *lon1*) with azimuth *azi1*
1081  and length *s12*. The default value of *outmask* is STANDARD, i.e.,
1082  the *lat1*, *lon1*, *azi1*, *lat2*, *lon2*, *azi2*, *s12*, *a12*
1083  entries are returned.
1084 
1085  """
1086 
1087  a12, lat2, lon2, azi2, s12, m12, M12, M21, S12 = self._GenDirect(
1088  lat1, lon1, azi1, False, s12, outmask)
1089  outmask &= Geodesic.OUT_MASK
1090  result = {'lat1': Math.LatFix(lat1),
1091  'lon1': lon1 if outmask & Geodesic.LONG_UNROLL else
1092  Math.AngNormalize(lon1),
1093  'azi1': Math.AngNormalize(azi1),
1094  's12': s12}
1095  result['a12'] = a12
1096  if outmask & Geodesic.LATITUDE: result['lat2'] = lat2
1097  if outmask & Geodesic.LONGITUDE: result['lon2'] = lon2
1098  if outmask & Geodesic.AZIMUTH: result['azi2'] = azi2
1099  if outmask & Geodesic.REDUCEDLENGTH: result['m12'] = m12
1100  if outmask & Geodesic.GEODESICSCALE:
1101  result['M12'] = M12; result['M21'] = M21
1102  if outmask & Geodesic.AREA: result['S12'] = S12
1103  return result
1104 
1105  def ArcDirect(self, lat1, lon1, azi1, a12,
1106  outmask = GeodesicCapability.STANDARD):
1107  """Solve the direct geodesic problem in terms of spherical arc length
1108 
1109  :param lat1: latitude of the first point in degrees
1110  :param lon1: longitude of the first point in degrees
1111  :param azi1: azimuth at the first point in degrees
1112  :param a12: spherical arc length from the first point to the second
1113  in degrees
1114  :param outmask: the :ref:`output mask <outmask>`
1115  :return: a :ref:`dict`
1116 
1117  Compute geodesic starting at (*lat1*, *lon1*) with azimuth *azi1*
1118  and arc length *a12*. The default value of *outmask* is STANDARD,
1119  i.e., the *lat1*, *lon1*, *azi1*, *lat2*, *lon2*, *azi2*, *s12*,
1120  *a12* entries are returned.
1121 
1122  """
1123 
1124  a12, lat2, lon2, azi2, s12, m12, M12, M21, S12 = self._GenDirect(
1125  lat1, lon1, azi1, True, a12, outmask)
1126  outmask &= Geodesic.OUT_MASK
1127  result = {'lat1': Math.LatFix(lat1),
1128  'lon1': lon1 if outmask & Geodesic.LONG_UNROLL else
1129  Math.AngNormalize(lon1),
1130  'azi1': Math.AngNormalize(azi1),
1131  'a12': a12}
1132  if outmask & Geodesic.DISTANCE: result['s12'] = s12
1133  if outmask & Geodesic.LATITUDE: result['lat2'] = lat2
1134  if outmask & Geodesic.LONGITUDE: result['lon2'] = lon2
1135  if outmask & Geodesic.AZIMUTH: result['azi2'] = azi2
1136  if outmask & Geodesic.REDUCEDLENGTH: result['m12'] = m12
1137  if outmask & Geodesic.GEODESICSCALE:
1138  result['M12'] = M12; result['M21'] = M21
1139  if outmask & Geodesic.AREA: result['S12'] = S12
1140  return result
1141 
1142  def Line(self, lat1, lon1, azi1,
1143  caps = GeodesicCapability.STANDARD |
1144  GeodesicCapability.DISTANCE_IN):
1145  """Return a GeodesicLine object
1146 
1147  :param lat1: latitude of the first point in degrees
1148  :param lon1: longitude of the first point in degrees
1149  :param azi1: azimuth at the first point in degrees
1150  :param caps: the :ref:`capabilities <outmask>`
1151  :return: a :class:`~geographiclib.geodesicline.GeodesicLine`
1152 
1153  This allows points along a geodesic starting at (*lat1*, *lon1*),
1154  with azimuth *azi1* to be found. The default value of *caps* is
1155  STANDARD | DISTANCE_IN, allowing direct geodesic problem to be
1156  solved.
1157 
1158  """
1159 
1160  from geographiclib.geodesicline import GeodesicLine
1161  return GeodesicLine(self, lat1, lon1, azi1, caps)
1162 
1163  def _GenDirectLine(self, lat1, lon1, azi1, arcmode, s12_a12,
1164  caps = GeodesicCapability.STANDARD |
1165  GeodesicCapability.DISTANCE_IN):
1166  """Private: general form of DirectLine"""
1167  from geographiclib.geodesicline import GeodesicLine
1168  # Automatically supply DISTANCE_IN if necessary
1169  if not arcmode: caps |= Geodesic.DISTANCE_IN
1170  line = GeodesicLine(self, lat1, lon1, azi1, caps)
1171  if arcmode:
1172  line.SetArc(s12_a12)
1173  else:
1174  line.SetDistance(s12_a12)
1175  return line
1176 
1177  def DirectLine(self, lat1, lon1, azi1, s12,
1178  caps = GeodesicCapability.STANDARD |
1179  GeodesicCapability.DISTANCE_IN):
1180  """Define a GeodesicLine object in terms of the direct geodesic
1181  problem specified in terms of spherical arc length
1182 
1183  :param lat1: latitude of the first point in degrees
1184  :param lon1: longitude of the first point in degrees
1185  :param azi1: azimuth at the first point in degrees
1186  :param s12: the distance from the first point to the second in
1187  meters
1188  :param caps: the :ref:`capabilities <outmask>`
1189  :return: a :class:`~geographiclib.geodesicline.GeodesicLine`
1190 
1191  This function sets point 3 of the GeodesicLine to correspond to
1192  point 2 of the direct geodesic problem. The default value of *caps*
1193  is STANDARD | DISTANCE_IN, allowing direct geodesic problem to be
1194  solved.
1195 
1196  """
1197 
1198  return self._GenDirectLine(lat1, lon1, azi1, False, s12, caps)
1199 
1200  def ArcDirectLine(self, lat1, lon1, azi1, a12,
1201  caps = GeodesicCapability.STANDARD |
1202  GeodesicCapability.DISTANCE_IN):
1203  """Define a GeodesicLine object in terms of the direct geodesic
1204  problem specified in terms of spherical arc length
1205 
1206  :param lat1: latitude of the first point in degrees
1207  :param lon1: longitude of the first point in degrees
1208  :param azi1: azimuth at the first point in degrees
1209  :param a12: spherical arc length from the first point to the second
1210  in degrees
1211  :param caps: the :ref:`capabilities <outmask>`
1212  :return: a :class:`~geographiclib.geodesicline.GeodesicLine`
1213 
1214  This function sets point 3 of the GeodesicLine to correspond to
1215  point 2 of the direct geodesic problem. The default value of *caps*
1216  is STANDARD | DISTANCE_IN, allowing direct geodesic problem to be
1217  solved.
1218 
1219  """
1220 
1221  return self._GenDirectLine(lat1, lon1, azi1, True, a12, caps)
1222 
1223  def InverseLine(self, lat1, lon1, lat2, lon2,
1224  caps = GeodesicCapability.STANDARD |
1225  GeodesicCapability.DISTANCE_IN):
1226  """Define a GeodesicLine object in terms of the invese geodesic problem
1227 
1228  :param lat1: latitude of the first point in degrees
1229  :param lon1: longitude of the first point in degrees
1230  :param lat2: latitude of the second point in degrees
1231  :param lon2: longitude of the second point in degrees
1232  :param caps: the :ref:`capabilities <outmask>`
1233  :return: a :class:`~geographiclib.geodesicline.GeodesicLine`
1234 
1235  This function sets point 3 of the GeodesicLine to correspond to
1236  point 2 of the inverse geodesic problem. The default value of *caps*
1237  is STANDARD | DISTANCE_IN, allowing direct geodesic problem to be
1238  solved.
1239 
1240  """
1241 
1242  from geographiclib.geodesicline import GeodesicLine
1243  a12, _, salp1, calp1, _, _, _, _, _, _ = self._GenInverse(
1244  lat1, lon1, lat2, lon2, 0)
1245  azi1 = Math.atan2d(salp1, calp1)
1246  if caps & (Geodesic.OUT_MASK & Geodesic.DISTANCE_IN):
1247  caps |= Geodesic.DISTANCE
1248  line = GeodesicLine(self, lat1, lon1, azi1, caps, salp1, calp1)
1249  line.SetArc(a12)
1250  return line
1251 
1252  def Polygon(self, polyline = False):
1253  """Return a PolygonArea object
1254 
1255  :param polyline: if True then the object describes a polyline
1256  instead of a polygon
1257  :return: a :class:`~geographiclib.polygonarea.PolygonArea`
1258 
1259  """
1260 
1261  from geographiclib.polygonarea import PolygonArea
1262  return PolygonArea(self, polyline)
1263 
1264  EMPTY = GeodesicCapability.EMPTY
1265  """No capabilities, no output."""
1266  LATITUDE = GeodesicCapability.LATITUDE
1267  """Calculate latitude *lat2*."""
1268  LONGITUDE = GeodesicCapability.LONGITUDE
1269  """Calculate longitude *lon2*."""
1270  AZIMUTH = GeodesicCapability.AZIMUTH
1271  """Calculate azimuths *azi1* and *azi2*."""
1272  DISTANCE = GeodesicCapability.DISTANCE
1273  """Calculate distance *s12*."""
1274  STANDARD = GeodesicCapability.STANDARD
1275  """All of the above."""
1276  DISTANCE_IN = GeodesicCapability.DISTANCE_IN
1277  """Allow distance *s12* to be used as input in the direct geodesic
1278  problem."""
1279  REDUCEDLENGTH = GeodesicCapability.REDUCEDLENGTH
1280  """Calculate reduced length *m12*."""
1281  GEODESICSCALE = GeodesicCapability.GEODESICSCALE
1282  """Calculate geodesic scales *M12* and *M21*."""
1283  AREA = GeodesicCapability.AREA
1284  """Calculate area *S12*."""
1285  ALL = GeodesicCapability.ALL
1286  """All of the above."""
1287  LONG_UNROLL = GeodesicCapability.LONG_UNROLL
1288  """Unroll longitudes, rather than reducing them to the range
1289  [-180d,180d].
1290 
1291  """
1292 
1293 Geodesic.WGS84 = Geodesic(Constants.WGS84_a, Constants.WGS84_f)
1294 """Instantiation for the WGS84 ellipsoid"""
def _Lambda12(self, sbet1, cbet1, dn1, sbet2, cbet2, dn2, salp1, calp1, slam120, clam120, diffp, C1a, C2a, C3a)
Definition: geodesic.py:630
#define max(a, b)
Definition: datatypes.h:20
def Direct(self, lat1, lon1, azi1, s12, outmask=GeodesicCapability.STANDARD)
Definition: geodesic.py:1069
#define min(a, b)
Definition: datatypes.h:19
def InverseLine(self, lat1, lon1, lat2, lon2, caps=GeodesicCapability.STANDARD|GeodesicCapability.DISTANCE_IN)
Definition: geodesic.py:1225
def _GenDirect(self, lat1, lon1, azi1, arcmode, s12_a12, outmask)
Definition: geodesic.py:1060
def _Lengths(self, eps, sig12, ssig1, csig1, dn1, ssig2, csig2, dn2, cbet1, cbet2, outmask, C1a, C2a)
Definition: geodesic.py:432
def ArcDirect(self, lat1, lon1, azi1, a12, outmask=GeodesicCapability.STANDARD)
Definition: geodesic.py:1106
def Polygon(self, polyline=False)
Definition: geodesic.py:1252
def Inverse(self, lat1, lon1, lat2, lon2, outmask=GeodesicCapability.STANDARD)
Definition: geodesic.py:1018
def _InverseStart(self, sbet1, cbet1, dn1, sbet2, cbet2, dn2, lam12, slam12, clam12, C1a, C2a)
Definition: geodesic.py:486
def ArcDirectLine(self, lat1, lon1, azi1, a12, caps=GeodesicCapability.STANDARD|GeodesicCapability.DISTANCE_IN)
Definition: geodesic.py:1202
def _GenInverse(self, lat1, lon1, lat2, lon2, outmask)
Definition: geodesic.py:704
def _C4f(self, eps, c)
Definition: geodesic.py:416
def DirectLine(self, lat1, lon1, azi1, s12, caps=GeodesicCapability.STANDARD|GeodesicCapability.DISTANCE_IN)
Definition: geodesic.py:1179
Definition: pytypes.h:1979
def Line(self, lat1, lon1, azi1, caps=GeodesicCapability.STANDARD|GeodesicCapability.DISTANCE_IN)
Definition: geodesic.py:1144
Double_ range(const Point2_ &p, const Point2_ &q)
size_t len(handle h)
Get the length of a Python object.
Definition: pytypes.h:2244
#define abs(x)
Definition: datatypes.h:17
def _C3f(self, eps, c)
Definition: geodesic.py:404
def _GenDirectLine(self, lat1, lon1, azi1, arcmode, s12_a12, caps=GeodesicCapability.STANDARD|GeodesicCapability.DISTANCE_IN)
Definition: geodesic.py:1165


gtsam
Author(s):
autogenerated on Tue Jul 4 2023 02:34:17