Coverage for pygeodesy/triaxials/bases.py: 91%

1 

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

3 

4u'''(INTERNAL) Base classes for I{ordered} triaxial ellipsoid classes L{Conformal}, L{Conformal3}, 

5L{Triaxial}, L{Triaxial3} and I{unordered} L{Triaxial_}. 

6 

7Transcoded to pure Python from I{Karney}'s GeographicLib 2.7 C++ classes U{Ellipsoid3<https:// 

8GeographicLib.SourceForge.io/C++/doc/classGeographicLib_1_1Triaxial_1_1Ellipsoid3.html>}, 

9U{Cartesian3<https://GeographicLib.SourceForge.io/C++/doc/classGeographicLib_1_1Triaxial_1_1Cartesian3.html>} and 

10U{Conformal3<https://GeographicLib.SourceForge.io/C++/doc/classGeographicLib_1_1Triaxial_1_1Conformal3.html>}. 

11 

12GeographicLib 2.5.2 C++ class U{JacobiConformal<https://GeographicLib.SourceForge.io/C++/doc/ 

13classGeographicLib_1_1JacobiConformal.html#details>}. 

14 

15Copyright (C) U{Charles Karney<mailto:Karney@Alum.MIT.edu>} (2008-2024, 2025) and licensed under the MIT/X11 License. 

16For more information, see the U{GeographicLib 2.5.2 and 2.7<https://GeographicLib.SourceForge.io/>} documentation. 

17 

18Enum-like C{Lat-/Longitude Kinds (LLK)}, see I{Karney}'s U{coord<https://GeographicLib.SourceForge.io/ 

19C++/doc/classGeographicLib_1_1Triaxial_1_1Cartesian3.html>}: 

20 

21@var LLK.CONFORMAL: Jacobi conformal X and Y projection 

22@var LLK.ELLIPSOIDAL: Ellipsoidal lat-, longitude and heading C{bet}, C{omg}, C{alp} (L{Ang}) 

23@var LLK.GEOCENTRIC: Geocentric lat-, longitude and heading C{phi}", C{lam}" and C{zet} (L{Ang}) 

24@var LLK.GEOCENTRIC_X: Geocentric with pole along major X axis 

25@var LLK.GEODETIC: Geodetic lat-, longitude and heading C{phi}, C{lam} and C{zet} (L{Ang}) 

26@var LLK.GEODETIC_X: Geodetic with pole along major X axis 

27@var LLK.GEODETIC_LON0: Geodetic lat-, longitude I{- lon0} and heading C{phi}, C{lam} and C{zet} (L{Ang}) 

28@var LLK.GEOGRAPHIC = LLK.GEODETIC 

29@var LLK.PARAMETRIC: Parametric lat-, longitude and heading C{phi}', C{lam}' and C{zet} (L{Ang}) 

30@var LLK.PARAMETRIC_X: Parametric with pole along major X axis 

31@var LLK.PLANETODETIC = LLK.GEODETIC 

32@var LLK.PLANETOCENTRIC = LLK.GEOCENTRIC 

33''' 

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

35from __future__ import division as _; del _ # noqa: E702 ; 

36 

37# from pygeodesy.angles import Ang, isAng # _MODS 

38# from pygeodesy.basics import map1 # from .namedTuples 

39from pygeodesy.constants import EPS, EPS0, EPS02, EPS4, INT0, NAN, PI_3, PI4, \ 

40 _EPS2e4, _isfinite, float0_, _1_over, _0_0, \ 

41 _1_0, _N_1_0, _3_0, _4_0 # PYCHOK used! 

42# from pygeodesy.ellipses import Ellipse, _isFlat # _MODS 

43# from pygeodesy.ellipsoids import Ellipsoid, _EWGS84 # _MODS 

44# from pygeodesy.elliptic import Elliptic # _MODS 

45# from pygeodesy.errors import _ValueError, _xkwds # from .utily 

46from pygeodesy.fmath import cbrt, fmean_, hypot, norm2, sqrt0, fabs, sqrt 

47from pygeodesy.fsums import _Fsumf_, fsumf_ 

48# from pygeodesy.internals import typename # _MODS 

49from pygeodesy.interns import _a_, _b_, _c_, _inside_, _not_, _NOTEQUAL_, _null_, \ 

50 _outside_, _scale_, _SPACE_, _spherical_, _x_, _y_, _z_ 

51from pygeodesy.lazily import _ALL_DOCS, _ALL_LAZY, _ALL_MODS as _MODS, _FOR_DOCS 

52from pygeodesy.named import _NamedEnum, _NamedEnumItem, _NamedTuple, _Pass # _MODS 

53from pygeodesy.namedTuples import Vector4Tuple, map1 

54from pygeodesy.props import Property_RO, property_doc_, property_RO, \ 

55 deprecated_method, deprecated_property_RO 

56# from pygeodesy.streprs import Fmt # _MODS 

57from pygeodesy.units import Degrees, Easting, Float, Height, Height_, Meter, \ 

58 Meter2, Meter3, Northing, Radius_, Scalar 

59from pygeodesy.utily import km2m, m2km, _ValueError, _xkwds 

60from pygeodesy.vector3d import _otherV3d, Vector3d 

61 

62# from math import fabs, sqrt # from .fmath 

63 

64__all__ = _ALL_LAZY.triaxials_bases 

65__version__ = '26.02.15' 

66 

67_bet_ = 'bet' # PYCHOK shared 

68_llk_ = 'llk' # PYCHOK shared 

69_KTpFlat = 1.5849625007 

70_MAXIT = 33 # 20 # PYCHOK shared 

71_not_ordered_ = _not_('ordered') 

72_omg_ = 'omg' # PYCHOK shared 

73 

74 

75class Conformal5Tuple(_NamedTuple): # see .Forward4Tuple 

76 '''5-Tuple C{(x, y, z, scale, llk)} with the easting C{x} and 

77 northing C{y} projection, C{scale} or C{NAN} I{but with} 

78 C{z=INT0} I{and kind} C{llk=LLK.CONFORMAL} I{always}. 

79 ''' 

80 _Names_ = (_x_, _y_, _z_, _scale_, _llk_) 

81 _Units_ = ( Easting, Northing, _Pass, Scalar, _Pass) 

82 

83 def __new__(cls, x, y, z=INT0, scale=NAN, llk=None, **kwds): # **iteration_name 

84 args = x, y, (z or INT0), scale, (llk or LLK.CONFORMAL) 

85 return _NamedTuple.__new__(cls, args, **kwds) 

86 

87 

88class _LLK(str): 

89 '''(INTERNAL) Lat-/Longitude Kind. 

90 ''' 

91 def __init__(self, llk): # aka C++ alt 

92 self._X = bool(llk.endswith('_X')) 

93 str.__init__(llk) 

94 

95 

96class LLK(object): 

97 '''Enum-like C{Lat-/Longitude Kinds (LLK)}, see U{coord<https://GeographicLib. 

98 SourceForge.io/C++/doc/classGeographicLib_1_1Triaxial_1_1Cartesian3.html>}. 

99 ''' 

100 CONFORMAL = _LLK('CONFORMAL') 

101 

102 ELLIPSOIDAL = _LLK('ELLIPSOIDAL') # bet, omg, alp 

103 GEOCENTRIC = _LLK('GEOCENTRIC') # phi2p, lam2p, zet 

104 GEOCENTRIC_X = _LLK('GEOCENTRIC_X') 

105 GEODETIC = _LLK('GEODETIC') # phi, lam, zet 

106 GEODETIC_LON0 = _LLK('GEODETIC_LON0') 

107 GEODETIC_X = _LLK('GEODETIC_X') 

108 GEOGRAPHIC = GEODETIC 

109 PARAMETRIC = _LLK('PARAMETRIC') # phi1p, lam1p, zet 

110 PARAMETRIC_X = _LLK('PARAMETRIC_X') 

111 PLANETODETIC = GEODETIC 

112 PLANETOCENTRIC = GEOCENTRIC 

113 

114 _CENTRICS = (GEOCENTRIC, GEOCENTRIC_X, PLANETOCENTRIC) 

115 _DETICS = (GEODETIC, GEODETIC_X, GEODETIC_LON0, GEOGRAPHIC, PLANETODETIC) 

116 _METRICS = (PARAMETRIC, PARAMETRIC_X) 

117 _NOIDAL = (None, ELLIPSOIDAL) 

118# _XCLUDE = (CONFORMAL, GEOGRAPHIC, PLANETOCENTRIC, PLANETODETIC) 

119 

120 def __getitem__(self, name): 

121 llk = self.get(name, None) 

122 if llk is None: 

123 t = _MODS.internals.typename(self) 

124 t = _MODS.streprs.Fmt.SQUARE(t, name) 

125 raise _ValueError(t, name) 

126 return llk 

127 

128 def get(self, name, dflt=None): 

129 '''Get an C{LLK} by C{name}. 

130 ''' 

131 llk = getattr(self, name, None) 

132 return llk if isinstance(llk, _LLK) else dflt 

133 

134 def items(self): 

135 '''Yield all C{LLK (name, value)} pairs. 

136 ''' 

137 for n, llk in LLK.__class__.__dict__.items(): 

138 if isinstance(llk, _LLK): 

139 yield n, llk 

140 

141 def keys(self): 

142 '''Yield all C{LLK} names. 

143 ''' 

144 for n, _ in self.items(): 

145 yield n 

146 

147 def values(self): 

148 '''Yield all C{LLK} values. 

149 ''' 

150 for _, llk in self.items(): 

151 yield llk 

152 

153if not _FOR_DOCS: # PYCHOK force epydoc 

154 LLK = LLK() # singleton 

155del _FOR_DOCS 

156 

157 

158def _HeightINT0(h): 

159 return h if h is INT0 else Height(h=h) 

160 

161 

162class _UnOrderedTriaxialBase(_NamedEnumItem): 

163 '''(INTERNAL) Base class for all I{unordered} triaxial classes. 

164 ''' 

165 _ijk = _kji = None 

166 _unordered = True 

167 

168 def __init__(self, a_triaxial, b=None, c=None, **name): 

169 '''New I{unordered} C{Triaxial_}. 

170 

171 @arg a_triaxial: Large, C{X} semi-axis (C{scalar}, conventionally in 

172 C{meter}) or an other L{Triaxial}, L{Triaxial_} or 

173 L{TriaxialB} instance. 

174 @kwarg b: Middle, C{Y} semi-axis (C{meter}, same units as B{C{a}}), 

175 required if C{B{a_triaxial} is scalar}, ignored otherwise. 

176 @kwarg c: Small, C{Z} semi-axis (C{meter}, like B{C{b}}). 

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

178 

179 @raise TriaxialError: Invalid semi-axis or -axes. 

180 ''' 

181 try: 

182 try: 

183 a = a_triaxial 

184 t = a._abc3 

185 name = _xkwds(name, name=a.name) 

186 except AttributeError: 

187 t = Radius_(a=a), Radius_(b=b), Radius_(c=c) 

188 except (TypeError, ValueError) as x: 

189 raise TriaxialError(a=a, b=b, c=c, cause=x) 

190 if name: 

191 self.name = name 

192 

193 a, b, c = self._abc3 = t 

194 if self._unordered: # == not isinstance(self, Triaxial) 

195 s, _, t = sorted(t) 

196 if not (_isfinite(t) and _isfinite(s) and s > 0): 

197 raise TriaxialError(a=a, b=b, c=c) # txt=_invalid_ 

198 elif not (_isfinite(a) and a >= b >= c > 0): # see TriaxialB 

199 raise TriaxialError(a=a, b=b, c=c, txt=_not_ordered_) 

200 elif not (a > c and self._a2c2 > 0 and self.e2ac > 0): 

201 raise TriaxialError(a=a, c=c, e2ac=self.e2ac, txt=_spherical_) 

202 

203# def __str__(self): # in _NamedEnumItem 

204# return self.toStr() 

205 

206 @Property_RO 

207 def a(self): 

208 '''Get the C{largest, x} semi-axis (C{meter}, conventionally). 

209 ''' 

210 a, _, _ = self._abc3 

211 return a 

212 

213 @Property_RO 

214 def a2(self): 

215 '''Get C{a**2}. 

216 ''' 

217 return self.a**2 

218 

219 @Property_RO 

220 def _a2b2(self): 

221 '''(INTERNAL) Get C{a**2 - b**2} == E_sub_e**2. 

222 ''' 

223 a, b, _ = self._abc3 

224 d = a - b 

225 return (d * (a + b)) if d else _0_0 

226 

227 @Property_RO 

228 def _a2_b2(self): 

229 '''(INTERNAL) Get C{(a / b)**2}. 

230 ''' 

231 a, b, _ = self._abc3 

232 return (a / b)**2 if a != b else _1_0 

233 

234 @Property_RO 

235 def abc3(self): # in geed3solve._a12d 

236 '''Get the semi-axes as 3-tuple C{(a, b, c)}. 

237 ''' 

238 return self._abc3 

239 

240 @Property_RO 

241 def _a2b2c23(self): 

242 '''(INTERNAL) Get 3-tuple C{(a**2, b**2, c**2)}. 

243 ''' 

244 a, b, c = self._abc3 

245 return self.a2, self.b2, self.c2 

246 

247 @Property_RO 

248 def _a2c2(self): 

249 '''(INTERNAL) Get C{a**2 - c**2} == E_sub_x**2. 

250 ''' 

251 a, _, c = self._abc3 

252 d = a - c 

253 return (d * (a + c)) if d else _0_0 

254 

255 @Property_RO 

256 def area(self): 

257 '''Get the surface area (C{meter} I{squared}). 

258 ''' 

259 return self.areaKT(_KTpFlat) if self.isFlat else self.areaRG 

260 

261 def areaKT(self, *p): 

262 '''I{Approximate} the surface area using U{Knud Thomson's 

263 <https://WikiPedia.org/wiki/Ellipsoid#Approximate_formula>} 

264 formula (C{meter} I{squared}). 

265 

266 @arg p: Exponent (C{scalar} > 0), 1.6075 for near-spherical 

267 or 1.5849625007 for "near-flat" triaxials. 

268 ''' 

269 a, b, c = self._abc3 

270 if a == b == c: 

271 a *= a 

272 else: 

273 _p = pow 

274 p = p[0] if p else (_KTpFlat if self.isFlat else 1.6075) 

275 a = _p(fmean_(_p(a * b, p), _p(a * c, p), _p(b * c, p)), _1_over(p)) 

276 return Meter2(areaKT=a * PI4) 

277 

278 @deprecated_method 

279 def area_p(self, p=1.6075): 

280 '''DEPRECATED on 2026-02-15, use method L{areaKT<Triaxial_.areaKT>}.''' 

281 return Meter2(area_p=self.areaKT(p)) 

282 

283 @Property_RO 

284 def areaRG(self): 

285 '''Get the surface area using U{Carlson's symmetric RG 

286 <https://WikiPedia.org/wiki/Ellipsoid#Surface_Area>} 

287 form (C{meter} I{squared}). 

288 ''' 

289 t = self.a2, self.b2, self.c2 

290 r = _MODS.elliptic._rG3(*map(_1_over, t)) 

291 return Meter2(areaRG=self.volume * r * _3_0) 

292 

293 @Property_RO 

294 def b(self): 

295 '''Get the C{middle, y} semi-axis (C{meter}, same units as B{C{a}}). 

296 ''' 

297 _, b, _ = self._abc3 

298 return b 

299 

300 @Property_RO 

301 def b2(self): 

302 '''Get C{b**2}. 

303 ''' 

304 return self.b**2 

305 

306 @Property_RO 

307 def _b2_a2(self): 

308 '''(INTERNAL) Get C{(b / a)**2}. 

309 ''' 

310 a, b, _ = self._abc3 

311 return (b / a)**2 if a != b else _1_0 

312 

313 @Property_RO 

314 def _b2c2(self): 

315 '''(INTERNAL) Get C{b**2 - c**2} == E_sub_y**2. 

316 ''' 

317 _, b, c = self._abc3 

318 d = b - c 

319 return (d * (b + c)) if d else _0_0 

320 

321 @Property_RO 

322 def c(self): 

323 '''Get the C{smallest, z} semi-axis (C{meter}, same units as B{C{a}}). 

324 ''' 

325 _, _, c = self._abc3 

326 return c 

327 

328 @Property_RO 

329 def c2(self): 

330 '''Get C{c**2}. 

331 ''' 

332 return self.c**2 

333 

334 @Property_RO 

335 def _c2_a2(self): 

336 '''(INTERNAL) Get C{(c / a)**2}. 

337 ''' 

338 a, _, c = self._abc3 

339 return (c / a)**2 if a != c else _1_0 

340 

341 @Property_RO 

342 def _c2_b2(self): 

343 '''(INTERNAL) Get C{(c / b)**2}. 

344 ''' 

345 _, b, c = self._abc3 

346 return (c / b)**2 if b != c else _1_0 

347 

348 @Property_RO 

349 def e2ab(self): 

350 '''Get the C{ab} ellipse' I{(1st) eccentricity squared} (C{scalar}), M{1 - (b/a)**2}. 

351 ''' 

352 return Float(e2ab=(_1_0 - self._b2_a2) or _0_0) 

353 

354# _1e2ab = _b2_a2 # == C{1 - e2ab} == C{(b/a)**2} 

355 

356 @Property_RO 

357 def e2ac(self): 

358 '''Get the C{ac} ellipse' I{(1st) eccentricity squared} (C{scalar}), M{1 - (c/a)**2}. 

359 ''' 

360 return Float(e2ac=(_1_0 - self._c2_a2) or _0_0) 

361 

362# _1e2ac = _c2_a2 # == C{1 - e2ac} == C{(c/a)**2} 

363 

364 @Property_RO 

365 def e2bc(self): 

366 '''Get the C{bc} ellipse' I{(1st) eccentricity squared} (C{scalar}), M{1 - (c/b)**2}. 

367 ''' 

368 return Float(e2bc=(_1_0 - self._c2_b2) or _0_0) 

369 

370# _1e2bc = _c2_b2 # == C{1 - e2bc} == C{(c/b)**2} 

371 

372 def hartzell4(self, pov, los=False, **name): 

373 '''Compute the intersection of this triaxial's surface with a Line-Of-Sight 

374 from a Point-Of-View in space. 

375 

376 @see: Function L{hartzell4<triaxials.triaxial5.hartzell4>} for further details. 

377 ''' 

378 return _MODS.triaxials.hartzell4(pov, los=los, tri_biax=self, **name) 

379 

380 def height4(self, x_xyz, y=None, z=None, normal=True, eps=EPS, **name): 

381 '''Compute the projection on and the height above or below this triaxial's surface. 

382 

383 @see: Function L{height4<triaxials.triaxial5.height4>} for further details. 

384 ''' 

385 return _MODS.triaxials.height4(x_xyz, y=y, z=z, tri_biax=self, normal=normal, eps=eps, **name) 

386 

387 @Property_RO 

388 def isFlat(self): 

389 '''Is this triaxial "flat", too pro-/oblate (C{bool})? 

390 ''' 

391 _f = _MODS.ellipses._isFlat 

392 c, b, a = sorted(self._abc3) 

393 return _f(a, c) or _f(b, c) or _f(a, b) 

394 

395 @Property_RO 

396 def isOblate(self): 

397 '''Is this triaxial oblate (C{bool})? 

398 ''' 

399 return not (self.isProlate or self.isSpherical) 

400 

401 @Property_RO 

402 def isOrdered(self): 

403 '''Is this triaxial I{ordered} and I{not spherical} (C{bool})? 

404 ''' 

405 a, b, c = self._abc3 

406 return bool(a >= b > c) # b > c! 

407 

408 @Property_RO 

409 def isProlate(self): 

410 '''Is this triaxial prolate (C{bool})? 

411 ''' 

412 a, b, c = self._abc3 

413 return a < b or b < c or a < c 

414 

415 @Property_RO 

416 def isSpherical(self): 

417 '''Is this triaxial I{spherical} (C{Radius} or INT0)? 

418 ''' 

419 a, b, c = self._abc3 

420 return a if a == b == c else INT0 

421 

422 def _norm2(self, s, c, *a): 

423 '''(INTERNAL) Normalize C{s} and C{c} iff not already. 

424 ''' 

425 if fabs(_hypot2_1(s, c)) > EPS02: 

426 s, c = norm2(s, c) 

427 if a: 

428 s, c = norm2(s * self.b, c * a[0]) 

429 return float0_(s, c) 

430 

431 def normal3d(self, x_xyz, y=None, z=None, length=_1_0): 

432 '''Get a 3-D vector I{on and perpendicular to} this triaxial's surface. 

433 

434 @arg x_xyz: X component (C{scalar}) or a cartesian (C{Cartesian}, 

435 L{Ecef9Tuple}, L{Vector3d}, L{Vector3Tuple} or L{Vector4Tuple}). 

436 @kwarg y: Y component (C{scalar}), required if B{C{x_xyz}} if C{scalar}, ignored 

437 otherwise. 

438 @kwarg z: Z component (C{scalar}), like B{C{y}}. 

439 @kwarg length: Optional, signed length in out-/inward direction (C{scalar}). 

440 

441 @return: A C{Vector3d(x_, y_, z_)} normalized to B{C{length}}, pointing out- 

442 or inward for postive respectively negative B{C{length}}. 

443 

444 @raise TriaxialError: Zero length cartesian or vector. 

445 

446 @note: Cartesian C{(B{x}, B{y}, B{z})} I{must be on} this triaxial's surface, 

447 use method L{Triaxial.sideOf} to validate. 

448 

449 @see: Methods L{Triaxial.height4} and L{Triaxial.sideOf}. 

450 ''' 

451 # n = 2 * (x / a2, y / b2, z / c2) 

452 # == 2 * (x, y * a2 / b2, z * a2 / c2) / a2 # iff ordered 

453 # == 2 * (x, y / _b2_a2, z / _c2_a2) / a2 

454 # == unit(x, y / _b2_a2, z / _c2_a2).times(length) 

455 x, y, z = _otherV3d_(x_xyz, y, z).xyz3 

456 n = Vector3d(x, y / self._b2_a2, 

457 z / self._c2_a2, name__=self.normal3d) 

458 u = n.length 

459 if u < EPS0: 

460 raise TriaxialError(x=x_xyz, y=y, z=z, txt=_null_) 

461 return n.times(length / u) 

462 

463 def normal4(self, x_xyz, y=None, z=None, height=0, normal=True): 

464 '''Compute a cartesian at a height above or below this triaxial's surface. 

465 

466 @arg x_xyz: X component (C{scalar}) or a cartesian (C{Cartesian}, L{Ecef9Tuple}, 

467 L{Vector3d}, L{Vector3Tuple} or L{Vector4Tuple}). 

468 @kwarg y: Y component (C{scalar}), required if B{C{x_xyz}} if C{scalar}, ignored 

469 otherwise. 

470 @kwarg z: Z component (C{scalar}), like B{C{y}}. 

471 

472 @kwarg normal: If C{True}, the B{C{height}} is I{perpendicular, plumb} to the 

473 triaxial's surface, otherwise C{radially} to the center of this 

474 triaxial (C{bool}). 

475 @kwarg name: Optional C{B{name}="normal4"} (C{str}). 

476 

477 @return: L{Vector4Tuple}C{(x, y, z, h)} with the cartesian coordinates C{x}, 

478 C{y} and C{z} and C{h} the I{signed, normal distance} to the triaxial's 

479 surface in C{meter}, conventionally. Positive C{h} indicates, the 

480 cartesian is outside the triaxial, negative C{h} means inside. 

481 

482 @raise TriaxialError: Zero length cartesian or vector. 

483 

484 @note: Cartesian C{(B{x}, B{y}, B{z})} I{must be on} this triaxial's surface, 

485 use method L{Triaxial.sideOf} to validate. 

486 

487 @see: Methods L{Triaxial.normal3d} and L{Triaxial.height4}. 

488 ''' 

489 v, h = _otherV3d_(x_xyz, y, z), Height_(height, low=None) 

490 if h: 

491 if v.length < EPS0: 

492 raise TriaxialError(x=x_xyz, y=y, z=z, txt=_null_) 

493 if normal: 

494 n = self.normal3d(v, length=h) 

495 h = n.length 

496 n += v 

497 else: 

498 h = h / v.length 

499 n = v.times(h + _1_0) 

500 else: 

501 n = v 

502 return Vector4Tuple(n.x, n.y, n.z, h, name__=self.normal4) 

503 

504 def _order3(self, *abc, **reverse): # reverse=False 

505 '''(INTERNAL) Un-/Order C{a}, C{b} and C{c}. 

506 

507 @return: 3-Tuple C{(a, b, c)} ordered by or un-ordered 

508 (reverse-ordered) C{ijk} if C{B{reverse}=True}. 

509 ''' 

510 ijk = self._order_ijk(**reverse) 

511 return _getitems(abc, *ijk) if ijk else abc 

512 

513 def _order3d(self, v, **reverse): # reverse=False 

514 '''(INTERNAL) Un-/Order a C{Vector3d}. 

515 

516 @return: Vector3d(x, y, z) un-/ordered. 

517 ''' 

518 ijk = self._order_ijk(**reverse) 

519 return v.classof(*_getitems(v.xyz3, *ijk)) if ijk else v 

520 

521 @Property_RO 

522 def _ordered4(self): 

523 '''(INTERNAL) Helper for C{_hartzell3} and C{_plumbTo5}. 

524 ''' 

525 def _order2(reverse, a, b, c): 

526 '''(INTERNAL) Un-Order C{a}, C{b} and C{c}. 

527 

528 @return: 2-Tuple C{((a, b, c), ijk)} with C{a} >= C{b} >= C{c} 

529 and C{ijk} a 3-tuple with the initial indices. 

530 ''' 

531 i, j, k = range(3) 

532 if a < b: 

533 a, b, i, j = b, a, j, i 

534 if a < c: 

535 a, c, i, k = c, a, k, i 

536 if b < c: 

537 b, c, j, k = c, b, k, j 

538 # reverse (k, j, i) since (a, b, c) is reversed-sorted 

539 ijk = (k, j, i) if reverse else (None if i < j < k else (i, j, k)) 

540 return (a, b, c), ijk 

541 

542 abc, T = self._abc3, self 

543 if not self.isOrdered: 

544 abc, ijk = _order2(False, *abc) 

545 if ijk: 

546 _, kji = _order2(True, *ijk) 

547 T = _UnOrderedTriaxialBase(*abc) 

548 T._ijk, T._kji = ijk, kji 

549 return abc + (T,) 

550 

551 def _order_ijk(self, reverse=False): 

552 '''(INTERNAL) Get the un-/order indices. 

553 ''' 

554 return self._kji if reverse else self._ijk 

555 

556 @deprecated_property_RO 

557 def perimeter4ab(self): 

558 '''DEPRECATED on 2026.02.09, use property L{Ellipse<pygeodesy.Ellipse>}C{(a, b).perimeter2k_}.''' 

559 a, b, _ = self._abc3 

560 return Meter(perimeter4ab=_MODS.ellipses.Ellipse(a, b).perimeter2k_) 

561 

562 @deprecated_property_RO 

563 def perimeter4ac(self): 

564 '''DEPRECATED on 2026.02.09, use property L{Ellipse<pygeodesy.Ellipse>}C{(a, c).perimeter2k_}.''' 

565 a, _, c = self._abc3 

566 return Meter(perimeter4ac=_MODS.ellipses.Ellipse(a, c).perimeter2k_) 

567 

568 @deprecated_property_RO 

569 def perimeter4bc(self): 

570 '''DEPRECATED on 2026.02.09, use property L{Ellipse<pygeodesy.Ellipse>}C{(b, c).perimeter2k_}.''' 

571 _, b, c = self._abc3 

572 return Meter(perimeter4bc=_MODS.ellipses.Ellipse(b, c).perimeter2k_) 

573 

574 @Property_RO 

575 def R2(self): 

576 '''Get the I{authalic} earth radius (C{meter}), M{sqrt(area / PI4)}. 

577 ''' 

578 r = self.isSpherical 

579 return Meter(R2=r if r else sqrt(self.area / PI4)) # Radius 

580 

581 Rauthalic = R2 

582 

583 @Property_RO 

584 def R3(self): 

585 '''Get the I{volumetric} earth radius (C{meter}), M{(a * b * c)**(1/3)}. 

586 ''' 

587 a, b, c = self._abc3 

588 return Meter(R3=a if a == b == c else cbrt(a * b * c)) # Radius 

589 

590 Rvolumetric = R3 

591 

592 def _radialTo3(self, sbeta, cbeta, somega, comega): 

593 '''(INTERNAL) I{Unordered} helper for C{.height4}. 

594 ''' 

595 def _rphi(a, b, sphi, cphi): 

596 # <https://WikiPedia.org/wiki/Ellipse#Polar_form_relative_to_focus> 

597 # polar form: radius(phi) = a * b / hypot(a * sphi, b * cphi) 

598 return (b / hypot(sphi, b / a * cphi)) if a > b else ( 

599 (a / hypot(cphi, a / b * sphi)) if a < b else a) 

600 

601 sa, ca = self._norm2(sbeta, cbeta) 

602 sb, cb = self._norm2(somega, comega) 

603 

604 a, b, c = self._abc3 

605 if a != b: 

606 a = _rphi(a, b, sb, cb) 

607 if a != c: 

608 c = _rphi(a, c, sa, ca) 

609 t = c * ca 

610 return (t * cb), (t * sb), (c * sa) 

611 

612 def sideOf(self, x_xyz, y=None, z=None, eps=EPS4): 

613 '''Is a cartesian on, above or below the surface of this triaxial? 

614 

615 @arg x_xyz: X component (C{scalar}) or a cartesian (C{Cartesian}, 

616 L{Ecef9Tuple}, L{Vector3d}, L{Vector3Tuple} or L{Vector4Tuple}). 

617 @kwarg y: Y component (C{scalar}), required if B{C{x_xyz}} is C{scalar}, 

618 ignored otherwise. 

619 @kwarg z: Z component (C{scalar}), like B{C{y}}. 

620 @kwarg eps: On-surface tolerance (C{scalar}, distance I{squared}). 

621 

622 @return: C{INT0} if C{(B{x}, B{y}, B{z})} is near this triaxial's surface 

623 within tolerance B{C{eps}}, otherwise the signed, radial distance 

624 I{squared} (C{float}), nega-/positive for in- respectively outside 

625 this triaxial. 

626 

627 @see: Methods L{Triaxial.height4} and L{Triaxial.normal3d}. 

628 ''' 

629 v = _otherV3d_(x_xyz, y, z) 

630 s = fsumf_(_N_1_0, *map(_over02, v.xyz3, self._abc3)) 

631 return INT0 if fabs(s) < eps else s 

632 

633 def _sideOn(self, v, eps=_EPS2e4): 

634 s = self.sideOf(v.xyz, eps=eps) 

635 if s: # PYCHOK no cover 

636 t = _SPACE_((_inside_ if s < 0 else _outside_), self.toRepr()) 

637 raise TriaxialError(eps=eps, sideOf=s, x=v.x, y=v.y, z=v.z, txt=t) 

638 return s 

639 

640 def toEllipsoid(self, **name): 

641 '''Convert this triaxial to a I{biaxial} L{Ellipsoid}, provided 2 axes match. 

642 

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

644 

645 @return: An L{Ellipsoid} with north along this C{Z} axis if C{a == b}, 

646 this C{Y} axis if C{a == c} or this C{X} axis if C{b == c}. 

647 

648 @raise TriaxialError: This C{a != b}, C{b != c} and C{c != a}. 

649 

650 @see: Method L{Ellipsoid.toTriaxial}. 

651 ''' 

652 a, b, c = self._abc3 

653 if a == b: 

654 b = c # N = c-Z 

655 elif b == c: # N = a-X 

656 a, b = b, a 

657 elif a != c: # N = b-Y 

658 t = _SPACE_(_a_, _NOTEQUAL_, _b_, _NOTEQUAL_, _c_) 

659 raise TriaxialError(a=a, b=b, c=c, txt=t) 

660 return _MODS.ellipsoids.Ellipsoid(a, b=b, name=self._name__(name)) 

661 

662 toBiaxial = toEllipsoid 

663 

664 def toStr(self, prec=9, **name): # PYCHOK signature 

665 '''Return this C{Triaxial} as a string. 

666 

667 @kwarg prec: Precision, number of decimal digits (0..9). 

668 @kwarg name: Optional name (C{str}), to override or C{None} 

669 to exclude this triaxial's name. 

670 

671 @return: This C{Triaxial}'s attributes (C{str}). 

672 ''' 

673 T = _UnOrderedTriaxialBase 

674 m = _MODS.triaxials 

675 C = m.Triaxial3B 

676 if isinstance(self, C): 

677 t = T.b, C.e2, C.k2, C.kp2 

678 else: 

679 t = T.a, # props 

680 C = m.ConformalSphere 

681 t += (C.ab, C.bc) if isinstance(self, C) else (T.b, T.c) 

682 C = _Triaxial3Base 

683 t += (C.k2, C.kp2) if isinstance(self, C) else \ 

684 (T.e2ab, T.e2bc, T.e2ac) 

685 for C in (m.Conformal, m.Conformal3): 

686 if isinstance(self, C): 

687 t += C.xyQ2, 

688 break 

689 t += T.volume, T.area, T.R2 

690 return self._instr(prec=prec, props=t, **name) 

691 

692 @Property_RO 

693 def unOrdered(self): 

694 '''Is this triaxial I{un-ordered} and I{not spherical} (C{bool})? 

695 ''' 

696 return not (self.isOrdered or bool(self.isSpherical)) 

697 

698 @Property_RO 

699 def volume(self): 

700 '''Get the volume (C{meter**3}), M{4 / 3 * PI * a * b * c}. 

701 ''' 

702 a, b, c = self._abc3 

703 return Meter3(volume=a * b * c * PI_3 * _4_0) 

704 

705 

706class _OrderedTriaxialBase(_UnOrderedTriaxialBase): 

707 '''(INTERNAL) Base class for all I{ordered} triaxial classes. 

708 ''' 

709 _unordered = False 

710 

711 def __init__(self, a_triaxial, b=None, c=None, **name): 

712 '''New I{ordered} L{Triaxial}, L{Triaxial3}, L{Conformal} or L{Conformal3}. 

713 

714 @arg a_triaxial: Largest semi-axis (C{scalar}, conventionally in C{meter}) 

715 or an other L{Triaxial} or L{Triaxial_} instance. 

716 @kwarg b: Middle semi-axis (C{meter}, same units as B{C{a}}), required 

717 if C{B{a_triaxial} is scalar}, ignored otherwise. 

718 @kwarg c: Smallest semi-axis (C{meter}, like B{C{b}}). 

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

720 

721 @note: The semi-axes must be ordered as C{B{a} >= B{b} >= B{c} > 0} and 

722 must be ellipsoidal, C{B{a} > B{c}}. 

723 

724 @raise TriaxialError: Semi-axes unordered, spherical or invalid. 

725 ''' 

726 _UnOrderedTriaxialBase.__init__(self, a_triaxial, b=b, c=c, **name) 

727 

728 @Property_RO 

729 def _a2b2_a2c2(self): 

730 '''@see: Methods C{.forwardBetaOmega} and property C{._k2_kp2E}. 

731 ''' 

732 s = self._a2c2 

733 if s: 

734 s = self._a2b2 / s 

735 return s or _0_0 

736 

737 @Property_RO 

738 def area(self): 

739 '''Get the surface area (C{meter} I{squared}). 

740 

741 @see: U{Surface area<https://WikiPedia.org/wiki/Ellipsoid#Surface_area>}. 

742 ''' 

743 a, b, c = self._abc3 

744 if a != b and self.e2ac > 0: 

745 a = _UnOrderedTriaxialBase(self).area 

746 else: # a == b > c 

747 a = _MODS.ellipsoids.Ellipsoid(a, b=c).areax 

748 return a 

749 

750# @Property_RO 

751# def area2k1J(self): 

752# '''Get the surface area using elliptic integrals (C{meter} I{squared}). 

753# ''' 

754# a, b, c = self._abc3 

755# if a != b and self.e2ac > 0: 

756# k2, kp2 = t = self._k2_kp2E 

757# if min(t) < 0 or max(t) > _1_0 or self.isFlat: 

758# a = self.areaKT() 

759# else: # k2, kp2 swapped! 

760# # from pygeodesy import asin1, fdot_, PI2 

761# aE = _MODS.elliptic.Elliptic(kp2, _0_0, k2, _1_0) 

762# s = sqrt(self.e2ac) # sin(phi)**2 == 1 - (c/a)**2 

763# t = self._c2_a2 / s # cos(phi)**2 == (c/a)**2 

764# r = asin1(s) # phi = atan2(c/a, s) = asin(s) 

765# c = fdot_(aE.fE(r), s, c / a, c / b, 

766# aE.fF(r), t) 

767# a = Meter2(area2k1J=a * b * c * PI2) 

768# else: # a == b > c 

769# a = _MODS.ellipsoids.Ellipsoid(a, b=c).areax 

770# return a 

771 

772 @Property_RO 

773 def _k2_kp2E(self): 

774 '''(INTERNAL) Get elliptic C{k2} and C{kp2} for C{._xE}, C{._yE} and C{.areaE}. 

775 ''' 

776 # k2 = a2b2 / a2c2 * c2_b2 

777 # kp2 = b2c2 / a2c2 * a2_b2 

778 # b2 = b**2 

779 # xE = Elliptic(k2, -a2b2 / b2, kp2, a2_b2) 

780 # yE = Elliptic(kp2, +b2c2 / b2, k2, c2_b2) 

781 # aE = Elliptic(kp2, 0, k2, 1) 

782 k2 = (self._c2_b2 * self._a2b2_a2c2) or _0_0 

783 kp2 = (self._a2_b2 * self._b2c2 / self._a2c2) if k2 else _1_0 

784 return k2, kp2 

785 

786 def _radialTo3(self, sbeta, cbeta, somega, comega): 

787 '''(INTERNAL) Convert I{ellipsoidal} lat- C{beta} and longitude 

788 C{omega} to a cartesian I{on this triaxial's surface}, also 

789 I{ordered} helper for C{.height4 with normal=False}. 

790 ''' 

791 sa, ca = self._norm2(sbeta, cbeta) 

792 sb, cb = self._norm2(somega, comega) 

793 

794 b2_a2 = self._b2_a2 # == (b/a)**2 

795 c2_a2 = -self._c2_a2 # == -(c/a)**2 

796 a2c2_a2 = self. e2ac # (a**2 - c**2) / a**2 == 1 - (c/a)**2 

797 

798 x2 = _Fsumf_(_1_0, -b2_a2 * sa**2, c2_a2 * ca**2).fover(a2c2_a2) 

799 z2 = _Fsumf_(c2_a2, sb**2, b2_a2 * cb**2).fover(a2c2_a2) 

800 

801 x, y, z = self._abc3 

802 x *= cb * _sqrt0(x2) 

803 y *= ca * sb 

804 z *= sa * _sqrt0(z2) 

805 return x, y, z 

806 

807 

808class _Triaxial3Base(_OrderedTriaxialBase): 

809 '''(INTERNAL) Base class for I{unordered} triaxialC{3} classes. 

810 ''' 

811 _e2_k2_kp2 = None 

812 _Lon0 = None 

813 

814 @Property_RO 

815 def e2(self): 

816 '''Get the I{squared eccentricity} (C{scalar}), M{(a**2 - c**2) / b**2}. 

817 ''' 

818 if self._e2_k2_kp2: 

819 e2, _, _ = self._e2_k2_kp2 

820 else: 

821 e2 = self._a2c2 / self.b2 

822 return Float(e2=e2) 

823