Coverage for pygeodesy/auxilats/auxily.py: 92%
114 statements
« prev ^ index » next coverage.py v7.2.2, created at 2023-10-04 12:08 -0400
« prev ^ index » next coverage.py v7.2.2, created at 2023-10-04 12:08 -0400
2# -*- coding: utf-8 -*-
4u'''(INTERNAL) Private I{Auxiliary} base classes, constants and functions.
6Class L{AuxAngle} transcoded to Python from I{Karney}'s C++ class U{AuxAngle
7<https://GeographicLib.SourceForge.io/C++/doc/classGeographicLib_1_1AuxAngle.html>}
8in I{GeographicLib version 2.2+}.
10Copyright (C) U{Charles Karney<mailto:Karney@Alum.MIT.edu>} (2022-2023) and licensed
11under the MIT/X11 License. For more information, see the U{GeographicLib
12<https://GeographicLib.SourceForge.io>} documentation.
13'''
14# make sure int/int division yields float quotient, see .basics
15from __future__ import division as _; del _ # PYCHOK semicolon
17from pygeodesy.constants import INF, NAN, isinf, isnan, _0_0, _0_5, \
18 _1_0, _copysign_1_0, _over, _1_over
19from pygeodesy.errors import AuxError, NN
20from pygeodesy.fmath import hypot1 as _sc, hypot2_
21# from pygeodesy.interns import NN # from .errors
22from pygeodesy.karney import _ALL_DOCS, _Dict, _MODS # PYCHOK used!
23# from pygeodesy.lazily import _ALL_DOCS, _ALL_MODS as _MODS # from .karney
24# from pygeodesy.named import _Dict # from .karney
25from pygeodesy.utily import atan1
27from math import asinh, copysign
29__all__ = ()
30__version__ = '23.09.26'
33class Aux(object):
34 '''Enum-style Aux names.
35 '''
36 GEOGRAPHIC = PHI = GEODETIC = 0
37 PARAMETRIC = BETA = REDUCED = 1
38 GEOCENTRIC = THETA = 2 # all ...
39 RECTIFYING = MU = 3 # use n^2
40 CONFORMAL = CHI = 4 # use n
41 AUTHALIC = XI = 5 # use n
42 N = 6
43 N2 = 36
45 def __index__(self, aux):
46 # throws KeyError, not IndexError
47 return _MODS.auxilats.auxLat._Aux2Greek[aux]
49 def __len__(self):
50 return Aux.N
52 def _1d(self, auxout, auxin):
53 '''Get the 1-d index into N^2 coeffs.
54 '''
55 N = Aux.N
56 if 0 <= auxout < N and 0 <= auxin < N:
57 return N * auxout + auxin
58 raise AuxError(auxout=auxout, auxin=auxin, N=N) # == AssertionError
60 def Greek(self, aux):
61 '''Get an angle's name (C{str}).
62 '''
63 return _Aux2Greek.get(aux, NN)
65 def len(self, ALorder): # PYCHOK no cover
66 aL = ALorder # aka Lmax
67 mu = Aux.MU * (Aux.MU + 1)
68 ns = Aux.N2 - Aux.N
69 return (mu * (aL * (aL + 3) - 2 * (aL // 2)) // 4 +
70 (ns - mu) * (aL * (aL + 1)) // 2)
72 def power(self, auxout, auxin):
73 '''Get the C{convert} exponent (C{int} or C{None}).
74 '''
75 self._1d(auxout, auxin) # validate
76 return (auxout - auxin) if max(auxin, auxout) < Aux.MU else None
78 def use_n2(self, aux):
79 return aux not in (Aux.CHI, Aux.XI)
81Aux = Aux() # PYCHOK singleton
83_Aux2Greek = {Aux.AUTHALIC: 'Xi',
84 Aux.CONFORMAL: 'Chi',
85 Aux.GEOCENTRIC: 'Theta',
86 Aux.GEODETIC: 'Phi', # == .GEOGRAPHIC
87 Aux.PARAMETRIC: 'Beta',
88 Aux.RECTIFYING: 'Mu'}
89_Greek2Aux = dict(map(reversed, _Aux2Greek.items())) # PYCHOK exported
90# _Greek2Aux.update((_g.upper(), _x) for _g, _x in _Greek2Aux.items())
93class _Coeffs(_Dict):
94 '''(INTERNAL) With C{items keys} string-ified.
95 '''
96 def items(self):
97 for n, v in _Dict.items(self):
98 yield str(n), v
101class _Ufloats(dict): # in .auxilats.auxily
102 '''(INTERNAL) "Uniquify" floats.
103 '''
104 n = 0 # total number of floats
106 def __call__(self, *fs):
107 '''Return a tuple of "uniquified" floats.
108 '''
109 self.n += len(fs)
110 _f = self.setdefault
111 return tuple(_f(f, f) for f in map(float, fs)) # PYCHOK as attr
113 def _Coeffs(self, ALorder, coeffs):
114 '''Return C{coeffs} (C{_Coeffs}, I{embellished}).
115 '''
116# if True:
117# n = 0
118# for d in coeffs.values():
119# for t in d.values():
120# n += len(t)
121# assert n == self.n
122 Cx = _Coeffs(coeffs)
123 Cx.set_(ALorder=ALorder, # used in .auxilats.__main__
124 n=self.n, # total number of floats
125 u=len(self.keys())) # unique floats
126 return Cx
129def _Dasinh(x, y):
130 d = y - x
131 if isinf(d): # PYCHOK no cover
132 r = _0_0
133 elif isnan(d): # PYCHOK no cover
134 r = NAN
135 elif d:
136 xy = x * y
137 if xy > 0:
138 hx, hy = _sc(x), _sc(y)
139 if xy < 1:
140 hx, hy = hy, hx
141 else:
142 x = _1_0 / x
143 y = _1_0 / y
144 r = _over(x + y, hx * x + hy * y)
145 r = asinh(r * d) / d
146 else:
147 r = (asinh(y) - asinh(x)) / d
148 else:
149 r = _1_over(_sc(x))
150 return r
153def _Datan(x, y):
154 xy = x * y
155 r = xy + _1_0
156 if isnan(r): # PYCHOK no cover
157 pass
158 elif x == y:
159 r = _1_over(r)
160 elif x > 0 and isinf(xy): # PYCHOK no cover
161 r = _0_0
162 else:
163 d = y - x
164 if (r + xy) > 0:
165 r = atan1(d, r) / d # atan(d / r) / d
166 else:
167 r = (atan1(y) - atan1(x)) / d
168 return r
171def _Dh(x, y):
172 r = x + y
173 if isnan(r):
174 pass # N.B. NAN for inf-inf
175 elif isinf(x): # PYCHOK no cover
176 r = copysign(_0_5, x)
177 elif isinf(y): # PYCHOK no cover
178 r = copysign(_0_5, y)
179 else:
180 snx, sny = _sn(x), _sn(y)
181 dy = sny * y
182 dx = snx * x
183 d = dy + dx
184 if (d * _0_5):
185 if (x * y) > 0:
186 r *= hypot2_(snx / _sc(y), snx * sny,
187 sny / _sc(x)) / (d + d)
188 else:
189 r = _over((dy - dx) * _0_5, y - x)
190 else: # underflow and x == y == d == 0
191 r *= _0_5 # PYCHOK no cover
192 return r
195def _Dlam(x, y): # Chi1.tan, Chi2.tan
196 # I{Divided difference} of the isometric latitude
197 # with respect to the conformal latitude
198 if isnan(x) or isnan(y): # PYCHOK no cover
199 r = NAN
200 elif isinf(x) or isinf(y): # PYCHOK no cover
201 r = INF
202 elif x == y:
203 r = _sc(x)
204 else:
205 r = _over(_Dasinh(x, y), _Datan(x, y))
206 return r
209def _Dm(X, Y, s): # in .auxDLat, .auxDST
210 # Return M{(X - Y) * s}, inplace X
211 X -= Y
212 X *= s
213 return X # Fsum
216def _Dp0Dpsi(x, y): # Chi1.tan, Chi2.tan
217 # I{Divided difference} of the spherical rhumb area
218 # term with respect to the isometric latitude
219 r = x + y
220 if isnan(r): # PYCHOK no cover
221 pass # NAN for inf-inf
222 elif isinf(x): # PYCHOK no cover
223 r = _copysign_1_0(x)
224 elif isinf(y): # PYCHOK no cover
225 r = _copysign_1_0(y)
226 elif x == y:
227 r = _sn(x)
228 else:
229 r = _Dasinh(_h(x), _h(y))
230 r = _over(_Dh(x, y) * r, _Dasinh(x, y))
231 return r
234def _h(tx):
235 '''(INTERNAL) M{h(tan(x)) = tan(x) * sin(x) / 2}
236 '''
237 if tx:
238 tx *= _sn(tx) * _0_5
239 return tx # preserve signed-0
242def _sn(tx):
243 '''(INTERNAL) M{tx / sqrt(tx**2 + 1)}, converting tan to sin.
244 '''
245 if tx:
246 tx = _copysign_1_0(tx) if isinf(tx) else (
247 NAN if isnan(tx) else (tx / _sc(tx)))
248 return tx # preserve signed-0
251__all__ += _ALL_DOCS(Aux.__class__)
253# **) MIT License
254#
255# Copyright (C) 2023-2023 -- mrJean1 at Gmail -- All Rights Reserved.
256#
257# Permission is hereby granted, free of charge, to any person obtaining a
258# copy of this software and associated documentation files (the "Software"),
259# to deal in the Software without restriction, including without limitation
260# the rights to use, copy, modify, merge, publish, distribute, sublicense,
261# and/or sell copies of the Software, and to permit persons to whom the
262# Software is furnished to do so, subject to the following conditions:
263#
264# The above copyright notice and this permission notice shall be included
265# in all copies or substantial portions of the Software.
266#
267# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
268# OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
269# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
270# THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR
271# OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE,
272# ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR
273# OTHER DEALINGS IN THE SOFTWARE.