Coverage for pygeodesy/rhumbx.py: 98%
243 statements
« prev ^ index » next coverage.py v7.2.2, created at 2023-10-11 16:04 -0400
« prev ^ index » next coverage.py v7.2.2, created at 2023-10-11 16:04 -0400
2# -*- coding: utf-8 -*-
4u'''A pure Python version of I{Karney}'s C++ classes U{Rhumb
5<https://GeographicLib.SourceForge.io/C++/doc/classGeographicLib_1_1Rhumb.html>} and U{RhumbLine
6<https://GeographicLib.SourceForge.io/C++/doc/classGeographicLib_1_1RhumbLine.html>} from
7I{GeographicLib version 2.0}.
9Class L{RhumbLine} has been enhanced with methods C{intersection2} and C{nearestOn4} to iteratively
10find the intersection of two rhumb lines, respectively the nearest point on a rumb line along a
11geodesic or perpendicular rhumb line from an other point.
13For more details, see the C++ U{GeographicLib<https://GeographicLib.SourceForge.io/C++/doc/index.html>}
14documentation, especially the U{Class List<https://GeographicLib.SourceForge.io/C++/doc/annotated.html>},
15the background information on U{Rhumb lines<https://GeographicLib.SourceForge.io/C++/doc/rhumb.html>},
16the utily U{RhumbSolve<https://GeographicLib.SourceForge.io/C++/doc/RhumbSolve.1.html>} and U{Online
17rhumb line calculations<https://GeographicLib.SourceForge.io/cgi-bin/RhumbSolve>}.
19Copyright (C) U{Charles Karney<mailto:Karney@Alum.MIT.edu>} (2014-2022) and licensed under the MIT/X11
20License. For more information, see the U{GeographicLib<https://GeographicLib.SourceForge.io>} documentation.
21'''
22# make sure int/int division yields float quotient
23from __future__ import division as _; del _ # PYCHOK semicolon
25from pygeodesy.basics import copysign0, neg, _zip
26from pygeodesy.constants import PI_2, _0_0s, _0_0, _0_5, _1_0, \
27 _2_0, _4_0, _720_0, _over, _1_over
28# from pygeodesy.ellipsoids import _EWGS84 # from .karney
29from pygeodesy.errors import itemsorted, RhumbError, _Xorder
30from pygeodesy.fmath import hypot, hypot1
31# from pygeodesy.fsums import fsum1f_ # _MODS
32from pygeodesy.interns import NN, _COMMASPACE_
33from pygeodesy.karney import Caps, _GTuple, _EWGS84
34from pygeodesy.ktm import KTransverseMercator, _Xs, \
35 _AlpCoeffs, _BetCoeffs # PYCHOK used!
36from pygeodesy.lazily import _ALL_DOCS, _ALL_LAZY, _ALL_MODS as _MODS
37from pygeodesy.props import deprecated_method, Property, Property_RO, property_RO
38from pygeodesy.rhumbBase import RhumbBase, RhumbLineBase, Int, pairs, \
39 _update_all_rls
40# from pygeodesy.streprs import pairs # from .rhumbBase
41# from pygeodesy.units import Int # from .rhumbBase
42from pygeodesy.utily import atan1, sincos2_
44from math import asinh, atan, cos, cosh, fabs, radians, sin, sinh, sqrt, tan
46__all__ = _ALL_LAZY.rhumbx
47__version__ = '23.09.29'
50class Rhumb(RhumbBase):
51 '''Class to solve the I{direct} and I{inverse rhumb} problems, based on
52 I{elliptic functions} or I{Krüger} series expansion.
54 @see: The U{Detailed Description<https://GeographicLib.SourceForge.io/C++/doc/
55 classGeographicLib_1_1Rhumb.html>} of I{Karney}'s C++ C{Rhumb Class}.
56 '''
57 _mRA = 6 # see .RAorder
59 def __init__(self, a_earth=_EWGS84, f=None, exact=True, name=NN, **RA_TMorder):
60 '''New C{rhumbx.Rhumb}.
62 @kwarg a_earth: This rhumb's earth model (L{Ellipsoid}, L{Ellipsoid2},
63 L{a_f2Tuple}, L{Datum}, 2-tuple C{(a, f)}) or the
64 (equatorial) radius (C{scalar}).
65 @kwarg f: The ellipsoid's flattening (C{scalar}), iff B{C{a_earth}} is
66 a C{scalar}, ignored otherwise.
67 @kwarg exact: If C{True}, use an addition theorem for elliptic integrals
68 to compute I{Divided differences}, otherwise use the I{Krüger}
69 series expansion (C{bool} or C{None}), see also properties
70 C{exact} and C{TMorder}.
71 @kwarg name: Optional name (C{str}).
72 @kwarg RA_TMorder: Optional keyword arguments B{C{RAorder}} and B{C{TMorder}}
73 to set the respective C{order}, see properties C{RAorder}
74 and C{TMorder} and method C{orders}.
76 @raise RhumbError: Invalid B{C{a_earth}}, B{C{f}} or B{C{RA_TMorder}}.
77 '''
78 RhumbBase.__init__(self, a_earth, f, exact, name)
79 if RA_TMorder:
80 self.orders(**RA_TMorder)
82 @Property_RO
83 def _A2(self): # Conformal2RectifyingCoeffs
84 m = self.TMorder
85 return _Xs(_AlpCoeffs, m, self.ellipsoid), m
87 @Property_RO
88 def _B2(self): # Rectifying2ConformalCoeffs
89 m = self.TMorder
90 return _Xs(_BetCoeffs, m, self.ellipsoid), m
92 def _DConformal2Rectifying(self, x, y): # radians
93 return _1_0 + (_sincosSeries(True, x, y, *self._A2) if self.f else _0_0)
95 @deprecated_method
96 def Direct7(self, lat1, lon1, azi12, s12, outmask=Caps.LATITUDE_LONGITUDE_AREA):
97 '''DEPRECATED, use method L{Rhumb.Direct8}.
99 @return: A I{DEPRECATED} L{Rhumb7Tuple}.
100 '''
101 return self.Direct8(lat1, lon1, azi12, s12, outmask=outmask)._to7Tuple()
103 def _DIsometrict(self, phix, phiy, tphix, tphiy, _Dtan_phix_phiy):
104 E = self.ellipsoid
105 return _Dtan_phix_phiy * _Dasinh(tphix, tphiy) - \
106 _Dsin(phix, phiy) * _DeatanhE(sin(phix), sin(phiy), E)
108 def _DIsometric2Rectifyingd(self, psix, psiy): # degrees
109 if self.exact:
110 E = self.ellipsoid
111 phix, phiy, tphix, tphiy = _Eaux4(E.auxIsometric, psix, psiy)
112 t = _Dtant(phix - phiy, tphix, tphiy)
113 r = _over(self._DRectifyingt( tphix, tphiy, t),
114 self._DIsometrict(phix, phiy, tphix, tphiy, t))
115 else:
116 x, y = radians(psix), radians(psiy)
117 r = self._DConformal2Rectifying(_gd(x), _gd(y)) * _Dgd(x, y)
118 return r
120 def _DRectifyingt(self, tphix, tphiy, _Dtan_phix_phiy):
121 E = self.ellipsoid
122 tbetx = E.f1 * tphix
123 tbety = E.f1 * tphiy
124 return (E.f1 * _Dtan_phix_phiy * E.b * PI_2
125 * _DfEt( tbetx, tbety, self._eF)
126 * _Datan(tbetx, tbety)) / E.L
128 def _DRectifying2Conformal(self, x, y): # radians
129 return _1_0 - (_sincosSeries(True, x, y, *self._B2) if self.f else _0_0)
131 def _DRectifying2Isometricd(self, mux, muy): # degrees
132 E = self.ellipsoid
133 phix, phiy, tphix, tphiy = _Eaux4(E.auxRectifying, mux, muy)
134 if self.exact:
135 t = _Dtant(phix - phiy, tphix, tphiy)
136 r = _over(self._DIsometrict(phix, phiy, tphix, tphiy, t),
137 self._DRectifyingt( tphix, tphiy, t))
138 else:
139 r = self._DRectifying2Conformal(radians(mux), radians(muy)) * \
140 _Dgdinv(E.es_taupf(tphix), E.es_taupf(tphiy))
141 return r
143 @Property_RO
144 def _eF(self):
145 '''(INTERNAL) Get the ellipsoid's elliptic function.
146 '''
147 # .k2 = 0.006739496742276434
148 return self._E._elliptic_e12 # _MODS.elliptic.Elliptic(-self._E._e12)
150 def _Inverse4(self, lon12, r, outmask):
151 '''(INTERNAL) See method C{RhumbBase.Inverse}.
152 '''
153 E, Cs = self.ellipsoid, Caps
154 psi1 = E.auxIsometric(r.lat1)
155 psi2 = E.auxIsometric(r.lat2)
156 psi12 = psi2 - psi1 # degrees
157 if (outmask & Cs.DISTANCE):
158 a = s = hypot(lon12, psi12)
159 if a:
160 a *= self._DIsometric2Rectifyingd(psi2, psi1)
161 s = self._mpd * a # == E._Lpd
162 a = copysign0(a, s)
163 r.set_(a12=a, s12=s)
165 if ((outmask | self._debug) & Cs._DEBUG_INVERSE): # PYCHOK no cover
166 r.set_(a=E.a, f=E.f, f1=E.f1, L=E.L,
167 b=E.b, e=E.e, e2=E.e2, k2=self._eF.k2,
168 lon12=lon12, psi1=psi1, exact=self.exact,
169 psi12=psi12, psi2=psi2)
170 return lon12, psi12, psi1, psi2
172 @deprecated_method
173 def Inverse7(self, lat1, lon1, azi12, s12, outmask=Caps.AZIMUTH_DISTANCE_AREA):
174 '''DEPRECATED, use method L{Rhumb.Inverse8}.
176 @return: A I{DEPRECATED} L{Rhumb7Tuple}.
177 '''
178 return self.Inverse8(lat1, lon1, azi12, s12, outmask=outmask)._to7Tuple()
180 @Property_RO
181 def _mpd(self): # meter per degree
182 return self.ellipsoid._Lpd
184 @deprecated_method
185 def orders(self, RAorder=None, TMorder=None): # PYCHOK expected
186 '''DEPRECATED, use properties C{RAorder} and/or C{TMorder}.
188 Get and set the I{RAorder} and/or I{TMorder}.
190 @kwarg RAorder: I{Rhumb Area} order (C{int}, 4, 5, 6, 7
191 or 8).
192 @kwarg TMorder: I{Transverse Mercator} order (C{int}, 4,
193 5, 6, 7 or 8).
195 @return: L{RhumbOrder2Tuple}C{(RAorder, TMorder)} with
196 the previous C{RAorder} and C{TMorder} setting.
197 '''
198 t = RhumbOrder2Tuple(self.RAorder, self.TMorder)
199 if RAorder not in (None, t.RAorder): # PYCHOK attr
200 self.RAorder = RAorder
201 if TMorder not in (None, t.TMorder): # PYCHOK attr
202 self.TMorder = TMorder
203 return t
205 @Property_RO
206 def _RA2(self):
207 # for WGS84: (0, -0.0005583633519275459, -3.743803759172812e-07, -4.633682270824446e-10,
208 # RAorder 6: -7.709197397676237e-13, -1.5323287106694307e-15, -3.462875359099873e-18)
209 m = self.RAorder
210 return _Xs(_RACoeffs, m, self.ellipsoid, RA=True), m
212 @Property
213 def RAorder(self):
214 '''Get the I{Rhumb Area} order (C{int}, 4, 5, 6, 7 or 8).
215 '''
216 return self._mRA
218 @RAorder.setter # PYCHOK setter!
219 def RAorder(self, order):
220 '''Set the I{Rhumb Area} order (C{int}, 4, 5, 6, 7 or 8).
221 '''
222 m = _Xorder(_RACoeffs, RhumbError, RAorder=order)
223 if self._mRA != m:
224 _update_all_rls(self)
225 self._mRA = m
227# _RhumbLine = RhumbLine # see further below
229 def _S12d(self, psi1, psi2, lon12): # degrees
230 '''(INTERNAL) Compute the area C{S12}.
231 '''
232 S = (self.ellipsoid.areax if self.exact else
233 self.ellipsoid.area) * lon12 / _720_0
234 if S:
235 x, y = radians(psi1), radians(psi2) # _meanSinXi(x, y)
236 s = _Dlog(cosh(x), cosh(y)) * _Dcosh(x, y)
237 if self.f:
238 s += _sincosSeries(False, _gd(x), _gd(y), *self._RA2) * _Dgd(x, y)
239 S *= s
240 return S
242 def toStr(self, prec=6, sep=_COMMASPACE_, **unused): # PYCHOK signature
243 '''Return this C{Rhumb} as string.
245 @kwarg prec: The C{float} precision, number of decimal digits (0..9).
246 Trailing zero decimals are stripped for B{C{prec}} values
247 of 1 and above, but kept for negative B{C{prec}} values.
248 @kwarg sep: Separator to join (C{str}).
250 @return: Tuple items (C{str}).
251 '''
252 d = dict(ellipsoid=self.ellipsoid, RAorder=self.RAorder,
253 exact=self.exact, TMorder=self.TMorder)
254 return sep.join(pairs(itemsorted(d, asorted=False), prec=prec))
257class RhumbLine(RhumbLineBase):
258 '''Compute one or several points on a single rhumb line.
260 Class C{RhumbLine} facilitates the determination of points on
261 a single rhumb line. The starting point (C{lat1}, C{lon1})
262 and the azimuth C{azi12} are specified once.
264 Method C{RhumbLine.Position} returns the location of an other
265 point at distance C{s12} along and the area C{S12} under the
266 rhumb line.
268 Method C{RhumbLine.intersection2} finds the intersection between
269 two rhumb lines.
271 Method C{RhumbLine.nearestOn4} computes the nearest point on and
272 the distance to a rhumb line in different ways.
273 '''
274 _Rhumb = Rhumb # rhumbx.Rhumb
276 def __init__(self, rhumb, lat1=0, lon1=0, azi12=None, **caps_name): # PYCHOK signature
277 '''New C{rhumbx.RhumbLine}.
279 @arg rhumb: The rhumb reference (C{rhumbx.Rhumb}).
280 @kwarg lat1: Latitude of the start point (C{degrees90}).
281 @kwarg lon1: Longitude of the start point (C{degrees180}).
282 @kwarg azi12: Azimuth of this rhumb line (compass C{degrees}).
283 @kwarg caps_name: Optional keyword arguments C{B{name}=NN} and
284 C{B{caps}=0}, a bit-or'ed combination of L{Caps}
285 values specifying the required capabilities. Include
286 C{Caps.LINE_OFF} if updates to the B{C{rhumb}} should
287 I{not} be reflected in this rhumb line.
288 '''
289 RhumbLineBase.__init__(self, rhumb, lat1, lon1, azi12, **caps_name)
291 @Property_RO
292 def _dpm12(self): # PYCHOK no cover
293 '''(INTERNAL) Get this rhumb line's parallel I{circle radius} (C{degree per meter}).
294 '''
295 r = self._salp
296 if r:
297 r = _over(r, self.ellipsoid.circle4(self.lat1).radius)
298 return r
300 @Property_RO
301 def _mu1(self):
302 '''(INTERNAL) Get the I{rectifying auxiliary} latitude (C{degrees}).
303 '''
304 return self.ellipsoid.auxRectifying(self.lat1)
306 def _mu2lat(self, mu):
307 '''(INTERNAL) Get the inverse I{rectifying auxiliary} latitude (C{degrees}).
308 '''
309 return self.ellipsoid.auxRectifying(mu, inverse=True)
311 def _Position4(self, unused, mu2, s12, mu12):
312 '''(INTERNAL) See method C{RhumbLineBase._Position}.
313 '''
314 psi1 = psi2 = self._psi1
315 if mu12: # self._calp != 0
316 lat2 = self._mu2lat(mu2)
317 psi12 = self.rhumb._DRectifying2Isometricd(mu2, self._mu1) * mu12
318 lon2 = self._talp * psi12
319 psi2 += psi12
320 else: # meridional
321 lat2 = self.lat1
322 lon2 = self._dpm12 * s12
323 return lat2, lon2, psi1, psi2
325 @Property_RO
326 def _psi1(self):
327 '''(INTERNAL) Get the I{isometric auxiliary} latitude C{psi} (C{degrees}).
328 '''
329 return self.ellipsoid.auxIsometric(self.lat1)
331 @property_RO
332 def RAorder(self):
333 '''Get this rhumb line's I{Rhumb Area} order (C{int}, 4, 5, 6, 7 or 8).
334 '''
335 return self.rhumb.RAorder
337Rhumb._RhumbLine = RhumbLine # PYCHOK see RhumbBase._RhumbLine
340class RhumbOrder2Tuple(_GTuple):
341 '''2-Tuple C{(RAorder, TMorder)} with a I{Rhumb Area} and
342 I{Transverse Mercator} order, both C{int}, DEPRECATED.
343 '''
344 _Names_ = (Rhumb.RAorder.name, Rhumb.TMorder.name)
345 _Units_ = ( Int, Int)
348# Use I{Divided Differences} to determine (mu2 - mu1) / (psi2 - psi1) accurately.
349# Definition: _Df(x,y,d) = (f(x) - f(y)) / (x - y), @see W. M. Kahan & R. J.
350# Fateman, "Symbolic computation of Divided Differences", SIGSAM Bull. 33(3),
351# 7-28 (1999). U{ACM<https://DL.ACM.org/doi/pdf/10.1145/334714.334716> and @see
352# U{UCB<https://www.CS.Berkeley.edu/~fateman/papers/divdiff.pdf>}, Dec 8, 1999.
354def _Dasinh(x, y):
355 hx = hypot1(x)
356 d = x - y
357 if d:
358 hx *= y
359 hy = x * hypot1(y)
360 t = (d * (x + y) / (hy + hx)) if (x * y) > 0 else (hy - hx)
361 r = asinh(t) / d
362 else:
363 r = _1_0 / hx
364 return r
367def _Datan(x, y):
368 xy = x * y
369 r = xy + _1_0
370 d = x - y
371 if d: # 2 * xy > -1 == 2 * xy + 1 > 0 == xy + r > 0 == xy > -r
372 r = (atan1(d, r) if xy > -r else (atan1(x) - atan1(y))) / d
373 else:
374 r = _1_over(r)
375 return r
378def _Dcosh(x, y):
379 return _Dsincos(x, y, sinh, sinh)
382def _DeatanhE(x, y, E): # see .albers._Datanhee
383 # Deatanhe(x, y) = eatanhe((x - y) / (1 - e^2 * x * y)) / (x - y)
384 e = _1_0 - E.e2 * x * y
385 if e: # assert not isnear0(e)
386 d = x - y
387 e = (E._es_atanh(d / e) / d) if d else (E.e2 / e)
388 return e
391def _DfEt(tx, ty, eF): # tangents
392 # eF = Elliptic(-E.e12) # -E.e2 / (1 - E.e2)
393 r, x, y, = _1_0, atan(tx), atan(ty)
394 d = x - y
395 if (x * y) > 0:
396 # See U{DLMF<https://DLMF.NIST.gov/19.11>}: 19.11.2 and 19.11.4
397 # letting theta -> x, phi -> -y, psi -> z
398 # (E(x) - E(y)) / d = E(z)/d - k2 * sin(x) * sin(y) * sin(z)/d
399 # tan(z/2) = (sin(x)*Delta(y) - sin(y)*Delta(x)) / (cos(x) + cos(y))
400 # = d * Dsin(x,y) * (sin(x) + sin(y))/(cos(x) + cos(y)) /
401 # (sin(x)*Delta(y) + sin(y)*Delta(x))
402 # = t = d * Dt
403 # sin(z) = 2*t/(1+t^2); cos(z) = (1-t^2)/(1+t^2)
404 # Alt (this only works for |z| <= pi/2 -- however, this conditions
405 # holds if x*y > 0):
406 # sin(z) = d * Dsin(x,y) * (sin(x) + sin(y)) /
407 # (sin(x)*cos(y)*Delta(y) + sin(y)*cos(x)*Delta(x))
408 # cos(z) = sqrt((1-sin(z))*(1+sin(z)))
409 sx, cx, sy, cy = sincos2_(x, y)
410 D = (cx + cy) * (eF.fDelta(sy, cy) * sx +
411 eF.fDelta(sx, cx) * sy)
412 D = (sx + sy) * _Dsin(x, y) / D
413 t = D * d
414 t2 = _1_0 + t**2
415 D *= _2_0 / t2
416 s = D * d
417 if s:
418 c = (t + _1_0) * (_1_0 - t) / t2
419 r = eF.fE(s, c, eF.fDelta(s, c)) / s
420 r = D * (r - eF.k2 * sx * sy)
421 elif d:
422 r = (eF.fE(x) - eF.fE(y)) / d
423 return r
426def _Dgd(x, y):
427 return _Datan(sinh(x), sinh(y)) * _Dsinh(x, y)
430def _Dgdinv(x, y): # x, y are tangents
431 return _Dasinh(x, y) / _Datan(x, y)
434def _Dlog(x, y):
435 d = (x - y) * _0_5
436 # Changed atanh(t / (x + y)) to asinh(t / (2 * sqrt(x*y))) to
437 # avoid taking atanh(1) when x is large and y is 1. This also
438 # fixes bogus results being returned for the area when an endpoint
439 # is at a pole. N.B. this routine is invoked with positive x
440 # and y, so the sqrt is always taken of a positive quantity.
441 return (asinh(d / sqrt(x * y)) / d) if d else _1_over(x)
444def _Dsin(x, y):
445 return _Dsincos(x, y, sin, cos)
448def _Dsincos(x, y, sin_, cos_):
449 r = cos_((x + y) * _0_5)
450 d = (x - y) * _0_5
451 if d:
452 r *= sin_(d) / d
453 return r
456def _Dsinh(x, y):
457 return _Dsincos(x, y, sinh, cosh)
460def _Dtan(x, y): # PYCHOK no cover
461 return _Dtant(x - y, tan(x), tan(y))
464def _Dtant(dxy, tx, ty):
465 txy = tx * ty
466 r = txy + _1_0
467 if dxy: # 2 * txy > -1 == 2 * txy + 1 > 0 == txy + r > 0 == txy > -r
468 r = ((tan(dxy) * r) if txy > -r else (tx - ty)) / dxy
469 return r
472def _Eaux4(E_aux, mu_psi_x, mu_psi_y): # degrees
473 # get inverse auxiliary lats in radians and tangents
474 phix = radians(E_aux(mu_psi_x, inverse=True))
475 phiy = radians(E_aux(mu_psi_y, inverse=True))
476 return phix, phiy, tan(phix), tan(phiy)
479def _gd(x):
480 return atan(sinh(x))
483def _sincosSeries(sinp, x, y, C, n):
484 # N.B. C[] has n+1 elements of which
485 # C[0] is ignored and n >= 0
486 # Use Clenshaw summation to evaluate
487 # m = (g(x) + g(y)) / 2 -- mean value
488 # s = (g(x) - g(y)) / (x - y) -- average slope
489 # where
490 # g(x) = sum(C[j] * SC(2 * j * x), j = 1..n)
491 # SC = sinp ? sin : cos
492 # CS = sinp ? cos : sin
493 # ...
494 d, _neg = (x - y), neg
495 sp, cp, sd, cd = sincos2_(x + y, d)
496 sd = (sd / d) if d else _1_0
497 s = _neg(sp * sd) # negative
498 # 2x2 matrices in row-major order
499 a1 = s * d**2
500 a2 = s * _4_0
501 a0 = a3 = _2_0 * cp * cd # m
502 b2 = b1 = _0_0s(4)
503 if n > 0:
504 b1 = C[n], _0_0, _0_0, C[n]
506 _fsum = _MODS.fsums.fsum1f_
507 for j in range(n - 1, 0, -1): # C[0] unused
508 b1, b2, Cj = b2, b1, C[j]
509 # b1 = a * b2 - b1 + C[j] * I
510 m0, m1, m2, m3 = b2
511 n0, n1, n2, n3 = map(_neg, b1)
512 b1 = (_fsum(a0 * m0, a1 * m2, n0, Cj),
513 _fsum(a0 * m1, a1 * m3, n1),
514 _fsum(a2 * m0, a3 * m2, n2),
515 _fsum(a2 * m1, a3 * m3, n3, Cj))
516 # Here are the full expressions for m and s
517 # f01, f02, f11, f12 = (0, 0, cd * sp, 2 * sd * cp) if sinp else \
518 # (1, 0, cd * cp, -2 * sd * sp)
519 # m = -b2[1] * f02 + (C[0] - b2[0]) * f01 + b1[0] * f11 + b1[1] * f12
520 # s = -b2[2] * f01 + (C[0] - b2[3]) * f02 + b1[2] * f11 + b1[3] * f12
521 cd *= b1[2]
522 sd *= b1[3] * _2_0
523 s = _fsum(cd * sp, sd * cp) if sinp else \
524 _fsum(cd * cp, _neg(sd * sp), _neg(b2[2]))
525 return s
528_RACoeffs = { # Generated by Maxima on 2015-05-15 08:24:04-04:00
529 4: ( # GEOGRAPHICLIB_RHUMBAREA_ORDER == 4
530 691, 7860, -20160, 18900, 0, 56700, # R[0]/n^0, polynomial(n), order 4
531 1772, -5340, 6930, -4725, 14175, # R[1]/n^1, polynomial(n), order 3
532 -1747, 1590, -630, 4725, # PYCHOK R[2]/n^2, polynomial(n), order 2
533 104, -31, 315, # R[3]/n^3, polynomial(n), order 1
534 -41, 420), # PYCHOK R[4]/n^4, polynomial(n), order 0, count = 20
535 5: ( # GEOGRAPHICLIB_RHUMBAREA_ORDER == 5
536 -79036, 22803, 259380, -665280, 623700, 0, 1871100, # PYCHOK R[0]/n^0, polynomial(n), order 5
537 41662, 58476, -176220, 228690, -155925, 467775, # PYCHOK R[1]/n^1, polynomial(n), order 4
538 18118, -57651, 52470, -20790, 155925, # PYCHOK R[2]/n^2, polynomial(n), order 3
539 -23011, 17160, -5115, 51975, # PYCHOK R[3]/n^3, polynomial(n), order 2
540 5480, -1353, 13860, # PYCHOK R[4]/n^4, polynomial(n), order 1
541 -668, 5775), # PYCHOK R[5]/n^5, polynomial(n), order 0, count = 27
542 6: ( # GEOGRAPHICLIB_RHUMBAREA_ORDER == 6
543 128346268, -107884140, 31126095, 354053700, -908107200, 851350500, 0, 2554051500, # R[0]/n^0, polynomial(n), order 6
544 -114456994, 56868630, 79819740, -240540300, 312161850, -212837625, 638512875, # PYCHOK R[1]/n^1, polynomial(n), order 5
545 51304574, 24731070, -78693615, 71621550, -28378350, 212837625, # R[2]/n^2, polynomial(n), order 4
546 1554472, -6282003, 4684680, -1396395, 14189175, # R[3]/n^3, polynomial(n), order 3
547 -4913956, 3205800, -791505, 8108100, # PYCHOK R[4]/n^4, polynomial(n), order 2
548 1092376, -234468, 2027025, # R[5]/n^5, polynomial(n), order 1
549 -313076, 2027025), # PYCHOK R[6]/n^6, polynomial(n), order 0, count = 35
550 7: ( # GEOGRAPHICLIB_RHUMBAREA_ORDER == 7
551 -317195588, 385038804, -323652420, 93378285, 1062161100, -2724321600, 2554051500, 0, 7662154500, # PYCHOK R[0]/n^0, polynomial(n), order 7
552 258618446, -343370982, 170605890, 239459220, -721620900, 936485550, -638512875, 1915538625, # PYCHOK R[1]/n^1, polynomial(n), order 6
553 -248174686, 153913722, 74193210, -236080845, 214864650, -85135050, 638512875, # PYCHOK R[2]/n^2, polynomial(n), order 5
554 114450437, 23317080, -94230045, 70270200, -20945925, 212837625, # PYCHOK R[3]/n^3, polynomial(n), order 4
555 15445736, -103193076, 67321800, -16621605, 170270100, # PYCHOK R[4]/n^4, polynomial(n), order 3
556 -27766753, 16385640, -3517020, 30405375, # PYCHOK R[4]/n^4, polynomial(n), order 3
557 4892722, -939228, 6081075, # PYCHOK R[4]/n^4, polynomial(n), order 3
558 -3189007, 14189175), # PYCHOK R[7]/n^7, polynomial(n), order 0, count = 44
559 8: ( # GEOGRAPHICLIB_RHUMBAREA_ORDER == 8
560 71374704821, -161769749880, 196369790040, -165062734200, 47622925350, 541702161000, -1389404016000, 1302566265000, 0, 3907698795000, # R[0]/n^0, polynomial(n), order 8
561 -13691187484, 65947703730, -87559600410, 43504501950, 61062101100, -184013329500, 238803815250, -162820783125, 488462349375, # PYCHOK R[1]/n^1, polynomial(n), order 7
562 30802104839, -63284544930, 39247999110, 18919268550, -60200615475, 54790485750, -21709437750, 162820783125, # R[2]/n^2, polynomial(n), order 6
563 -8934064508, 5836972287, 1189171080, -4805732295, 3583780200, -1068242175, 10854718875, # PYCHOK R[3]/n^3, polynomial(n), order 5
564 50072287748, 3938662680, -26314234380, 17167059000, -4238509275, 43418875500, # R[4]/n^4, polynomial(n), order 4
565 359094172, -9912730821, 5849673480, -1255576140, 10854718875, # R[5]/n^5, polynomial(n), order 3
566 -16053944387, 8733508770, -1676521980, 10854718875, # PYCHOK R[6]/n^6, polynomial(n), order 2
567 930092876, -162639357, 723647925, # R[7]/n^7, polynomial(n), order 1
568 -673429061, 1929727800) # PYCHOK R[8]/n^8, polynomial(n), order 0, count = 54
569}
571__all__ += _ALL_DOCS(Caps, Rhumb, RhumbLine)
573if __name__ == '__main__':
575 from pygeodesy.lazily import printf
577 def _re(fmt, r3, x3):
578 e3 = []
579 for r, x in _zip(r3, x3): # strict=True
580 e = fabs(r - x) / fabs(x)
581 e3.append('%.g' % (e,))
582 printf((fmt % r3) + ' rel errors: ' + ', '.join(e3))
584 # <https://GeographicLib.SourceForge.io/cgi-bin/RhumbSolve> version 2.0
585 rhumb = Rhumb(exact=True) # WGS84 default
586 printf('# %r\n', rhumb)
587 r = rhumb.Direct8(40.6, -73.8, 51, 5.5e6) # from JFK about NE
588 _re('# JFK NE lat2=%.8f, lon2=%.8f, S12=%.1f', (r.lat2, r.lon2, r.S12), (71.68889988, 0.25551982, 44095641862956.148438))
589 r = rhumb.Inverse8(40.6, -73.8, 51.6, -0.5) # JFK to LHR
590 _re('# JFK-LHR azi12=%.8f, s12=%.3f S12=%.1f', (r.azi12, r.s12, r.S12), (77.76838971, 5771083.383328, 37395209100030.367188))
591 r = rhumb.Inverse8(40.6, -73.8, 35.8, 140.3) # JFK to Tokyo Narita
592 _re('# JFK-NRT azi12=%.8f, s12=%.3f S12=%.1f', (r.azi12, r.s12, r.S12), (-92.388887981699639, 12782581.0676841792, -63760642939072.492))
594# % python3 -m pygeodesy.rhumbx
596# Rhumb(RAorder=6, TMorder=6, ellipsoid=Ellipsoid(name='WGS84', a=6378137, b=6356752.31424518, f_=298.25722356, f=0.00335281, f2=0.00336409, n=0.00167922, e=0.08181919, e2=0.00669438, e21=0.99330562, e22=0.0067395, e32=0.00335843, A=6367449.14582341, L=10001965.72931272, R1=6371008.77141506, R2=6371007.18091847, R3=6371000.79000916, Rbiaxial=6367453.63451633, Rtriaxial=6372797.5559594), exact=True)
598# JFK NE lat2=71.68889988, lon2=0.25551982, S12=44095641862956.1 rel errors: 4e-11, 2e-08, 5e-16
599# JFK-LHR azi12=77.76838971, s12=5771083.383 S12=37395209100030.4 rel errors: 3e-12, 5e-15, 0
600# JFK-NRT azi12=-92.38888798, s12=12782581.068 S12=-63760642939072.5 rel errors: 2e-16, 3e-16, 0
602# **) MIT License
603#
604# Copyright (C) 2022-2023 -- mrJean1 at Gmail -- All Rights Reserved.
605#
606# Permission is hereby granted, free of charge, to any person obtaining a
607# copy of this software and associated documentation files (the "Software"),
608# to deal in the Software without restriction, including without limitation
609# the rights to use, copy, modify, merge, publish, distribute, sublicense,
610# and/or sell copies of the Software, and to permit persons to whom the
611# Software is furnished to do so, subject to the following conditions:
612#
613# The above copyright notice and this permission notice shall be included
614# in all copies or substantial portions of the Software.
615#
616# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
617# OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
618# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
619# THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR
620# OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE,
621# ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR
622# OTHER DEALINGS IN THE SOFTWARE.