Coverage for pygeodesy/auxilats/auxDLat.py: 95%
148 statements
« prev ^ index » next coverage.py v7.2.2, created at 2023-08-07 07:28 -0400
« prev ^ index » next coverage.py v7.2.2, created at 2023-08-07 07:28 -0400
1# -*- coding: utf-8 -*-
3u'''Class L{AuxDLat} transcoded to Python from I{Karney}'s C++ class U{DAuxLatitude
4<https://GeographicLib.SourceForge.io/C++/doc/classGeographicLib_1_1DAuxLatitude.html>}
5in I{GeographicLib version 2.2+}.
7Copyright (C) U{Charles Karney<mailto:Charles@Karney.com>} (2022-2023) and licensed
8under the MIT/X11 License. For more information, see the U{GeographicLib
9<https://GeographicLib.SourceForge.io>} documentation.
10'''
11# make sure int/int division yields float quotient, see .basics
12from __future__ import division as _; del _ # PYCHOK semicolon
14from pygeodesy.auxilats.auxily import Aux, _Datan, _Dasinh, _sc, _sn, AuxError
15from pygeodesy.auxilats.auxLat import AuxLat, _ALL_DOCS
16from pygeodesy.basics import map1, _reverange
17from pygeodesy.constants import INF, NAN, isfinite, isinf, isnan, _over, \
18 _0_0, _0_5, _1_0, _2_0, _N_2_0, _3_0
19from pygeodesy.elliptic import Elliptic as _Ef, Fsum
20# from pygeodesy.errors import AuxError # from .auxilats.auxily
21# from pygeodesy.fsums import Fsum # from .elliptic
22# from pygeodesy.lazily import _ALL_DOCS # from .auxilats.auxLat
24from math import atan, atan2, cos, sin, sqrt
26__all__ = ()
27__version__ = '23.08.05'
30class AuxDLat(AuxLat):
31 '''Class to compute C{Divided Differences} of I{Auxiliary}
32 latitudes and other C{Divided Differences} needed for
33 L{RhumbAux} and L{RhumbLineAux} calculations.
34 '''
36 def _Datanhee(self, x, y):
37 # atan(e*sn(tphi))/e:
38 # Datan(e*sn(x),e*sn(y))*Dsn(x,y)/Datan(x,y)
39 # asinh(e1*sn(fm1*tphi)):
40 # Dasinh(e1*sn(fm1*x)), e1*sn(fm1*y)) *
41 # e1 * Dsn(fm1*x, fm1*y) *fm1 / (e * Datan(x,y))
42 # = Dasinh(e1*sn(fm1*x)), e1*sn(fm1*y)) *
43 # Dsn(fm1*x, fm1*y) / Datan(x,y)
44 if self.f < 0:
45 e = self._e
46 r = _Datan(e * _sn(x), e * _sn(y))
47 else:
48 x *= self._fm1
49 y *= self._fm1
50 e1 = self._e1
51 r = _Dasinh(e1 * _sn(x), e1 * _sn(y))
52 return _Dsn(x, y) * r
54 def Dconvert(self, auxout, Zeta1, Zeta2):
55 '''I{Divided Difference} of one auxiliary latitude wrt another.
56 '''
57 auxin = Zeta1._AUX
58 # assert Zeta2._AUX == auxin
59 try:
60 if auxin != auxout:
61 cs = self._coeffs(auxout, auxin)
62 # assert len(cs) == self.ALorder
63 r = _DClenshaw(True, Zeta1, Zeta2, cs, self.ALorder)
64 else:
65 r = _1_0
66 except AuxError: # no _coeffs
67 r = NAN
68 return r
70 def DE(self, X, Y):
71 # We assume that X and Y are in [-90d, 90d] and
72 # have the same sign. If not we would include
73 # if (Xn.y() * Yn.y() < 0)
74 # return d != 0 ? (E(X) - E(Y)) / d : 1
75 # The general formula fails for x = y = 0d and
76 # x = y = 90d. Probably this is fixable (the
77 # formula works for other x = y. But let's
78 # also stipulate that x != y.
80 # Make both y positive, so we can do the swap a <-> b trick
81 sx, cx, x = X._yxr_normalized(True)
82 sy, cy, y = Y._yxr_normalized(True)
83 Dt, k2, d = _0_0, -self._e12, (y - x)
84 # Switch prolate to oblate, then use formulas for k2 < 0
85 if self.f < 0: # XXX and False?
86 sx, cx = cx, sx
87 sy, cy = cy, sy
88 d, k2 = -d, self._e2
89 # See DLMF: Eqs (19.11.2) and (19.11.4) letting
90 if sx and sy:
91 t = _sxk2y(sx, sy, k2) + _sxk2y(sy, sx, k2)
92 Dt = _over(_Dsin(x, y) * (sx + sy), t * (cx + cy))
93 t = d * Dt
94 t2 = _1_0 + t**2
95 Dt *= _2_0 / t2
96 sk2 = (d * Dt)**2 * k2
97 d2 = _1_0 - sk2
98 c2 = ((_1_0 - t) * (_1_0 + t) / t2)**2 if t else _1_0
99 # E(z)/sin(z)
100 E_s = (_Ef.fRF(c2, d2, _1_0) -
101 _Ef.fRD(c2, d2, _1_0, _3_0) * sk2)
102 Dt *= E_s - k2 * sx * sy
103 return Dt
105 def DIsometric(self, Phi1, Phi2):
106 '''I{Divided Difference} of the isometric wrt the geographic latitude.
107 '''
108 tx, ty = Phi1.tan, Phi2.tan
109 if isnan(ty) or isnan(tx): # PYCHOK no cover
110 r = NAN
111 elif isinf(ty) or isinf(tx): # PYCHOK no cover
112 r = INF
113 else: # psi = asinh(tan(Phi)) - e^2 * atanhee(tan(Phi))
114 r = self._Datanhee(tx, ty) * self._e2
115 r = _over(_Dasinh(tx, ty) - r, _Datan(tx, ty))
116 return r
118 def DParametric(self, Phi1, Phi2):
119 '''I{Divided Difference} of the parametric wrt the geographic latitude.
120 '''
121 fm1, e2m1 = self._fm1, self._e2m1
122 tx, ty = Phi1.tan, Phi2.tan
123 # DbetaDphi = Datan(fm1*tx, fm1*ty) * fm1 / Datan(tx, ty)
124 # Datan(x, y) = 1/(1 + x^2), for x = y
125 # = (atan(y) - atan(x)) / (y-x), for x*y < 0
126 # = atan( (y-x) / (1 + x*y) ) / (y-x), for x*y > 0
127 txy = tx * ty
128 if txy < 0 or (isinf(ty) and not tx):
129 _a = atan
130 r = _over(_a(fm1 * ty) - _a(fm1 * tx), _a(ty) - _a(tx))
131 elif tx == ty: # includes tx = ty = inf
132 if txy > 1: # == tx**2
133 txy = _1_0 / txy
134 r = txy + e2m1
135 else:
136 r = txy * e2m1 + _1_0
137 r = _over(fm1 * (txy + _1_0), r)
138 else:
139 if txy > 1:
140 tx = _1_0 / tx
141 ty = _1_0 / ty
142 txy = tx * ty
143 t = txy + e2m1
144 else:
145 t = txy * e2m1 + _1_0
146 r = ty - tx
147 r = _over(atan2(r * fm1, t), atan2(r, _1_0 + txy))
148 return r
150 def DParametricZ(self, Zeta1, Zeta2):
151 '''Short for C{.Dconvert(Aux.BETA, Zeta1, Zeta2)}.
152 '''
153 return self.Dconvert(Aux.BETA, Zeta1, Zeta2)
155 def DRectifying(self, Phi1, Phi2):
156 '''I{Divided Difference} of the rectifying wrt the geographic latitude.
157 '''
158 # Stipulate that Phi1 and Phi2 are in [-90d, 90d]
159 x, y = Phi1.toRadians, Phi2.toRadians
160 if y == x: # isnear0
161 Mu1 = self.Rectifying(Phi1, diff=True)
162 tphi1, r = Phi1.tan, Mu1.diff
163 if isfinite(tphi1):
164 r *= _over(_sc(tphi1), _sc(Mu1.tan))**2
165 else: # PYCHOK no cover
166 r = _over(_1_0, r)
167 elif (x * y) < 0:
168 r = _over(self.Rectifying(Phi2).toRadians -
169 self.Rectifying(Phi1).toRadians, y - x)
170 else:
171 r = _over(self.b, self.RectifyingRadius(True))
172 r *= self.DE(*map1(self.Parametric, Phi1, Phi2))
173 r *= self.DParametric(Phi1, Phi2)
174 return r # or INF or NAN
176 def DRectifyingZ(self, Zeta1, Zeta2):
177 '''Short for C{.Dconvert(Aux.MU, Zeta1, Zeta2)}.
178 '''
179 return self.Dconvert(Aux.MU, Zeta1, Zeta2)
182def _DClenshaw(sinp, Zeta1, Zeta2, cs, K):
183 '''(INTERNAL) I{Divided Difference} of L{AuxLat._Clenshaw}.
185 @return: C{Fsum} if sinp otherwise a C{float}.
186 '''
187 s1, c1, r1 = Zeta1._yxr_normalized(False)
188 s2, c2, r2 = Zeta2._yxr_normalized(False)
189 Delta = r2 - r1
190 # Evaluate (Clenshaw(sinp, szeta2, czeta2, cs, K) -
191 # Clenshaw(sinp, szeta1, czeta1, cs, K)) / Delta
192 # or f = sin if sinp else cos
193 # sum(cs[k] * (f((2*k+2) * Zeta2) -
194 # f((2*k+2) * Zeta2))) / Delta
195 #
196 # Delta is EITHER 1, giving the plain difference OR (Zeta2 - Zeta1)
197 # in radians, giving the I{Divided Difference}. Other values will
198 # produce nonsense.
199 #
200 # Suffices a and b denote [1,1], [2,1] elements of matrix/vector
201 cp = cm = c2 * c1
202 t = s2 * s1
203 cp -= t # not +
204 cm += t # not -
206 sp = s2 * c1
207 t = c2 * s1
208 smd = ((sin(Delta) / Delta) if Delta != _1_0 else
209 (sp - t)) if Delta else _1_0
210 sp += t
212 xa = cp * cm * _2_0
213 xb = sp * smd * _N_2_0
214 xD = xb * Delta**2
216 if isfinite(xD) and isfinite(xb) and isfinite(xa):
217 U0a, U1a = Fsum(), Fsum()
218 U0b, U1b = Fsum(), Fsum()
219 else: # XXX avoid Fsum(NAN) exceptions
220 U0a = U1a = U0b = U1b = _0_0
221 for k in _reverange(K): # assert len(cs) == K
222 # t = x . U0 - U1 + cs[k] * I
223 U1a -= U0a * xa + U0b * xD + cs[k]
224 U1b -= U0a * xb + U0b * xa
225 U1a, U0a = U0a, -U1a
226 U1b, U0b = U0b, -U1b
227 # F0a = (sp if sinp else cp) * cm
228 # F0b = (cp if sinp else -sp) * smd
229 # Fm1a = 0 if sinp else 1 # Fm1b = 0
230 # return (U0b * F0a + U0a * F0b - U1b * Fm1a) * 2
231 if sinp:
232 U1b = _0_0
233 else:
234 sp, cp = cp, -sp
235 U0b *= sp * cm
236 U0a *= cp * smd
237 U0a += U0b
238 U0a -= U1b
239 U0a *= _2_0
240 return float(U0a) if sinp else U0a # Fsum
243def _Dsin(x, y): # see also .rhumbx._Dsin
244 r = cos((x + y) * _0_5)
245 d = (x - y) * _0_5
246 if d:
247 r *= sin(d) / d
248 return r
251def _Dsn(x, y):
252 # (sn(y) - sn(x)) / (y - x)
253 if x != y:
254 snx, sny = map1(_sn, x, y)
255 if (x * y) > 0:
256 scx, scy = map1(_sc, x, y)
257 r = _over((snx / scy) + (sny / scx),
258 (snx + sny) * scy * scx)
259 else:
260 r = (sny - snx) / (y - x)
261 elif x:
262 r = _1_0 / (_sc(x) * (x**2 + _1_0)) # == 1 / sqrt3(x**2 + 1)
263 else:
264 r = _1_0
265 return r
268def _sxk2y(sx, sy, k2):
269 # .DE helper
270 if sx:
271 try:
272 sx *= sqrt(_1_0 - sy**2 * k2)
273 except ValueError: # domain error
274 sx = NAN
275 return sx
278__all__ += _ALL_DOCS(AuxDLat)
280# **) MIT License
281#
282# Copyright (C) 2023-2023 -- mrJean1 at Gmail -- All Rights Reserved.
283#
284# Permission is hereby granted, free of charge, to any person obtaining a
285# copy of this software and associated documentation files (the "Software"),
286# to deal in the Software without restriction, including without limitation
287# the rights to use, copy, modify, merge, publish, distribute, sublicense,
288# and/or sell copies of the Software, and to permit persons to whom the
289# Software is furnished to do so, subject to the following conditions:
290#
291# The above copyright notice and this permission notice shall be included
292# in all copies or substantial portions of the Software.
293#
294# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
295# OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
296# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
297# THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR
298# OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE,
299# ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR
300# OTHER DEALINGS IN THE SOFTWARE.