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

2

3 u'''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>}

5 in I{GeographicLib version 2.2+}.

6

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

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

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

13

14 from pygeodesy.auxilats.auxily import Aux, _Datan, _Dasinh, _Dm, _sc, _sn, \

15 atan1, AuxError

16 from pygeodesy.auxilats.auxLat import AuxLat, _ALL_DOCS

17 from pygeodesy.basics import map1, _reverange

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

19 _1_0, _2_0, _N_2_0, _naninf, _over, _1_over

20 from 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# from pygeodesy.utily import atan1 # from .auxilats.auxily

25

26 from math import atan2, cos, sin, sqrt

27

28__all__ = ()

29__version__ = '23.09.26'

30

31

32 class AuxDLat(AuxLat):

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

34 latitudes and other C{Divided Differences} needed for

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

36 '''

37

38 def CParametric(self, Zeta1, Zeta2):

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

40 '''

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

42

43 def CRectifying(self, Zeta1, Zeta2):

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

45 '''

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

47

48 def _Datanhee(self, x, y):

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

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

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

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

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

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

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

56 if self.f < 0:

57 e = self._e

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

59 else:

60 x *= self._fm1

61 y *= self._fm1

62 e1 = self._e1

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

64 return _Dsn(x, y) * r

65

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

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

68 '''

69 auxin = Zeta1._AUX

70 # assert Zeta2._AUX == auxin

71 try:

72 if auxout != auxin:

73 cs = self._coeffs(auxout, auxin)

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

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

76 else:

77 r = _1_0

78 except AuxError: # no _coeffs

79 r = NAN

80 return r

81

82 def DE(self, X, Y):

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

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

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

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

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

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

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

90 # also stipulate that x != y.

91

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

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

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

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

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

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

98 sx, cx = cx, sx

99 sy, cy = cy, sy

100 d, k2 = -d, self._e2

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

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

103 if Dt:

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

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

106 t = d * Dt

107 t2 = _1_0 + t**2

108 Dt *= _2_0 / t2

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

110 d2 = _1_0 - sk2

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

112 # E(z)/sin(z)

113 Dt *= _Ef._RFRD(c2, d2, _1_0, sk2) - k2 * sx * sy

114 return Dt

115

116 def DIsometric(self, Phi1, Phi2):

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

118 '''

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

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

121 r = NAN

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

123 r = INF

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

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

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

127 return r

128

129 def DParametric(self, Phi1, Phi2):

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

131 '''

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

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

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

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

136 # = (atan1(y) - atan1(x)) / (y-x) if x*y < 0

137 # = atan1(y-x, x*y + 1) / (y-x) if x*y > 0

138 txy = tx * ty

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

140 _a = atan1

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

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

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

144 txy = _1_over(txy)

145 r = txy + e2m1

146 else:

147 r = txy * e2m1 + _1_0

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

149 else:

150 if txy > 1:

151 tx = _1_over(tx)

152 ty = _1_over(ty)

153 txy = _1_over(txy)

154 t = txy + e2m1

155 else:

156 t = txy * e2m1 + _1_0

157 r = ty - tx

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

159 return r

160

161 def DRectifying(self, Phi1, Phi2):

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

163 '''

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

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

166 if y == x: # isnear0

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

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

169 if isfinite(tphi1):

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

171 else: # PYCHOK no cover

172 r = _1_over(r)

173 elif (x * y) < 0:

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

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

176 else:

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

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

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

180 return r # or INF or NAN

181

182

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

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

185

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

187 '''

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

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

190 Delta = r2 - r1

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

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

193 # or f = sin if sinp else cos

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

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

196 #

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

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

199 # produce nonsense.

200 #

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

202 cp = cm = c2 * c1

203 t = s2 * s1

204 cp -= t # not +

205 cm += t # not -

206

207 sp = s2 * c1

208 t = c2 * s1

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

210 (sp - t)) if Delta else _1_0

211 sp += t

212

213 xa = cp * cm * _2_0

214 xb = sp * smd * _N_2_0

215 xD = xb * Delta**2

216

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

218 U0a, U1a = Fsum(), Fsum()

219 U0b, U1b = Fsum(), Fsum()

220 for k in _reverange(K): # assert len(cs) == K

221 # t = x . U0 - U1 + cs[k] * I

222 U1a -= U0a * xa + U0b * xD + cs[k]

223 U1b -= U0a * xb + U0b * xa

224 U1a, U0a = U0a, -U1a

225 U1b, U0b = U0b, -U1b

226 # F0a = (sp if sinp else cp) * cm

227 # F0b = (cp if sinp else -sp) * smd

228 # Fm1a = 0 if sinp else 1 # Fm1b = 0

229 # return (U0b * F0a + U0a * F0b - U1b * Fm1a) * 2

230 if sinp:

231 U1b = _0_0

232 else:

233 sp, cp = cp, -sp

234 U0b *= sp * cm

235 U0a *= cp * smd

236 U0a += U0b

237 U0a = _Dm(U0a, U1b, _2_0)

238 r = float(U0a) if sinp else U0a # Fsum

239 else:

240 r = _naninf(xD, xb, xa)

241 return r

242

243

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

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

246 d = (x - y) * _0_5

247 if d:

248 r *= sin(d) / d

249 return r

250

251

252 def _Dsn(x, y):

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

254 if x != y:

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

256 if (x * y) > 0:

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

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

259 (snx + sny) * scy * scx)

260 else:

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

262 elif x:

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

264 else:

265 r = _1_0

266 return r

267

268

269 def _sxk2y(sx, sy, k2):

270 # .DE helper

271 sy *= sy * k2

272 if sy:

273 try:

274 sx *= sqrt(_1_0 - sy)

275 except ValueError: # domain error

276 sx = NAN

277 return sx

278

279

280__all__ += _ALL_DOCS(AuxDLat)

281

282# **) MIT License

283#

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

285#

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

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

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

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

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

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

292#

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

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

295#

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

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

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

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

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

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

302# OTHER DEALINGS IN THE SOFTWARE.