1

2# -*- coding: utf-8 -*-

3

4 u'''(INTERNAL) Private I{Auxiliary} base classes, constants and functions.

5

6 Class L{AuxAngle} transcoded to Python from I{Karney}'s C++ class U{AuxAngle

7<https://GeographicLib.SourceForge.io/C++/doc/classGeographicLib_1_1AuxAngle.html>}

8 in I{GeographicLib version 2.2+}.

9

10 Copyright (C) U{Charles Karney<mailto:Karney@Alum.MIT.edu>} (2022-2023) and licensed

11 under 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

15 from __future__ import division as _; del _ # PYCHOK semicolon

16

17 from pygeodesy.constants import INF, NAN, isinf, isnan, _0_0, _0_5, \

18 _1_0, _copysign_1_0, _over, _1_over

19 from pygeodesy.errors import AuxError, NN

20 from pygeodesy.fmath import hypot1 as _sc, hypot2_

21# from pygeodesy.interns import NN # from .errors

22 from 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

25

26 from math import asinh, atan, copysign

27

28__all__ = ()

29__version__ = '23.08.20'

30

31

32 class Aux(object):

33 '''Enum-style Aux names.

34 '''

35 GEOGRAPHIC = PHI = GEODETIC = 0

36 PARAMETRIC = BETA = REDUCED = 1

37 GEOCENTRIC = THETA = 2 # all ...

38 RECTIFYING = MU = 3 # use n^2

39 CONFORMAL = CHI = 4 # use n

40 AUTHALIC = XI = 5 # use n

41 N = 6

42 N2 = 36

43

44 def __index__(self, aux):

45 # throws KeyError, not IndexError

46 return _MODS.auxilats.auxLat._Aux2Greek[aux]

47

48 def __len__(self):

49 return Aux.N

50

51 def _1d(self, auxout, auxin):

52 '''Get the 1-d index into N^2 coeffs.

53 '''

54 N = Aux.N

55 if 0 <= auxout < N and 0 <= auxin < N:

56 return N * auxout + auxin

57 raise AuxError(auxout=auxout, auxin=auxin, N=N) # == AssertionError

58

59 def Greek(self, aux):

60 '''Get an angle's name (C{str}).

61 '''

62 return _Aux2Greek.get(aux, NN)

63

64 def len(self, ALorder): # PYCHOK no cover

65 aL = ALorder # aka Lmax

66 mu = Aux.MU * (Aux.MU + 1)

67 ns = Aux.N2 - Aux.N

68 return (mu * (aL * (aL + 3) - 2 * (aL // 2)) // 4 +

69 (ns - mu) * (aL * (aL + 1)) // 2)

70

71 def power(self, auxout, auxin):

72 '''Get the C{convert} exponent (C{int} or C{None}).

73 '''

74 self._1d(auxout, auxin) # validate

75 return (auxout - auxin) if max(auxin, auxout) < Aux.MU else None

76

77 def use_n2(self, aux):

78 return aux not in (Aux.CHI, Aux.XI)

79

80 Aux = Aux() # PYCHOK singleton

81

82_Aux2Greek = {Aux.AUTHALIC: 'Xi',

83 Aux.CONFORMAL: 'Chi',

84 Aux.GEOCENTRIC: 'Theta',

85 Aux.GEODETIC: 'Phi', # == .GEOGRAPHIC

86 Aux.PARAMETRIC: 'Beta',

87 Aux.RECTIFYING: 'Mu'}

88_Greek2Aux = dict(map(reversed, _Aux2Greek.items())) # PYCHOK exported

89

90

91 class _Coeffs(_Dict):

92 '''(INTERNAL) With C{items keys} string-ified.

93 '''

94 def items(self):

95 for n, v in _Dict.items(self):

96 yield str(n), v

97

98

99 class _Ufloats(dict): # in .auxilats.auxily

100 '''(INTERNAL) "Uniquify" floats.

101 '''

102 n = 0 # total number of floats

103

104 def __call__(self, *fs):

105 '''Return a tuple of "uniquified" floats.

106 '''

107 self.n += len(fs)

108 _f = self.setdefault

109 return tuple(_f(f, f) for f in map(float, fs)) # PYCHOK as attr

110

111 def _Coeffs(self, ALorder, coeffs):

112 '''Return C{coeffs} (C{_Coeffs}, I{embellished}).

113 '''

114# if True:

115# n = 0

116# for d in coeffs.values():

117# for t in d.values():

118# n += len(t)

119# assert n == self.n

120 Cx = _Coeffs(coeffs)

121 Cx.set_(ALorder=ALorder, # used in .auxilats.__main__

122 n=self.n, # total number of floats

123 u=len(self.keys())) # unique floats

124 return Cx

125

126

127 def _Dasinh(x, y):

128 d = y - x

129 if isinf(d): # PYCHOK no cover

130 r = _0_0

131 elif isnan(d): # PYCHOK no cover

132 r = NAN

133 elif d:

134 xy = x * y

135 if xy > 0:

136 hx, hy = _sc(x), _sc(y)

137 if xy < 1:

138 hx, hy = hy, hx

139 else:

140 x = _1_0 / x

141 y = _1_0 / y

142 r = _over(x + y, hx * x + hy * y)

143 r = asinh(r * d) / d

144 else:

145 r = (asinh(y) - asinh(x)) / d

146 else:

147 r = _1_over(_sc(x))

148 return r

149

150

151 def _Datan(x, y):

152 xy = x * y

153 r = xy + _1_0

154 if isnan(r): # PYCHOK no cover

155 pass

156 elif x == y:

157 r = _1_over(r)

158 elif x > 0 and isinf(xy): # PYCHOK no cover

159 r = _0_0

160 else:

161 d = y - x

162 if (r + xy) > 0:

163 r = atan(_over(d, r)) / d # atan2(d, r) / d

164 else:

165 r = (atan(y) - atan(x)) / d

166 return r

167

168

169 def _Dh(x, y):

170 r = x + y

171 if isnan(r):

172 pass # N.B. NAN for inf-inf

173 elif isinf(x): # PYCHOK no cover

174 r = copysign(_0_5, x)

175 elif isinf(y): # PYCHOK no cover

176 r = copysign(_0_5, y)

177 else:

178 snx, sny = _sn(x), _sn(y)

179 d = sny * y + snx * x

180 if (d * _0_5):

181 if (x * y) > 0:

182 r *= hypot2_(snx / _sc(y), snx * sny,

183 sny / _sc(x)) / (d + d)

184 else:

185 r = _over(_h(y) - _h(x), y - x)

186 else: # underflow and x == y == d == 0

187 r *= _0_5 # PYCHOK no cover

188 return r

189

190

191 def _Dlam(x, y): # Chi1.tan, Chi2.tan

192 # I{Divided difference} of the isometric latitude

193 # with respect to the conformal latitude

194 if isnan(x) or isnan(y): # PYCHOK no cover

195 r = NAN

196 elif isinf(x) or isinf(y): # PYCHOK no cover

197 r = INF

198 elif x == y:

199 r = _sc(x)

200 else:

201 r = _over(_Dasinh(x, y), _Datan(x, y))

202 return r

203

204

205 def _Dp0Dpsi(x, y): # Chi1.tan, Chi2.tan

206 # I{Divided difference} of the spherical rhumb area

207 # term with respect to the isometric latitude

208 r = x + y

209 if isnan(r): # PYCHOK no cover

210 pass # NAN for inf-inf

211 elif isinf(x): # PYCHOK no cover

212 r = _copysign_1_0(x)

213 elif isinf(y): # PYCHOK no cover

214 r = _copysign_1_0(y)

215 elif x == y:

216 r = _sn(x)

217 else:

218 r = _Dasinh(_h(x), _h(y))

219 r = _over(_Dh(x, y) * r, _Dasinh(x, y))

220 return r

221

222

223 def _h(tx):

224 '''(INTERNAL) M{h(tan(x)) = tan(x) * sin(x) / 2}

225 '''

226 return (_sn(tx) * tx * _0_5) if tx else _0_0

227

228

229 def _sn(tx):

230 '''(INTERNAL) M{tx / sqrt(tx**2 + 1)}, converting tan to sin.

231 '''

232 if tx:

233 tx = _copysign_1_0(tx) if isinf(tx) else (

234 NAN if isnan(tx) else (tx / _sc(tx)))

235 return tx # preserve signed-0

236

237

238__all__ += _ALL_DOCS(Aux.__class__)

239

240# **) MIT License

241#

242# Copyright (C) 2023-2023 -- mrJean1 at Gmail -- All Rights Reserved.

243#

244# Permission is hereby granted, free of charge, to any person obtaining a

245# copy of this software and associated documentation files (the "Software"),

246# to deal in the Software without restriction, including without limitation

247# the rights to use, copy, modify, merge, publish, distribute, sublicense,

248# and/or sell copies of the Software, and to permit persons to whom the

249# Software is furnished to do so, subject to the following conditions:

250#

251# The above copyright notice and this permission notice shall be included

252# in all copies or substantial portions of the Software.

253#

254# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS

255# OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,

256# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL

257# THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR

258# OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE,

259# ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR

260# OTHER DEALINGS IN THE SOFTWARE.