Coverage for pygeodesy/auxilats/auxDLat.py: 95%

151 statements  

« prev     ^ index     » next       coverage.py v7.2.2, created at 2023-08-20 14:00 -0400

1# -*- coding: utf-8 -*- 

2 

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+}. 

6 

7Copyright (C) U{Charles Karney<mailto:Karney@Alum.MIT.edu>} (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 

13 

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, _0_0, \ 

18 _0_5, _1_0, _2_0, _N_2_0, _3_0, _naninf, \ 

19 _over, _1_over 

20from pygeodesy.elliptic import Elliptic as _Ef, Fsum 

21# from pygeodesy.errors import AuxError # from .auxilats.auxily 

22# from pygeodesy.fsums import Fsum # from .elliptic 

23# from pygeodesy.lazily import _ALL_DOCS # from .auxilats.auxLat 

24 

25from math import atan, atan2, cos, sin, sqrt 

26 

27__all__ = () 

28__version__ = '23.08.19' 

29 

30 

31class AuxDLat(AuxLat): 

32 '''Class to compute C{Divided Differences} of I{Auxiliary} 

33 latitudes and other C{Divided Differences} needed for 

34 L{RhumbAux} and L{RhumbLineAux} calculations. 

35 ''' 

36 

37 def CParametric(self, Zeta1, Zeta2): 

38 '''Short for C{.Dconvert(Aux.BETA, B{Zeta1}, B{Zeta2})}. 

39 ''' 

40 return self.Dconvert(Aux.BETA, Zeta1, Zeta2) 

41 

42 def CRectifying(self, Zeta1, Zeta2): 

43 '''Short for C{.Dconvert(Aux.MU, B{Zeta1}, B{Zeta2})}. 

44 ''' 

45 return self.Dconvert(Aux.MU, Zeta1, Zeta2) 

46 

47 def _Datanhee(self, x, y): 

48 # atan(e*sn(tphi))/e: 

49 # Datan(e*sn(x),e*sn(y))*Dsn(x,y)/Datan(x,y) 

50 # asinh(e1*sn(fm1*tphi)): 

51 # Dasinh(e1*sn(fm1*x)), e1*sn(fm1*y)) * 

52 # e1 * Dsn(fm1*x, fm1*y) *fm1 / (e * Datan(x,y)) 

53 # = Dasinh(e1*sn(fm1*x)), e1*sn(fm1*y)) * 

54 # Dsn(fm1*x, fm1*y) / Datan(x,y) 

55 if self.f < 0: 

56 e = self._e 

57 r = _Datan(e * _sn(x), e * _sn(y)) 

58 else: 

59 x *= self._fm1 

60 y *= self._fm1 

61 e1 = self._e1 

62 r = _Dasinh(e1 * _sn(x), e1 * _sn(y)) 

63 return _Dsn(x, y) * r 

64 

65 def Dconvert(self, auxout, Zeta1, Zeta2): 

66 '''I{Divided Difference} of one auxiliary latitude wrt another. 

67 ''' 

68 auxin = Zeta1._AUX 

69 # assert Zeta2._AUX == auxin 

70 try: 

71 if auxout != auxin: 

72 cs = self._coeffs(auxout, auxin) 

73 # assert len(cs) == self.ALorder 

74 r = _DClenshaw(True, Zeta1, Zeta2, cs, self.ALorder) 

75 else: 

76 r = _1_0 

77 except AuxError: # no _coeffs 

78 r = NAN 

79 return r 

80 

81 def DE(self, X, Y): 

82 # We assume that X and Y are in [-90d, 90d] and 

83 # have the same sign. If not we would include 

84 # if (Xn.y() * Yn.y() < 0) 

85 # return d != 0 ? (E(X) - E(Y)) / d : 1 

86 # The general formula fails for x = y = 0d and 

87 # x = y = 90d. Probably this is fixable (the 

88 # formula works for other x = y. But let's 

89 # also stipulate that x != y. 

90 

91 # Make both y positive, so we can do the swap a <-> b trick 

92 sx, cx, x = X._yxr_normalized(True) 

93 sy, cy, y = Y._yxr_normalized(True) 

94 k2, d = -self._e12, (y - x) 

95 # Switch prolate to oblate, then use formulas for k2 < 0 

96 if self.f < 0: # XXX and False? 

97 sx, cx = cx, sx 

98 sy, cy = cy, sy 

99 d, k2 = -d, self._e2 

100 # See DLMF: Eqs (19.11.2) and (19.11.4) letting 

101 Dt = _Dsin(x, y) * (sx + sy) 

102 if Dt: 

103 t = _sxk2y(sx, sy, k2) + _sxk2y(sy, sx, k2) 

104 Dt = _over(Dt, t * (cx + cy)) 

105 t = d * Dt 

106 t2 = _1_0 + t**2 

107 Dt *= _2_0 / t2 

108 sk2 = (d * Dt)**2 * k2 

109 d2 = _1_0 - sk2 

110 c2 = ((_1_0 - t) * (_1_0 + t) / t2)**2 if t else _1_0 

111 # E(z)/sin(z) 

112 E_s = (_Ef.fRF(c2, d2, _1_0) - 

113 _Ef.fRD(c2, d2, _1_0, _3_0) * sk2) 

114 Dt *= E_s - k2 * sx * sy 

115 return Dt 

116 

117 def DIsometric(self, Phi1, Phi2): 

118 '''I{Divided Difference} of the isometric wrt the geographic latitude. 

119 ''' 

120 tx, ty = Phi1.tan, Phi2.tan 

121 if isnan(ty) or isnan(tx): # PYCHOK no cover 

122 r = NAN 

123 elif isinf(ty) or isinf(tx): # PYCHOK no cover 

124 r = INF 

125 else: # psi = asinh(tan(Phi)) - e^2 * atanhee(tan(Phi)) 

126 r = self._Datanhee(tx, ty) * self._e2 

127 r = _over(_Dasinh(tx, ty) - r, _Datan(tx, ty)) 

128 return r 

129 

130 def DParametric(self, Phi1, Phi2): 

131 '''I{Divided Difference} of the parametric wrt the geographic latitude. 

132 ''' 

133 fm1, e2m1 = self._fm1, self._e2m1 

134 tx, ty = Phi1.tan, Phi2.tan 

135 # DbetaDphi = Datan(fm1*tx, fm1*ty) * fm1 / Datan(tx, ty) 

136 # Datan(x, y) = 1 / (1 + x^2), if x == y 

137 # = (atan(y) - atan(x)) / (y-x), if x*y < 0 

138 # = atan((y-x) / (1 + x*y)) / (y-x), if x*y > 0 

139 txy = tx * ty 

140 if txy < 0 or (isinf(ty) and not tx): 

141 _a = atan 

142 r = _over(_a(fm1 * ty) - _a(fm1 * tx), _a(ty) - _a(tx)) 

143 elif tx == ty: # includes tx = ty = inf 

144 if txy > 1: # == tx**2 

145 txy = _1_over(txy) 

146 r = txy + e2m1 

147 else: 

148 r = txy * e2m1 + _1_0 

149 r = _over(fm1 * (txy + _1_0), r) 

150 else: 

151 if txy > 1: 

152 tx = _1_over(tx) 

153 ty = _1_over(ty) 

154 txy = tx * ty 

155 t = txy + e2m1 

156 else: 

157 t = txy * e2m1 + _1_0 

158 r = ty - tx 

159 r = _over(atan2(r * fm1, t), atan2(r, _1_0 + txy)) 

160 return r 

161 

162 def DRectifying(self, Phi1, Phi2): 

163 '''I{Divided Difference} of the rectifying wrt the geographic latitude. 

164 ''' 

165 # Stipulate that Phi1 and Phi2 are in [-90d, 90d] 

166 x, y = Phi1.toRadians, Phi2.toRadians 

167 if y == x: # isnear0 

168 Mu1 = self.Rectifying(Phi1, diff=True) 

169 tphi1, r = Phi1.tan, Mu1.diff 

170 if isfinite(tphi1): 

171 r *= _over(_sc(tphi1), _sc(Mu1.tan))**2 

172 else: # PYCHOK no cover 

173 r = _1_over(r) 

174 elif (x * y) < 0: 

175 r = _over(self.Rectifying(Phi2).toRadians - 

176 self.Rectifying(Phi1).toRadians, y - x) 

177 else: 

178 r = _over(self.b, self.RectifyingRadius(True)) 

179 r *= self.DE(*map1(self.Parametric, Phi1, Phi2)) 

180 r *= self.DParametric(Phi1, Phi2) 

181 return r # or INF or NAN 

182 

183 

184def _DClenshaw(sinp, Zeta1, Zeta2, cs, K): 

185 '''(INTERNAL) I{Divided Difference} of L{AuxLat._Clenshaw}. 

186 

187 @return: C{Fsum} if B{C{sinp}} otherwise a C{float}. 

188 ''' 

189 s1, c1, r1 = Zeta1._yxr_normalized(False) 

190 s2, c2, r2 = Zeta2._yxr_normalized(False) 

191 Delta = r2 - r1 

192 # Evaluate (Clenshaw(sinp, szeta2, czeta2, cs, K) - 

193 # Clenshaw(sinp, szeta1, czeta1, cs, K)) / Delta 

194 # or f = sin if sinp else cos 

195 # sum(cs[k] * (f((2*k+2) * Zeta2) - 

196 # f((2*k+2) * Zeta2))) / Delta 

197 # 

198 # Delta is EITHER 1, giving the plain difference OR (Zeta2 - Zeta1) 

199 # in radians, giving the I{Divided Difference}. Other values will 

200 # produce nonsense. 

201 # 

202 # Suffices a and b denote [1,1], [2,1] elements of matrix/vector 

203 cp = cm = c2 * c1 

204 t = s2 * s1 

205 cp -= t # not + 

206 cm += t # not - 

207 

208 sp = s2 * c1 

209 t = c2 * s1 

210 smd = ((sin(Delta) / Delta) if Delta != _1_0 else 

211 (sp - t)) if Delta else _1_0 

212 sp += t 

213 

214 xa = cp * cm * _2_0 

215 xb = sp * smd * _N_2_0 

216 xD = xb * Delta**2 

217 

218 if isfinite(xD) and isfinite(xb) and isfinite(xa): 

219 U0a, U1a = Fsum(), Fsum() 

220 U0b, U1b = Fsum(), Fsum() 

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 r = float(U0a) if sinp else U0a # Fsum 

241 else: 

242 r = _naninf(xD, xb, xa) 

243 return r 

244 

245 

246def _Dsin(x, y): # see also .rhumbx._Dsin 

247 r = cos((x + y) * _0_5) 

248 d = (x - y) * _0_5 

249 if d: 

250 r *= sin(d) / d 

251 return r 

252 

253 

254def _Dsn(x, y): 

255 # (sn(y) - sn(x)) / (y - x) 

256 if x != y: 

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

258 if (x * y) > 0: 

259 scx, scy = _sc(x), _sc(y) 

260 r = _over((snx / scy) + (sny / scx), 

261 (snx + sny) * scy * scx) 

262 else: 

263 r = (sny - snx) / (y - x) 

264 elif x: 

265 r = _1_over(_sc(x) * (x**2 + _1_0)) # == 1 / sqrt3(x**2 + 1) 

266 else: 

267 r = _1_0 

268 return r 

269 

270 

271def _sxk2y(sx, sy, k2): 

272 # .DE helper 

273 sy *= sy * k2 

274 if sy: 

275 try: 

276 sx *= sqrt(_1_0 - sy) 

277 except ValueError: # domain error 

278 sx = NAN 

279 return sx 

280 

281 

282__all__ += _ALL_DOCS(AuxDLat) 

283 

284# **) MIT License 

285# 

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

287# 

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

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

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

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

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

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

294# 

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

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

297# 

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

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

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

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

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

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

304# OTHER DEALINGS IN THE SOFTWARE.