Coverage for pygeodesy/auxilats/auxLat.py: 93%

410 statements  

« prev     ^ index     » next       coverage.py v7.2.2, created at 2023-08-28 15:52 -0400

1 

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

3 

4u'''Class L{AuxLat} transcoded to Python from I{Karney}'s C++ class U{AuxLatitude 

5<https://GeographicLib.SourceForge.io/C++/doc/classGeographicLib_1_1AuxLatitude.html>} 

6in I{GeographicLib version 2.2+}. 

7 

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

9under the MIT/X11 License. For more information, see the U{GeographicLib 

10<https://GeographicLib.SourceForge.io>} documentation. 

11 

12@see: U{Auxiliary latitudes<https:#GeographicLib.SourceForge.io/C++/doc/auxlat.html>} 

13 U{On auxiliary latitudes<https:#ArXiv.org/abs/2212.05818>}. 

14''' 

15# make sure int/int division yields float quotient, see .basics 

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

17 

18from pygeodesy.auxilats.auxAngle import AuxAngle, AuxBeta, AuxChi, _AuxClass, \ 

19 AuxMu, AuxPhi, AuxTheta, AuxXi 

20from pygeodesy.auxilats.auxily import Aux, _sc, _sn, _Ufloats 

21from pygeodesy.basics import isscalar, _reverange, _xinstanceof 

22from pygeodesy.constants import INF, MAX_EXP, MIN_EXP, NAN, PI_2, PI_4, _EPSqrt, \ 

23 _0_0, _0_0s, _0_1, _0_25, _0_5, _1_0, _2_0, _3_0, \ 

24 _360_0, isfinite, isinf, isnan, _log2, _over 

25from pygeodesy.datums import _ellipsoidal_datum, _WGS84, Ellipsoid 

26# from pygeodesy.ellipsoids import Ellipsoid # from .datums 

27from pygeodesy.elliptic import Elliptic as _Ef 

28from pygeodesy.errors import AuxError, _xkwds, _xkwds_get, _Xorder 

29# from pygeodesy.fmath import cbrt # from .karney 

30from pygeodesy.fsums import Fsum, _sum 

31from pygeodesy.karney import _2cos2x, _ALL_DOCS, cbrt, _MODS 

32from pygeodesy.interns import NN, _DOT_, _UNDER_ # _earth_ 

33# from pygeodesy.lazily import _ALL_DOCS, _ALL_MODS as _MODS # from .karney 

34from pygeodesy.props import Property, Property_RO, _update_all 

35from pygeodesy.units import Degrees, Meter 

36# from pygeodesy.utily import _passarg # _MODS 

37 

38from math import asinh, atan, atan2, copysign, cosh, fabs, sin, sinh, sqrt 

39try: 

40 from math import exp2 as _exp2 

41except ImportError: # Python 3.11- 

42 

43 def _exp2(x): 

44 return pow(_2_0, x) 

45 

46__all__ = () 

47__version__ = '23.08.27' 

48 

49_TRIPS = 1024 # XXX 2 or 3? 

50 

51 

52class AuxLat(AuxAngle): 

53 '''Accurate conversion between I{Auxiliary} latitudes on an ellipsoid. 

54 

55 Latitudes are represented by instances of L{AuxAngle} in order to 

56 maintain precision close to the poles, I{Authalic} latitude C{Xi}, 

57 I{Conformal} C{Chi}, I{Geocentric} C{Theta}, I{Geographic} C{Phi}, 

58 I{Parametric} C{Beta} and I{Rectifying} C{Mu}. 

59 

60 @see: I{Karney}'s C++ class U{AuxLatitude 

61 <https://GeographicLib.SourceForge.io/C++/doc/classGeographicLib_1_1AuxLatitude.html>}. 

62 ''' 

63 _csc = dict() # global coeffs cache: [aL][k], upto max(k) * (4 + 6 + 8) floats 

64 _E = _WGS84.ellipsoid 

65 _mAL = 6 # 4, 6 or 8 aka Lmax 

66 

67 def __init__(self, a_earth=_WGS84, f=None, b=None, name=NN, **ALorder): 

68 '''New L{AuxLat} instance on an ellipsoid or datum. 

69 

70 @arg a_earth: Equatorial radius, semi-axis (C{meter}) or an 

71 ellipsoid or datum (L{Datum}, L{Ellipsoid}, 

72 L{Ellipsoid2} or L{a_f2Tuple}). 

73 @kwarg f: Flattening: M{(a - b) / a} (C{float}, near zero for 

74 spherical), ignored if B{C{a_earth}} is not scalar. 

75 @kwarg b: Optional polar radius, semi-axis (C{meter}, same 

76 units as B{C{a_earth}}), ignored if B{C{a_earth}} 

77 is not scalar. 

78 @kwarg ALorder: Optional keyword arguments B{C{ALorder}} to 

79 set L{AuxLat}'s C{order}, see property 

80 C{ALorder}. 

81 @kwarg name: Optional name (C{str}). 

82 ''' 

83 if a_earth is not _WGS84: 

84 n = name or AuxLat.__name__ 

85 try: 

86 if b is f is None: 

87 E = _ellipsoidal_datum(a_earth, name=n).ellipsoid # XXX raiser=_earth_ 

88 elif isscalar(a_earth): 

89 E = Ellipsoid(a_earth, f=f, b=b, name=_UNDER_(n)) 

90 else: 

91 raise ValueError() 

92 except Exception as x: 

93 raise AuxError(a_earth=a_earth, f=f, b=b, cause=x) 

94 self._E = E 

95 

96 if name: 

97 self.name = name 

98 if ALorder: 

99 self.ALorder = _xkwds_get(ALorder, ALorder=AuxLat._mAL) 

100 

101 @Property_RO 

102 def a(self): 

103 '''Get the equatorial radius (C{Meter}, conventionally). 

104 ''' 

105 return self.ellipsoid.a 

106 

107 equatoradius = a 

108 

109 @Property 

110 def ALorder(self): 

111 '''Get the I{AuxLat} order (C{int}, 4, 6 or 8). 

112 ''' 

113 return self._mAL 

114 

115 @ALorder.setter # PYCHOK setter! 

116 def ALorder(self, order): 

117 '''Set the I{AuxLat} order (C{int}, 4, 6 or 8). 

118 ''' 

119 m = _Xorder(_AR2Coeffs, AuxError, ALorder=order) 

120 if self._mAL != m: 

121 _update_all(self) 

122 self._mAL = m 

123 

124 def _atanhee(self, tphi): # see Ellipsoid._es_atanh, .albers._atanhee 

125 # atanh(e * sphi) = asinh(e' * sbeta) 

126 f = self.f 

127 s = _sn(self._fm1 * tphi) if f > 0 else _sn(tphi) 

128 if f: # atanh(e * sphi) = asinh(e' * sbeta) 

129 e = self._e 

130 s = (atan(e * s) if f < 0 else asinh(self._e1 * s)) / e 

131 return s 

132 

133 def Authalic(self, Phi, **diff_name): 

134 '''Convert I{Geographic} to I{Aunthalic} latitude. 

135 

136 @arg Phi: Geographic latitude (L{AuxAngle}). 

137 

138 @return: Parametric latitude, C{Beta} (L{AuxAngle}). 

139 ''' 

140 _xinstanceof(AuxAngle, Phi=Phi) 

141 # assert Phi._AUX == Aux.PHI 

142 tphi = fabs(Phi.tan) 

143 if isfinite(tphi) and tphi and self.f: 

144 y, x = Phi._yx_normalized 

145 q = self._q 

146 qv = self._qf(tphi) 

147 Dq2 = self._Dq(tphi) 

148 Dq2 *= (q + qv) / (fabs(y) + _1_0) # _Dq(-tphi) 

149 Xi = AuxXi(copysign(qv, Phi.y), x * sqrt(Dq2), 

150 name=_xkwds_get(diff_name, name=Phi.name)) 

151 

152 if _xkwds_get(diff_name, diff=False): 

153 if isnan(tphi): 

154 d = self._e2m1_sq2 

155 else: 

156 c = self.Parametric(Phi)._x_normalized 

157 d = _over(c, Xi._x_normalized)**3 

158 d *= _over(c, x) * _over(_2_0, q) 

159 Xi._diff = d 

160 else: 

161 Xi = AuxXi(*Phi._yx) # diff default 

162 # assert Xi._AUX == Aux.XI 

163 return Xi 

164 

165 def AuthalicRadius2(self, exact=False, f_max=_0_1): 

166 '''Get the I{Authalic} radius I{squared}. 

167 

168 @kwarg exact: If C{True}, use the exact expression, otherwise 

169 the I{Taylor} series. 

170 @kwarg f_max: C{Flattening} not to exceed (C{float}). 

171 

172 @return: Authalic radius I{squared} (C{meter} I{squared}, same 

173 units as the ellipsoid axes). 

174 

175 @raise AuxError: If C{B{exact}=False} and C{abs(flattening)} 

176 exceeds C{f_max}. 

177 ''' 

178 f = self.f 

179 if exact or not f: 

180 c2 = self.ellipsoid.b2 * self._q # == ellipsoid.c2x * 2 

181 elif fabs(f) < f_max: 

182 # Using a * (a + b) / 2 as the multiplying factor leads to a rapidly 

183 # converging series in n. Of course, using this series isn't really 

184 # necessary, since the exact expression is simple to evaluate. However, 

185 # we do it for consistency with RectifyingRadius and, presumably, the 

186 # roundoff error is smaller compared to that for the exact expression. 

187 m = self.ALorder 

188 c2 = _MODS.karney._polynomial(self._n, _AR2Coeffs[m], 0, m) 

189 c2 *= self.a * (self.a + self.b) 

190 else: 

191 raise AuxError(exact=exact, f=f, f_max=f_max) 

192 return c2 * _0_5 

193 

194 @Property_RO 

195 def b(self): 

196 '''Get the polar radius (C{Meter}, conventionally). 

197 ''' 

198 return self.ellipsoid.b # a * (_1_0 - f) 

199 

200 polaradius = b 

201 

202 def _coeffs(self, auxout, auxin): 

203 # Get the polynomial coefficients as 4-, 6- or 8-tuple 

204 aL = self.ALorder # aka Lmax 

205 k = Aux._1d(auxout, auxin) 

206 try: 

207 cs = AuxLat._csc[aL][k] 

208 except KeyError: 

209 if auxout == auxin: 

210 cs = _0_0s(aL) # not cached 

211 else: 

212 Cx = _CXcoeffs(aL) 

213 try: 

214 Cx = Cx[auxout][auxin] 

215 except KeyError as x: 

216 raise AuxError(auxout=auxout, auxin=auxin, cause=x) 

217 

218 d = x = n = self._n 

219 if Aux.use_n2(auxin) and Aux.use_n2(auxout): 

220 x = self._n2 

221 

222 def _m(m): 

223 return m // 2 

224 else: 

225 def _m(m): # PYCHOK expected 

226 return m 

227 

228 i = 0 

229 cs = [] 

230 _p = _MODS.karney._polynomial 

231 for r in _reverange(aL): 

232 j = i + _m(r) + 1 # order m = j - i - 1 

233 cs.append(_p(x, Cx, i, j) * d) 

234 d *= n 

235 i = j 

236 # assert i == len(Cx) and len(cs) == aL 

237 AuxLat._csc.setdefault(aL, {})[k] = tuple(cs) 

238 return cs 

239 

240 def Conformal(self, Phi, **diff_name): 

241 '''Convert I{Geographic} to I{Conformal} latitude. 

242 

243 @arg Phi: Geographic latitude (L{AuxAngle}). 

244 

245 @return: Conformal latitude, C{Chi} (L{AuxAngle}). 

246 ''' 

247 _xinstanceof(AuxAngle, Phi=Phi) 

248 # assert Phi._AUX == Aux.PHI 

249 tphi = tchi = fabs(Phi.tan) 

250 if isfinite(tphi) and tphi and self.f: 

251 sig = sinh(self._e2 * self._atanhee(tphi)) 

252 scsig = _sc(sig) 

253 scphi = _sc(tphi) 

254 if self.f > 0: 

255 # The general expression for tchi is 

256 # tphi * scsig - sig * scphi 

257 # This involves cancellation if f > 0, so change to 

258 # (tphi - sig) * (tphi + sig) / (tphi * scsig + sig * scphi) 

259 # To control overflow, write as (sigtphi = sig / tphi) 

260 # (tphi - sig) * (1 + sigtphi) / (scsig + sigtphi * scphi) 

261 sigtphi = sig / tphi 

262 if sig < (tphi * _0_5): 

263 t = tphi - sig 

264 else: 

265 def _asinh_2(x): 

266 return asinh(x) * _0_5 

267 # Still possibly dangerous cancellation in tphi - sig. 

268 # Write tphi - sig = (1 - e) * Dg(1, e) 

269 # Dg(x, y) = (g(x) - g(y)) / (x - y) 

270 # g(x) = sinh(x * atanh(sphi * x)) 

271 # Note sinh(atanh(sphi)) = tphi 

272 # Turn the crank on divided differences, substitute 

273 # sphi = tphi / _sc(tphi) 

274 # atanh(x) = asinh(x / sqrt(1 - x^2)) 

275 e = self._e 

276 em1 = self._e2m1 / (_1_0 + e) 

277 # assert em1 != 0 

278 scb = self._scbeta(tphi) 

279 scphib = scphi / scb # sec(phi) / sec(beta) 

280 atphib = _asinh_2(tphi * e / scb) # atanh(e * sphi) 

281 atphi = _asinh_2(tphi) # atanh(sphi) 

282 t = _asinh_2(em1 * (tphi * scphib)) / em1 

283 try: 

284 Dg = Fsum(atphi, atphib, t, e * t) 

285 except ValueError: # Fsum(NAN) exception 

286 Dg = _sum((atphi, atphib, t, e * t)) 

287 e *= atphib 

288 t = atphi - e 

289 if t: # sinh(0) == 0 

290 Dg *= sinh(t) / t * cosh(atphi + e) * em1 

291 t = float(Dg) # tphi - sig 

292 tchi = _over(t * (_1_0 + sigtphi), 

293 scsig + scphi * sigtphi) # if t else _0_0 

294 else: 

295 tchi = tphi * scsig - sig * scphi 

296 

297 n = _xkwds_get(diff_name, name=Phi.name) 

298 Chi = AuxChi(tchi, name=n).copyquadrant(Phi) 

299 

300 if _xkwds_get(diff_name, diff=False): 

301 if isinf(tphi): # PYCHOK np cover 

302 d = self._conformal_diff 

303 else: 

304 d = self.Parametric(Phi)._x_normalized 

305 if d: 

306 d = _over(d, Chi._x_normalized) * \ 

307 _over(d, Phi._x_normalized) * self._e2m1 

308 Chi._diff = d 

309 # assrt Chi._AUX == Aux.CHI 

310 return Chi 

311 

312 @Property_RO 

313 def _conformal_diff(self): # PYCHOK no cover 

314 '''(INTERNAL) Constant I{Conformal} diff. 

315 ''' 

316 e = self._e 

317 if self.f > 0: 

318 ss = sinh(asinh(self._e1) * e) 

319 d = _over(_1_0, _sc(ss) + ss) 

320 elif e: 

321 ss = sinh(-atan(e) * e) 

322 d = _sc(ss) - ss 

323 else: 

324 d = _1_0 

325 return d 

326 

327 def convert(self, auxout, Zeta_d, exact=False): 

328 # Convert I{Auxiliary} or I{scalar} latitude 

329 Z = d = Zeta_d 

330 if isinstance(Z, AuxAngle): 

331 A, auxin = _AuxClass(auxout), Z._AUX 

332 if auxin == auxout: 

333 pass 

334 elif not (isfinite(Z.tan) and Z.tan): # XXX 

335 Z = A(*Z._yx, aux=auxout, name=Z.name) 

336 elif exact: 

337 p = Aux.power(auxout, auxin) 

338 if p is None: 

339 P = self._fromAux(Z) # Phi 

340 Z = self._toAux(auxout, P) 

341 Z._iter = P.iteration 

342 else: 

343 y, x = Z._yx 

344 if p: 

345 y *= pow(self._fm1, p) 

346 Z = A(y, x, aux=auxout, name=Z.name) 

347 else: 

348 cs = self._coeffs(auxout, auxin) 

349 yx = Z._yx_normalized 

350 Z = A(*yx, aux=auxout, name=Z.name) 

351 # assert Z._yx == yx 

352 r = _Clenshaw(True, Z, cs, self.ALorder) 

353 Z += AuxAngle.fromRadians(r, aux=auxout) 

354 # assert Z._AUX == auxout 

355 return Z 

356 

357 elif isinstance(d, Degrees) or isscalar(d): 

358 Z = AuxPhi.fromDegrees(d) 

359 d = round((d - Z.toDegrees) / _360_0) * _360_0 

360 d += self.convert(auxout, Z, exact).toDegrees 

361 return Degrees(d, name=Aux.Greek(auxout)) 

362 

363 raise AuxError(auxout=auxout, Zeta_d=Zeta_d, exact=exact) 

364 

365 def _Dq(self, tphi): 

366 # I{Divided Difference} of (q(1) - q(sphi)) / (1 - sphi). 

367 sphi = _sn(tphi) 

368 if tphi > 0: 

369 scphi = _sc(tphi) 

370 d = _1_0 / (scphi**2 * (_1_0 + sphi)) # XXX - sphi 

371 if d: 

372 # General expression for _Dq(1, sphi) is 

373 # atanh(e * d / (1 - e2 * sphi)) / (e * d) + 

374 # (1 + e2 * sphi) / ((1 - e2 * sphi * sphi) * e2m1) 

375 # with atanh(e * d / (1 - e2 * sphi)) = 

376 # atanh(e * d * scphi / (scphi - e2 * tphi)) 

377 e2m1, ed = self._e2m1, (self._e * d) 

378 if e2m1 and ed: 

379 e2 = self._e2 

380 if e2 > 0: # assert self.f > 0 

381 scb = self._scbeta(tphi) 

382 q = scphib = scphi / scb 

383 q *= (scphi + tphi * e2) / (e2m1 * scb) 

384 q += asinh(self._e1 * d * scphib) / ed 

385 else: 

386 s2 = sphi * e2 

387 q = (_1_0 + s2) / ((_1_0 - sphi * s2) * e2m1) 

388 q += (atan2(ed, _1_0 - s2) / ed) if e2 < 0 else _1_0 

389 else: # PYCHOK no cover 

390 q = INF 

391 else: # PYCHOK no cover 

392 q = self._2_e2m12 

393 else: # not reached, open-coded in .Authalic 

394 q = _over(self._q - self._qf(tphi), _1_0 - sphi) 

395 return q 

396 

397 @Property_RO 

398 def _e(self): # unsigned, (1st) eccentricity 

399 return self.ellipsoid.e # sqrt(fabs(self._e2)) 

400 

401 @Property_RO 

402 def _e1(self): 

403 return sqrt(fabs(self._e12)) 

404 

405 @Property_RO 

406 def _e12(self): 

407 return _over(self._e2, _1_0 - self._e2) 

408 

409 @Property_RO 

410 def _e12p1(self): 

411 return _1_0 / self._e2m1 

412 

413 @Property_RO 

414 def _e2(self): # signed, (1st) eccentricity squared 

415 return self.ellipsoid.e2 

416 

417 @Property_RO 

418 def _e2m1(self): # 1 less 1st eccentricity squared 

419 return self.ellipsoid.e21 # == ._fm1**2 

420 

421 @Property_RO 

422 def _e2m1_sq2(self): 

423 return self._e2m1 * sqrt(self._q * _0_5) 

424 

425 @Property_RO 

426 def _2_e2m12(self): 

427 return _2_0 / self._e2m1**2 

428 

429 @Property_RO 

430 def _Ef_fRG_a2b2_PI_4(self): 

431 E = self.ellipsoid 

432 return _Ef.fRG(E.a2, E.b2) / PI_4 

433 

434 @Property_RO 

435 def ellipsoid(self): 

436 '''Get the ellipsoid (L{Ellipsoid}). 

437 ''' 

438 return self._E 

439 

440 @Property_RO 

441 def f(self): 

442 '''Get the flattening (C{float}). 

443 ''' 

444 return self.ellipsoid.f 

445 

446 flattening = f 

447 

448 @Property_RO 

449 def _fm1(self): # 1 - flattening 

450 return self.ellipsoid.f1 

451 

452 def _fromAux(self, Zeta, **name): 

453 '''Convert I{Auxiliary} to I{Geographic} latitude. 

454 

455 @arg Zeta: Auxiliary latitude (L{AuxAngle}). 

456 

457 @return: Geographic latitude, I{Phi} (L{AuxAngle}). 

458 ''' 

459 _xinstanceof(AuxAngle, Zeta=Zeta) 

460 aux = Zeta._AUX 

461 n = _xkwds_get(name, name=Zeta.name) 

462 f = self._fromAuxCase.get(aux, None) 

463 if f is None: 

464 Phi = AuxPhi(NAN, name=n) 

465 elif callable(f): 

466 t = self._fm1 

467 t *= f(t) 

468 Phi = _Newton(t, Zeta, self._toZeta(aux), name=n) 

469 else: # assert isscalar(f) 

470 y, x = Zeta._yx 

471 Phi = AuxPhi(y / f, x, name=n) 

472 # assert Phi._AUX == Aux.PHI 

473 return Phi 

474 

475 @Property_RO 

476 def _fromAuxCase(self): 

477 '''(INTERNAL) switch(auxin): ... 

478 ''' 

479 return {Aux.AUTHALIC: cbrt, 

480 Aux.CONFORMAL: _MODS.utily._passarg, 

481 Aux.GEOCENTRIC: self._e2m1, 

482 Aux.GEOGRAPHIC: _1_0, 

483 Aux.PARAMETRIC: self._fm1, 

484 Aux.RECTIFYING: sqrt} 

485 

486 def Geocentric(self, Phi, **diff_name): 

487 '''Convert I{Geographic} to I{Geocentric} latitude. 

488 

489 @arg Phi: Geographic latitude (L{AuxAngle}). 

490 

491 @return: Geocentric latitude, C{Phi} (L{AuxAngle}). 

492 ''' 

493 _xinstanceof(AuxAngle, Phi=Phi) 

494 # assert Phi._AUX == Aux.PHI 

495 Theta = AuxTheta(Phi.y * self._e2m1, Phi.x, 

496 name=_xkwds_get(diff_name, name=Phi.name)) 

497 if _xkwds_get(diff_name, diff=False): 

498 Theta._diff = self._e2m1 

499 return Theta 

500 

501 def Geodetic(self, Phi, **diff_name): # PYCHOK no cover 

502 '''Convert I{Geographic} to I{Geodetic} latitude. 

503 

504 @arg Phi: Geographic latitude (L{AuxAngle}). 

505 

506 @return: Geodetic latitude, C{Phi} (L{AuxAngle}). 

507 ''' 

508 _xinstanceof(AuxAngle, Phi=Phi) 

509 # assert Phi._AUX == Aux.PHI 

510 return AuxPhi(Phi, name=_xkwds_get(diff_name, name=Phi.name)) 

511 

512 @Property_RO 

513 def _n(self): # 3rd flattening 

514 return self.ellipsoid.n 

515 

516 @Property_RO 

517 def _n2(self): 

518 return self._n**2 

519 

520 def Parametric(self, Phi, **diff_name): 

521 '''Convert I{Geographic} to I{Parametric} latitude. 

522 

523 @arg Phi: Geographic latitude (L{AuxAngle}). 

524 

525 @return: Parametric latitude, C{Beta} (L{AuxAngle}). 

526 ''' 

527 _xinstanceof(AuxAngle, Phi=Phi) 

528 # assert Phi._AUX == Aux.PHI 

529 Beta = AuxBeta(Phi.y * self._fm1, Phi.x, 

530 name=_xkwds_get(diff_name, name=Phi.name)) 

531 if _xkwds_get(diff_name, diff=False): 

532 Beta._diff = self._fm1 

533 return Beta 

534 

535 Reduced = Parametric 

536 

537 @Property_RO 

538 def _q(self): # constant _q 

539 q, f = self._e12p1, self.f 

540 if f: 

541 e = self._e 

542 q += (asinh(self._e1) if f > 0 else atan(e)) / e 

543 else: 

544 q += _1_0 

545 return q 

546 

547 def _qf(self, tphi): 

548 # function _q: atanh(e * sphi) / e + sphi / (1 - (e * sphi)^2) 

549 scb = self._scbeta(tphi) 

550 return self._atanhee(tphi) + (tphi / scb) * (_sc(tphi) / scb) 

551 

552 def _qIntegrand(self, beta): 

553 # pbeta(beta) = integrate(q(beta), beta), with beta in radians 

554 # q(beta) = (1-f) * (sin(xi) - sin(chi)) / cos(phi) 

555 # = (1-f) * (cos(chi) - cos(xi)) / cos(phi) * 

556 # (cos(xi) + cos(chi)) / (sin(xi) + sin(chi)) 

557 # Fit q(beta)/cos(beta) with Fourier transform 

558 # q(beta)/cos(beta) = sum(c[k] * sin((2*k+1)*beta), k, 0, K-1) 

559 # then the integral is 

560 # pbeta = sum(d[k] * cos((2*k+2)*beta), k, 0, K-1) 

561 # where 

562 # d[k] = -1/(4*(k+1)) * (c[k] + c[k+1]) for k in 0..K-2 

563 # d[K-1] = -1/(4*K) * c[K-1] 

564 Beta = AuxBeta.fromRadians(beta) 

565 if Beta.x: # and self._fm1: 

566 Ax, _cv = Aux, self.convert 

567 Phi = _cv(Ax.PHI, Beta, exact=True) 

568 schi, cchi = _cv(Ax.CHI, Phi, exact=True)._yx_normalized 

569 sxi, cxi = _cv(Ax.XI, Phi, exact=True)._yx_normalized 

570 r = (sxi - schi) if fabs(schi) < fabs(cchi) else \ 

571 _over(_2cos2x(cchi, cxi), (sxi + schi) * _2_0) 

572 r *= _over(self._fm1, Phi._x_normalized * Beta._x_normalized) 

573 else: # beta == PI_2, PI3_2, ... 

574 r = _0_0 # XXX 0 avoids NAN summation exceptions 

575 return r 

576 

577 def Rectifying(self, Phi, **diff_name): 

578 '''Convert I{Geographic} to I{Rectifying} latitude. 

579 

580 @arg Phi: Geographic latitude (L{AuxAngle}). 

581 

582 @return: Rectifying latitude, C{Mu} (L{AuxAngle}). 

583 ''' 

584 Beta = self.Parametric(Phi) 

585 # assert Beta._AUX == Aux.BETA 

586 sb, cb = map(fabs, Beta._yx_normalized) 

587 a, ka, ka1 = _1_0, self._e2, self._e2m1 

588 b, kb, kb1 = self._fm1, -self._e12, self._e12p1 

589 if self.f < 0: 

590 a, b = b, a 

591 ka, kb = kb, ka 

592 ka1, kb1 = kb1, ka1 

593 sb, cb = cb, sb 

594 # now a, b = larger/smaller semiaxis 

595 # Beta measured from larger semiaxis 

596 # kb, ka = modulus-squared for distance from Beta = 0, pi/2 

597 # NB kb <= 0; 0 <= ka <= 1 

598 # sa = b*E(Beta, sqrt(kb)) 

599 # sb = a*E(Beta',sqrt(ka)) 

600 # 1 - ka * (1 - sb2) = 1 - ka + ka*sb2 

601 sb2 = sb**2 

602 cb2 = cb**2 

603 da2 = ka1 + ka * sb2 

604 db2 = _1_0 - kb * sb2 

605 # DLMF Eq. 19.25.9 

606 my = b * sb * (_Ef.fRF(cb2, db2, _1_0) - 

607 kb * sb2 * _Ef.fRD(cb2, db2, _1_0, _3_0)) 

608 # DLMF Eq. 19.25.10 with complementary angles 

609 mx = a * cb * (_Ef.fRF(sb2, da2, _1_0) * ka1 + 

610 ka * cb2 * _Ef.fRD(sb2, _1_0, da2, _3_0) * ka1 + 

611 ka * sb / sqrt(da2)) 

612 # my + mx = 2*_Ef.fRG(a*a, b*b) = a*E(e) = b*E(i*e') 

613 # mr = a*E(e)*(2/pi) = b*E(i*e')*(2/pi) 

614 if self.f < 0: 

615 a, b = b, a 

616 my, mx = mx, my 

617 mr = (my + mx) / PI_2 

618 if mr: 

619 my = sin(my / mr) 

620 mx = sin(mx / mr) # XXX zero? 

621 else: # zero Mu 

622 my, mx = _0_0, _1_0 

623 n = _xkwds_get(diff_name, name=Phi.name) 

624 Mu = AuxMu(my, mx, # normalized 

625 name=n).copyquadrant(Phi) 

626 

627 if _xkwds_get(diff_name, diff=False): 

628 d, x = _0_0, Beta._x_normalized 

629 if x and mr: 

630 if Mu.x and Phi.x and not isinf(Phi.tan): 

631 d = b / mr * (x / Mu.x)**2 \ 

632 * (x / Phi._x_normalized) 

633 else: 

634 d = mr / a 

635 Mu._diff = self._fm1 * d 

636 return Mu 

637 

638 def RectifyingRadius(self, exact=False): 

639 '''Get the I{Rectifying} radius. 

640 

641 @arg exact: If C{True}, use the exact expression, 

642 otherwise the I{Taylor} series. 

643 

644 @return: Rectifying radius (L{Meter}, same units 

645 as the ellipsoid axes). 

646 ''' 

647 r = self._Ef_fRG_a2b2_PI_4 if exact else self._RectifyingR 

648 return Meter(r, name=self.RectifyingRadius.__name__) 

649 

650 @Property_RO 

651 def _RectifyingR(self): 

652 m = self.ALorder 

653 d = _MODS.karney._polynomial(self._n2, _RRCoeffs[m], 0, m // 2) 

654 return d * (self.a + self.b) * _0_5 

655 

656 def _scbeta(self, tphi): 

657 return _sc(self._fm1 * tphi) 

658 

659 def _toAux(self, auxout, Phi, **diff_name): 

660 '''Convert I{Geographic} to I{Auxiliary} latitude. 

661 

662 @arg auxout: I{Auxiliary} kind (C{Aux.KIND}). 

663 @arg Phi: Geographic latitude (L{AuxLat}). 

664 

665 @return: Auxiliary latitude, I{Eta} (L{AuxLat}). 

666 ''' 

667 _xinstanceof(AuxAngle, Phi=Phi) 

668 # assert Phi._AUX == Aux.PHI 

669 n = _xkwds_get(diff_name, name=Phi.name) 

670 m = _toAuxCase.get(auxout, None) 

671 if m: # callable 

672 A = m(self, Phi, **_xkwds(diff_name, name=n)) 

673 elif auxout == Aux.GEODETIC: # == GEOGRAPHIC 

674 A = AuxPhi(Phi, name=n) 

675 else: # auxout? 

676 A = AuxPhi(NAN, name=n) 

677 # assert A._AUX == auxout 

678 return A 

679 

680 def _toZeta(self, zetaux): 

681 '''Return a (lean) function to create C{AuxPhi(tphi)} and 

682 convert that into C{AuxAngle} of (fixed) kind C{zetaux} 

683 for use only inside the C{_Newton} loop. 

684 ''' 

685 class _AuxPhy(AuxPhi): 

686 # lean C{AuxPhi} instance. 

687 # _diff = _1_0 

688 # _x = _1_0 

689 

690 def __init__(self, tphi): # PYCHOK signature 

691 self._y = tphi 

692 

693 m = _toAuxCase.get(zetaux, None) 

694 if m: # callable 

695 

696 def _toZeta(tphi): 

697 return m(self, _AuxPhy(tphi), diff=True) 

698 

699 elif zetaux == Aux.GEODETIC: # GEOGRAPHIC 

700 _toZeta = _AuxPhy 

701 

702 else: # zetaux? 

703 

704 def _toZeta(unused): # PYCHOK expected 

705 return _AuxPhy(NAN) 

706 

707 return _toZeta 

708 

709 

710# switch(auxout): ... 

711_toAuxCase = {Aux.AUTHALIC: AuxLat.Authalic, 

712 Aux.CONFORMAL: AuxLat.Conformal, 

713 Aux.GEOCENTRIC: AuxLat.Geocentric, 

714 Aux.PARAMETRIC: AuxLat.Parametric, 

715 Aux.RECTIFYING: AuxLat.Rectifying} 

716 

717 

718def _Clenshaw(sinp, Zeta, cs, K): 

719 sz, cz = Zeta._yx # isnormal 

720 # Evaluate sum(c[k] * sin((2*k+2) * Zeta)) if sinp else 

721 # sum(c[k] * cos((2*k+2) * Zeta)) 

722 x = _2cos2x(cz, sz) # 2 * cos(2*Zeta) 

723 if isfinite(x): 

724 U0, U1 = Fsum(), Fsum() 

725 # assert len(cs) == K 

726 for r in _reverange(K): 

727 U1 -= U0 * x + cs[r] 

728 U1, U0 = U0, -U1 

729 # u0*f0(Zeta) - u1*fm1(Zeta) 

730 # f0 = sin(2*Zeta) if sinp else cos(2*Zeta) 

731 # fm1 = 0 if sinp else 1 

732 if sinp: 

733 U0 *= sz * cz * _2_0 

734 else: 

735 U0 *= x * _0_5 

736 U0 -= U1 

737 x = float(U0) 

738 return x 

739 

740 

741def _CXcoeffs(aL): # PYCHOK in .auxilats.__main__ 

742 '''(INTERNAL) Get the C{CX_4}, C{_6} or C{_8} coefficients. 

743 ''' 

744 try: # from pygeodesy.auxilats._CX_x import _coeffs_x as _coeffs 

745 _CX_x = _DOT_(_MODS.auxilats.__name__, _UNDER_('_CX', aL)) 

746 _coeffs = _MODS.getattr(_CX_x, _UNDER_('_coeffs', aL)) 

747 except (AttributeError, ImportError, KeyError, TypeError) as x: 

748 raise AuxError(ALorder=aL, cause=x) 

749 # assert _coeffs.ALorder == aL 

750 # assert _coeffs.n == Aux.len(aL) 

751 return _coeffs 

752 

753 

754def _Newton(tphi, Zeta, _toZeta, **name): 

755 # Newton's method fro AuxLat._fromAux 

756 try: 

757 _lg2 = _log2 

758 _abs = fabs 

759 tz = _abs(Zeta.tan) 

760 tphi = tz / tphi # **) 

761 ltz = _lg2(tz) # **) 

762 ltphi = _lg2(tphi) # **) 

763 ltmin = min(ltphi, MIN_EXP) 

764 ltmax = max(ltphi, MAX_EXP) 

765# auxin = Zeta._AUX 

766 s, n, __2 = 0, 3, _0_5 # n = i + 2 

767 _TOL, _xp2 = _EPSqrt, _exp2 

768 for i in range(1, _TRIPS): # up to 1 Ki! 

769 # _toAux(auxin, AuxPhi(tphi), diff=True) 

770 Z = _toZeta(tphi) 

771 # assert Z._AUX == auxin 

772 t, s_ = Z.tan, s 

773 if t > tz: 

774 ltmax, s = ltphi, +1 

775 elif t < tz: 

776 ltmin, s = ltphi, -1 

777 else: 

778 break 

779 # derivative dtan(Z)/dtan(Phi) 

780 # to dlog(tan(Z))/dlog(tan(Phi)) 

781 d = (ltz - _lg2(t)) * t # **) 

782 if d: 

783 d = d / (Z.diff * tphi) # **) 

784 ltphi += d 

785 tphi = _xp2(ltphi) 

786 if _abs(d) < _TOL: 

787 i += 1 

788 # _toAux(auxin, AuxPhi(tphi), diff=True) 

789 Z = _toZeta(tphi) 

790 tphi -= _over(Z.tan - tz, Z.diff) 

791 break 

792 if (i > n and (s * s_) < 0) or not ltmin < ltphi < ltmax: 

793 s, n = 0, (i + 2) 

794 ltphi = (ltmin + ltmax) * __2 

795 tphi = _xp2(ltphi) 

796 else: 

797 i = _TRIPS 

798 Phi = AuxPhi(tphi, **name).copyquadrant(Zeta) 

799 Phi._iter = i 

800 except (ValueError, ZeroDivisionError): # **) zero t, tphi, tz or Z.diff 

801 Phi = AuxPhi(Zeta, **name) # diff as-as 

802 Phi._iter = 0 

803 # assert Phi._AUX == Aux.PHI 

804 return Phi 

805 

806 

807_f, _u = float, _Ufloats() 

808_1__f3 = -1 / _f(3) # XXX +1 / _f(3) 

809_AR2Coeffs = {4: _u(4 / _f(315), 4 / _f(105), 4 / _f(15), _1__f3), 

810 6: _u(4 / _f(1287), 4 / _f(693), 4 / _f(315), 4 / _f(105), 

811 4 / _f(15), _1__f3), 

812 8: _u(4 / _f(3315), 4 / _f(2145), 4 / _f(1287), 4 / _f(693), 

813 4 / _f(315), 4 / _f(105), 4 / _f(15), _1__f3)} 

814_RRCoeffs = {4: _u(1 / _f(64), _0_25), 

815 6: _u(1 / _f(256), 1 / _f(64), _0_25), 

816 8: _u(25 / _f(16384), 1 / _f(256), 1 / _f(64), _0_25)} # PYCHOK used! 

817del _f, _u, _Ufloats, _1__f3 

818# assert set(_AR2Coeffs.keys()) == set(_RRCoeffs.keys()) 

819 

820__all__ += _ALL_DOCS(AuxLat) 

821 

822# **) MIT License 

823# 

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

825# 

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

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

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

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

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

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

832# 

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

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

835# 

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

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

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

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

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

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

842# OTHER DEALINGS IN THE SOFTWARE.