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

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-2024) 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 _ # noqa: E702 ; 

17 

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

19 AuxMu, AuxPhi, AuxTheta, AuxXi 

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

21from pygeodesy.auxilats._CX_Rs import _Rdict, _Rkey, _Rtuple 

22from pygeodesy.basics import _isin, _reverange, _xinstanceof, _passarg 

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

24 _0_0, _0_0s, _0_1, _0_5, _1_0, _2_0, _3_0, _360_0, \ 

25 _log2, _over, isfinite, isinf, isnan 

26from pygeodesy.datums import _ellipsoidal_datum, _WGS84, \ 

27 Ellipsoid, _name__, _EWGS84 

28# from pygeodesy.ellipsoids import Ellipsoid, _EWGS84 # from .datums 

29from pygeodesy.elliptic import Elliptic as _Ef 

30from pygeodesy.errors import AuxError, _xkwds_not, _xkwds_pop2, _Xorder 

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

32from pygeodesy.fsums import Fsum, _Fsumf_, _sum 

33# from pygeodesy.internals import _passarg # from .basics 

34from pygeodesy.interns import NN, _not_scalar_, _UNDER_ 

35from pygeodesy.karney import _2cos2x, _polynomial, _ALL_DOCS, cbrt 

36# from pygeodesy.lazily import _ALL_DOCS # from .karney 

37# from pygeodesy.named import _name__ # from .datums 

38from pygeodesy.props import Property, Property_RO, _update_all 

39from pygeodesy.units import _isDegrees, _isRadius, Degrees, Meter 

40from pygeodesy.utily import atan1, atan2 

41 

42from math import asinh, copysign, cosh, fabs, sin, sinh, sqrt 

43try: 

44 from math import exp2 as _exp2 

45except ImportError: # Python 3.11- 

46 

47 def _exp2(x): 

48 return pow(_2_0, x) 

49 

50__all__ = () 

51__version__ = '26.11.08' 

52 

53_TRIPS = 1024 # XXX 2 or 3? 

54 

55 

56class AuxLat(AuxAngle): 

57 '''Base class for accurate conversion between I{Auxiliary} latitudes 

58 on an ellipsoid. 

59 

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

61 maintain precision near the poles, I{Authalic} latitude C{Xi}, 

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

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

64 

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

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

67 ''' 

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

69 _E = _EWGS84 

70# _Lmax = 0 # overwritten below 

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

72 

73 def __init__(self, a_earth=_EWGS84, f=None, b=None, **ALorder_name): # PYLINT W0231 

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

75 

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

77 datum (L{Datum}, L{Ellipsoid}, L{Ellipsoid2} or L{a_f2Tuple}). 

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

79 required if B{C{a_earth}} is C{scalar} and C{B{b}=None}. 

80 @kwarg b: Optional polar radius, semi-axis (C{meter}, required if B{C{a_earth}} 

81 is C{scalar} and C{B{f}=None}. 

82 @kwarg ALorder_name: Optional C{B{name}=NN} (C{str}) and optional order of 

83 this L{AuxLat} C{B{ALorder}=6}, see property C{ALorder}. 

84 ''' 

85 if ALorder_name: 

86 M = self._mAL 

87 m, name = _xkwds_pop2(ALorder_name, ALorder=M) 

88 if m != M: 

89 self.ALorder = m 

90 else: 

91 name = NN 

92 try: 

93 if not _isin(a_earth, _EWGS84, _WGS84): 

94 n = _name__(name, name__=AuxLat) 

95 if b is f is None: 

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

97 elif _isRadius(a_earth): 

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

99 else: 

100 raise ValueError(_not_scalar_) 

101 self._E = E 

102 elif not (b is f is None): 

103 # turn _UnexpectedError into AuxError 

104 name = _name__(name, **_xkwds_not(None, b=b, f=f)) 

105 

106 if name: 

107 self.name = name 

108 except Exception as x: 

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

110 

111 @Property_RO 

112 def a(self): 

113 '''Get the C{ellipsoid}'s equatorial radius (C{meter}, conventionally). 

114 ''' 

115 return self.ellipsoid.a 

116 

117 equatoradius = a 

118 

119 @Property 

120 def ALorder(self): 

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

122 ''' 

123 return self._mAL 

124 

125 @ALorder.setter # PYCHOK setter! 

126 def ALorder(self, order): 

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

128 ''' 

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

130 if self._mAL != m: 

131 _update_all(self) 

132 self._mAL = m 

133 

134 def _atanhee(self, tphi): # see _Taus.eatanhe, .albers._atanheE 

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

136 f = self.f 

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

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

139 e = self._e 

140 s = _over(atan1(e * s) if f < 0 else asinh(self._e1 * s), e) 

141 return s 

142 

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

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

145 

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

147 @kwarg diff_name: Use C{B{diff}=True} to set C{diff} 

148 optional C{B{name}=NN}. 

149 

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

151 ''' 

152 _xinstanceof(AuxAngle, Phi=Phi) 

153 # assert Phi._AUX == Aux.PHI 

154 tphi = fabs(Phi.tan) 

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

156 y, x = Phi._yx_normalized 

157 q = self._q 

158 qv = self._qf(tphi) 

159 Dq2 = self._Dq(tphi) 

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

161 d, n = _diff_name2(Phi, **diff_name) 

162 Xi = AuxXi(copysign(qv, Phi.y), x * sqrt(Dq2), name=n) 

163 

164 if d: 

165 if isnan(tphi): 

166 d = self._e2m1_sq2 

167 else: 

168 c = self.Parametric(Phi)._x_normalized 

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

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

171 Xi._diff = d 

172 else: 

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

174 # assert Xi._AUX == Aux.XI 

175 return Xi 

176 

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

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

179 

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

181 the I{Taylor} series. 

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

183 

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

185 units as the ellipsoid axes). 

186 

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

188 exceeds C{f_max}. 

189 ''' 

190 f = self.f 

191 if exact or not f: 

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

193 elif fabs(f) < f_max: 

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

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

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

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

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

199 m = self.ALorder 

200 c2 = _polynomial(self._n, _AR2Coeffs[m], 0, m) 

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

202 else: 

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

204 return c2 * _0_5 

205 

206 @Property_RO 

207 def b(self): 

208 '''Get the C{ellipsoid}'s polar radius (C{meter}, conventionally). 

209 ''' 

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

211 

212 polaradius = b 

213 

214 def _coeffs(self, auxout, auxin): 

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

216 aL = self.ALorder # aka Lmax 

217 if auxout == auxin: 

218 return _0_0s(aL) # uncached 

219 

220 k = Aux._1d(auxout, auxin) # PYLINT E1120 

221 try: # cached 

222 return AuxLat._csc[aL][k] 

223 except KeyError: 

224 pass 

225 

226 try: 

227 Cx = self._CXcoeffs[auxout][auxin] # _Rtuple! 

228 except KeyError as x: 

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

230 

231 d = x = n = self._n 

232 if Aux.use_n2(auxin) and Aux.use_n2(auxout): # PYLINT E1120 

233 x = self._n2 

234 

235 def _m(aL): 

236 for m in _reverange(aL): 

237 yield m // 2 

238 else: 

239 _m = _reverange # PYCHOK expected 

240 

241 i = 0 

242 cs = [] 

243 _c = cs.append 

244 _p = _polynomial 

245 for j in _m(aL): 

246 j += i + 1 # order m = j - i - 1 

247 _c(_p(x, Cx, i, j) * d) 

248 d *= n 

249 i = j 

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

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

252 return cs 

253 

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

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

256 

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

258 @kwarg diff_name: Use C{B{diff}=True} to set C{diff} 

259 and an optional C{B{name}=NN}. 

260 

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

262 ''' 

263 _xinstanceof(AuxAngle, Phi=Phi) 

264 # assert Phi._AUX == Aux.PHI 

265 tphi = tchi = fabs(Phi.tan) 

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

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

268 scsig = _sc(sig) 

269 scphi = _sc(tphi) 

270 if self.f > 0: 

271 # The general expression for tchi is 

272 # tphi * scsig - sig * scphi 

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

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

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

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

277 sigtphi = sig / tphi 

278 if sig < (tphi * _0_5): 

279 t = tphi - sig 

280 else: 

281 def _asinh_2(x): 

282 return asinh(x) * _0_5 

283 # Still possibly dangerous cancellation in tphi - sig. 

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

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

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

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

288 # Turn the crank on divided differences, substitute 

289 # sphi = tphi / _sc(tphi) 

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

291 e = self._e 

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

293 # assert em1 != 0 

294 scb = self._scbeta(tphi) 

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

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

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

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

299 try: 

300 Dg = _Fsumf_(atphi, atphib, t, e * t) 

301 except ValueError: # Fsum(NAN) exception 

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

303 e *= atphib 

304 t = atphi - e 

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

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

307 t = float(Dg) # tphi - sig 

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

309 scsig + scphi * sigtphi) if t else _0_0 

310 else: 

311 tchi = tphi * scsig - sig * scphi 

312 

313 d, n = _diff_name2(Phi, **diff_name) 

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

315 

316 if d: 

317 if isinf(tphi): # PYCHOK np cover 

318 d = self._conformal_diff 

319 else: 

320 d = self.Parametric(Phi)._x_normalized 

321 if d: 

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

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

324 Chi._diff = d 

325 # assrt Chi._AUX == Aux.CHI 

326 return Chi 

327 

328 @Property_RO 

329 def _conformal_diff(self): # PYCHOK no cover 

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

331 ''' 

332 e = self._e 

333 if self.f > 0: 

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

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

336 elif e: 

337 ss = sinh(-atan1(e) * e) 

338 d = _sc(ss) - ss 

339 else: 

340 d = _1_0 

341 return d 

342 

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

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

345 Z = d = Zeta_d 

346 if isinstance(Z, AuxAngle): 

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

348 if auxin == auxout: 

349 pass 

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

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

352 elif exact: 

353 p = Aux.power(auxout, auxin) # PYLINT E1120 

354 if p is None: 

355 P = self._fromAux(Z) # Phi 

356 Z = self._toAux(auxout, P) 

357 Z._iter = P.iteration 

358 else: 

359 y, x = Z._yx 

360 if p: 

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

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

363 else: 

364 cs = self._coeffs(auxout, auxin) 

365 yx = Z._yx_normalized 

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

367 # assert Z._yx == yx 

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

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

370 # assert Z._AUX == auxout 

371 return Z 

372 

373 elif _isDegrees(d): 

374 Z = AuxPhi.fromDegrees(d) 

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

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

377 return Degrees(d, name=Aux.Greek(auxout)) # PYLINT E1120 

378 

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

380 

381 @Property_RO 

382 def _CXcoeffs(self): # in .auxilats.__main__, .testAuxilats 

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

384 ''' 

385 return Aux._CXcoeffs(self.ALorder) # PYLINT E1120 

386 

387 def _Dq(self, tphi): 

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

389 sphi = _sn(tphi) 

390 if tphi > 0: 

391 scphi = _sc(tphi) 

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

393 if d: 

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

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

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

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

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

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

400 if e2m1 and ed: 

401 e2 = self._e2 

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

403 scb = self._scbeta(tphi) 

404 q = scphib = scphi / scb 

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

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

407 else: 

408 s2 = sphi * e2 

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

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

411 else: # PYCHOK no cover 

412 q = INF 

413 else: # PYCHOK no cover 

414 q = self._2_e2m12 

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

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

417 return q 

418 

419 @Property_RO 

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

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

422 

423 @Property_RO 

424 def _e1(self): 

425 return sqrt(fabs(self._e12)) 

426 

427 @Property_RO 

428 def _e12(self): 

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

430 

431 @Property_RO 

432 def _e12p1(self): 

433 return _1_0 / self._e2m1 

434 

435 @Property_RO 

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

437 return self.ellipsoid.e2 

438 

439 @Property_RO 

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

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

442 

443 @Property_RO 

444 def _e2m1_sq2(self): 

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

446 

447 @Property_RO 

448 def _2_e2m12(self): 

449 return _2_0 / self._e2m1**2 

450 

451 @Property_RO 

452 def _Ef_fRG_a2b2_PI_4(self): 

453 E = self.ellipsoid 

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

455 

456 @Property_RO 

457 def ellipsoid(self): 

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

459 ''' 

460 return self._E 

461 

462 @Property_RO 

463 def f(self): 

464 '''Get the C{ellipsoid}'s flattening (C{scalar}). 

465 ''' 

466 return self.ellipsoid.f 

467 

468 flattening = f 

469 

470 @Property_RO 

471 def _fm1(self): # 1 - flattening 

472 return self.ellipsoid.f1 

473 

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

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

476 

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

478 @kwarg name: Optional C{B{name}=NN} (C{str}). 

479 

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

481 ''' 

482 _xinstanceof(AuxAngle, Zeta=Zeta) 

483 aux = Zeta._AUX 

484 n = _name__(name, _or_nameof=Zeta) 

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

486 if f is None: 

487 Phi = AuxPhi(NAN, name=n) 

488 elif callable(f): 

489 t = self._fm1 

490 t *= f(t) 

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

492 else: # assert isscalar(f) 

493 y, x = Zeta._yx 

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

495 # assert Phi._AUX == Aux.PHI 

496 return Phi 

497 

498 @Property_RO 

499 def _fromAuxCase(self): 

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

501 ''' 

502 return {Aux.AUTHALIC: cbrt, 

503 Aux.CONFORMAL: _passarg, 

504 Aux.GEOCENTRIC: self._e2m1, 

505 Aux.GEOGRAPHIC: _1_0, 

506 Aux.PARAMETRIC: self._fm1, 

507 Aux.RECTIFYING: sqrt} 

508 

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

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

511 

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

513 @kwarg diff_name: Use C{B{diff}=True} to set C{diff} 

514 and an optional C{B{name}=NN}. 

515 

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

517 ''' 

518 _xinstanceof(AuxAngle, Phi=Phi) 

519 # assert Phi._AUX == Aux.PHI 

520 d, n = _diff_name2(Phi, **diff_name) 

521 Theta = AuxTheta(Phi.y * self._e2m1, Phi.x, name=n) 

522 if d: 

523 Theta._diff = self._e2m1 

524 return Theta 

525 

526 def Geodetic(self, Phi, **name): # PYCHOK no cover 

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

528 

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

530 @kwarg name: Optional C{B{name}=NN} (C{str}). 

531 

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

533 ''' 

534 _xinstanceof(AuxAngle, Phi=Phi) 

535 # assert Phi._AUX == Aux.PHI 

536 _, n = _diff_name2(Phi, **name) 

537 return AuxPhi(Phi, name=n) 

538 

539 @Property_RO 

540 def _n(self): # 3rd flattening 

541 return self.ellipsoid.n 

542 

543 @Property_RO 

544 def _n2(self): 

545 return self._n**2 

546 

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

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

549 

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

551 @kwarg diff_name: Use C{B{diff}=True} to set C{diff} 

552 and an optional C{B{name}=NN}. 

553 

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

555 ''' 

556 _xinstanceof(AuxAngle, Phi=Phi) 

557 # assert Phi._AUX == Aux.PHI 

558 d, n = _diff_name2(Phi, **diff_name) 

559 Beta = AuxBeta(Phi.y * self._fm1, Phi.x, name=n) 

560 if d: 

561 Beta._diff = self._fm1 

562 return Beta 

563 

564 Reduced = Parametric 

565 

566 @Property_RO 

567 def _q(self): # constant _q 

568 q, f = self._e12p1, self.f 

569 if f: 

570 e = self._e 

571 q += _over(asinh(self._e1) if f > 0 else atan1(e), e) 

572 else: 

573 q += _1_0 

574 return q 

575 

576 def _qf(self, tphi): 

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

578 scb = self._scbeta(tphi) 

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

580 

581 def _qIntegrand(self, beta): 

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

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

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

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

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

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

588 # then the integral is 

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

590 # where 

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

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

593 Beta = AuxBeta.fromRadians(beta) 

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

595 Ax, _cv = Aux, self.convert 

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

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

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

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

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

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

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

603 r = _0_0 # XXX 0 avoids NAN summation exceptions 

604 return r 

605 

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

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

608 

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

610 @kwarg diff_name: Use C{B{diff}=True} to set C{diff} 

611 and an optional C{B{name}=NN}. 

612 

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

614 ''' 

615 Beta = self.Parametric(Phi) 

616 # assert Beta._AUX == Aux.BETA 

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

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

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

620 if self.f < 0: 

621 a, b = b, a 

622 ka, kb = kb, ka 

623 ka1, kb1 = kb1, ka1 

624 sb, cb = cb, sb 

625 # now a, b = larger/smaller semi-axis 

626 # Beta measured from larger semi-axis 

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

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

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

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

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

632 sb2 = sb**2 

633 cb2 = cb**2 

634 da2 = ka1 + ka * sb2 

635 # DLMF Eq. 19.25.9 

636 my = b * sb * _Ef._fRF3RD(cb2, _1_0, kb * sb2) 

637 # DLMF Eq. 19.25.10 with complementary angles 

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

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

640 ka * sb / sqrt(da2)) 

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

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

643 if self.f < 0: 

644 a, b = b, a 

645 my, mx = mx, my 

646 mr = (my + mx) / PI_2 

647 if mr: 

648 my = sin(my / mr) 

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

650 else: # zero Mu 

651 my, mx = _0_0, _1_0 

652 d, n = _diff_name2(Phi, **diff_name) 

653 Mu = AuxMu(my, mx, # normalized 

654 name=n).copyquadrant(Phi) 

655 if d: 

656 d, x = _0_0, Beta._x_normalized 

657 if x and mr: 

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

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

660 * (x / Phi._x_normalized) 

661 else: 

662 d = mr / a 

663 Mu._diff = self._fm1 * d 

664 return Mu 

665 

666 def RectifyingRadius(self, exact=False): 

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

668 

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

670 otherwise the I{Taylor} series. 

671 

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

673 as the ellipsoid axes). 

674 ''' 

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

676 return Meter(r, name__=self.RectifyingRadius) 

677 

678 @Property_RO 

679 def _RectifyingR(self): 

680 m = self.ALorder 

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

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

683 

684 def _scbeta(self, tphi): 

685 return _sc(self._fm1 * tphi) 

686 

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

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

689 

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

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

692 @kwarg diff_name: Use C{B{diff}=True} to set C{diff} 

693 and an optional C{B{name}=NN}. 

694 

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

696 ''' 

697 _xinstanceof(AuxAngle, Phi=Phi) 

698 # assert Phi._AUX == Aux.PHI 

699 d, n = _diff_name2(Phi, **diff_name) 

700 m = _toAuxCase.get(auxout, None) 

701 if m: # callable 

702 A = m(self, Phi, diff=d, name=n) 

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

704 A = AuxPhi(Phi, name=n) 

705 else: # auxout? 

706 A = AuxPhi(NAN, name=n) 

707 # assert A._AUX == auxout 

708 return A 

709 

710 def _toZeta(self, zetaux): 

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

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

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

714 ''' 

715 class _AuxPhy(AuxPhi): 

716 # lean C{AuxPhi} instance. 

717 # _diff = _1_0 

718 # _x = _1_0 

719 

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

721 self._y = tphi 

722 

723 m = _toAuxCase.get(zetaux, None) 

724 if m: # callable 

725 

726 def _toZeta(tphi): 

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

728 

729 elif zetaux == Aux.GEODETIC: # GEOGRAPHIC 

730 _toZeta = _AuxPhy 

731 

732 else: # zetaux? 

733 

734 def _toZeta(unused): # PYCHOK expected 

735 return _AuxPhy(NAN) 

736 

737 return _toZeta 

738 

739 

740# switch(auxout): ... 

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

742 Aux.CONFORMAL: AuxLat.Conformal, 

743 Aux.GEOCENTRIC: AuxLat.Geocentric, 

744 Aux.PARAMETRIC: AuxLat.Parametric, 

745 Aux.RECTIFYING: AuxLat.Rectifying} 

746 

747 

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

749 sz, cz = Zeta._yx # isnormal 

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

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

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

753 if isfinite(x): 

754 U0, U1 = Fsum(), Fsum() 

755 # assert len(cs) == K 

756 for r in _reverange(K): 

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

758 U1, U0 = U0, -U1 

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

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

761 # fm1 = 0 if sinp else 1 

762 if sinp: 

763 U0 *= sz * cz * _2_0 

764 else: 

765 U0 *= x * _0_5 

766 U0 -= U1 

767 x = float(U0) 

768 return x 

769 

770 

771def _diff_name2(Phi, diff=False, **name): 

772 '''(INTERNAL) Get C{{Bdiff}=False} and C{B{name}=NN}. 

773 ''' 

774 n = _name__(name, _or_nameof=Phi) # if name else Phi.name 

775 return diff, n 

776 

777 

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

779 # Newton's method from AuxLat._fromAux 

780 try: 

781 _lg2 = _log2 

782 _abs = fabs 

783 tz = _abs(Zeta.tan) 

784 tphi = tz / tphi # **) 

785 ltz = _lg2(tz) # **) 

786 ltphi = _lg2(tphi) # **) 

787 ltmin = min(ltphi, MIN_EXP) 

788 ltmax = max(ltphi, MAX_EXP) 

789# auxin = Zeta._AUX 

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

791 _TOL, _xp2 = _EPSqrt, _exp2 

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

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

794 Z = _toZeta(tphi)