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

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 # _MODS 

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

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

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

41 _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 

53# from pygeodesy.named import _lazyNamedEnumItem as _lazy # _MODS 

54from pygeodesy.namedTuples import Ellipse5Tuple, Vector4Tuple, map1 

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

56 deprecated_method, deprecated_property_RO 

57# from pygeodesy.streprs import Fmt # _MODS 

58from pygeodesy.units import Degrees, Easting, Float, Height, Height_, _Lat0, \ 

59 Meter, Meter2, Meter3, Northing, Radius_, Scalar 

60from pygeodesy.utily import asin1, km2m, m2km, _ValueError, _xkwds 

61from pygeodesy.vector3d import _otherV3d, Vector3d 

62 

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

64 

65__all__ = _ALL_LAZY.triaxials_bases 

66__version__ = '26.03.12' 

67 

68_bet_ = 'bet' # PYCHOK shared 

69_llk_ = 'llk' # PYCHOK shared 

70_KTpFlat = 1.5849625007 

71_MAXIT = 33 # 20 # PYCHOK shared 

72_not_ordered_ = _not_('ordered') 

73_omg_ = 'omg' # PYCHOK shared 

74 

75 

76class Conformal5Tuple(_NamedTuple): # see .Forward4Tuple 

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

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

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

80 ''' 

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

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

83 

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

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

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

87 

88 

89class _LLK(str): 

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

91 ''' 

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

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

94 str.__init__(llk) 

95 

96 

97class LLK(object): 

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

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

100 ''' 

101 CONFORMAL = _LLK('CONFORMAL') 

102 

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

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

105 GEOCENTRIC_X = _LLK('GEOCENTRIC_X') 

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

107 GEODETIC_LON0 = _LLK('GEODETIC_LON0') 

108 GEODETIC_X = _LLK('GEODETIC_X') 

109 GEOGRAPHIC = GEODETIC 

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

111 PARAMETRIC_X = _LLK('PARAMETRIC_X') 

112 PLANETODETIC = GEODETIC 

113 PLANETOCENTRIC = GEOCENTRIC 

114 

115 _CENTRICS = (GEOCENTRIC, GEOCENTRIC_X, PLANETOCENTRIC) 

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

117 _METRICS = (PARAMETRIC, PARAMETRIC_X) 

118 _NOIDAL = (None, ELLIPSOIDAL) 

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

120 

121 def __getitem__(self, name): 

122 llk = self.get(name, None) 

123 if llk is None: 

124 t = _MODS.internals.typename(self) 

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

126 raise _ValueError(t, name) 

127 return llk 

128 

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

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

131 ''' 

132 llk = getattr(self, name, None) 

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

134 

135 def items(self): 

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

137 ''' 

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

139 if isinstance(llk, _LLK): 

140 yield n, llk 

141 

142 def keys(self): 

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

144 ''' 

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

146 yield n 

147 

148 def values(self): 

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

150 ''' 

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

152 yield llk 

153 

154if not _FOR_DOCS: # PYCHOK force epydoc 

155 LLK = LLK() # singleton 

156del _FOR_DOCS 

157 

158 

159def _HeightINT0(h): 

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

161 

162 

163class _UnOrderedTriaxialBase(_NamedEnumItem): 

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

165 ''' 

166 _ijk = _kji = None 

167 _unordered = True 

168 

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

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

171 

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

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

174 L{TriaxialB} instance. 

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

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

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

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

179 

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

181 ''' 

182 try: 

183 try: 

184 a = a_triaxial 

185 t = a._abc3 

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

187 except AttributeError: 

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

189 except (TypeError, ValueError) as x: 

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

191 if name: 

192 self.name = name 

193 

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

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

196 s, _, t = sorted(t) 

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

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

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

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

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

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

203 

204 def __repr__(self): 

205 '''Default C{repr(self)}. 

206 ''' 

207 return self.toRepr(terse=0) 

208 

209# def __str__(self): # in _NamedEnumItem 

210# return self.toStr() 

211 

212 @Property_RO 

213 def a(self): 

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

215 ''' 

216 a, _, _ = self._abc3 

217 return a 

218 

219 @Property_RO 

220 def a2(self): 

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

222 ''' 

223 return self.a**2 

224 

225 @Property_RO 

226 def _a2b2(self): 

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

228 ''' 

229 a, b, _ = self._abc3 

230 d = a - b 

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

232 

233 @Property_RO 

234 def _a2_b2(self): 

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

236 ''' 

237 a, b, _ = self._abc3 

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

239 

240 @Property_RO 

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

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

243 ''' 

244 return self._abc3 

245 

246 @Property_RO 

247 def _a2b2c23(self): # in .triaxials.triaxial3 

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

249 ''' 

250 return self.a2, self.b2, self.c2 

251 

252 @Property_RO 

253 def _a2c2(self): 

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

255 ''' 

256 a, _, c = self._abc3 

257 d = a - c 

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

259 

260 @Property_RO 

261 def area(self): 

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

263 ''' 

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

265 

266 def areaKT(self, *p): 

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

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

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

270 

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

272 or 1.5849625007 for "near-flat" triaxials. 

273 ''' 

274 a, b, c = self._abc3 

275 _p = pow 

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

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

278 return Meter2(areaKT=a * PI4) 

279 

280 @deprecated_method 

281 def area_p(self, p=1.6075): 

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

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

284 

285 @Property_RO 

286 def areaRG(self): 

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

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

289 form (C{meter} I{squared}), see also C{Elliptic.fRG} 

290 ''' 

291 t = sorted(self._a2b2c23) # all non-zero 

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

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

294 

295 @Property_RO 

296 def b(self): 

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

298 ''' 

299 _, b, _ = self._abc3 

300 return b 

301 

302 @Property_RO 

303 def b2(self): 

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

305 ''' 

306 return self.b**2 

307 

308 @Property_RO 

309 def _b2_a2(self): 

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

311 ''' 

312 a, b, _ = self._abc3 

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

314 

315 @Property_RO 

316 def _b2c2(self): 

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

318 ''' 

319 _, b, c = self._abc3 

320 d = b - c 

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

322 

323 @Property_RO 

324 def c(self): 

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

326 ''' 

327 _, _, c = self._abc3 

328 return c 

329 

330 @Property_RO 

331 def c2(self): 

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

333 ''' 

334 return self.c**2 

335 

336 @Property_RO 

337 def _c2_a2(self): 

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

339 ''' 

340 a, _, c = self._abc3 

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

342 

343 @Property_RO 

344 def _c2_b2(self): 

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

346 ''' 

347 _, b, c = self._abc3 

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

349 

350 @Property_RO 

351 def e2ab(self): 

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

353 ''' 

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

355 

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

357 

358 @Property_RO 

359 def e2ac(self): 

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

361 ''' 

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

363 

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

365 

366 @Property_RO 

367 def e2bc(self): 

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

369 ''' 

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

371 

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

373 

374 def ellipse5(self, lat): 

375 '''Get the equatorial or a parallel I{ellipse of lattitude}. 

376 

377 @arg lat: Geodetic latitude (C{degrees90}, C{str} or C{Ang}). 

378 

379 @return: An L{Ellipse5Tuple}C{(a, b, height, lat, beta)} with C{a}, 

380 C{b} and C{height} measured along this trixial's semi-axis 

381 C{a}, C{b} and C{c}, respectively. 

382 

383 @see: Method L{Ellipsoid.circle4<pygeodesy.Ellipsoid.circle4>} for 

384 further details. 

385 ''' 

386 a, b, c = self._abc3 

387 lat = _Lat0(lat) 

388 if lat and c > 0: 

389 E = _MODS.ellipsoids.Ellipsoid 

390 if a > b: 

391 r, z, lat, B = E(a, b=c).circle4(lat) 

392 b *= r / a 

393 a = r 

394 elif b > a: 

395 r, z, lat, B = E(b, b=c).circle4(lat) 

396 a *= r / b 

397 b = r 

398 else: # a == b 

399 r, z, lat, B = E(a, b=c).circle4(lat) 

400 a = b = r 

401 else: # equatorial or "flat" 

402 z = lat = B = _0_0 

403 return Ellipse5Tuple(a, b, z, lat, B) 

404 

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

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

407 from a Point-Of-View in space. 

408 

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

410 ''' 

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

412 

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

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

415 

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

417 ''' 

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

419 

420 @Property_RO 

421 def isFlat(self): 

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

423 ''' 

424 _f = _MODS.ellipses._isFlat 

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

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

427 

428 @Property_RO 

429 def isOblate(self): 

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

431 ''' 

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

433 

434 @Property_RO 

435 def isOrdered(self): 

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

437 ''' 

438 a, b, c = self._abc3 

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

440 

441 @Property_RO 

442 def isProlate(self): 

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

444 ''' 

445 a, b, c = self._abc3 

446 return a < b or b < c or a < c 

447 

448 @Property_RO 

449 def isSpherical(self): 

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

451 ''' 

452 a, b, c = self._abc3 

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

454 

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

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

457 ''' 

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

459 s, c = norm2(s, c) 

460 if a: 

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

462 return float0_(s, c) 

463 

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

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

466 

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

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

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

470 otherwise. 

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

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

473 

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

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

476 

477 @raise TriaxialError: Zero length cartesian or vector. 

478 

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

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

481 

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

483 ''' 

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

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

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

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

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

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

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

491 u = n.length 

492 if u < EPS0: 

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

494 return n.times(length / u) 

495 

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

497 '''Compute a cartesian at a B{C{height}} above or below this triaxial's surface. 

498 

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

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

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

502 otherwise. 

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

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

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

506 triaxial (C{bool}). 

507 

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

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

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

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

512 

513 @raise TriaxialError: Zero length cartesian or vector. 

514 

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

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

517 

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

519 ''' 

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

521 if h: 

522 if v.length < EPS0: 

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

524 if normal: 

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

526 h = n.length 

527 n += v 

528 else: 

529 h = h / v.length # /= chokes PyChecker 

530 n = v.times(h + _1_0) 

531 else: 

532 n = v 

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

534 

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

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

537 

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

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

540 ''' 

541 ijk = self._order_ijk(**reverse) 

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

543 

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

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

546 

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

548 ''' 

549 ijk = self._order_ijk(**reverse) 

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

551 

552 @Property_RO 

553 def _ordered4(self): 

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

555 ''' 

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

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

558 

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

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

561 ''' 

562 i, j, k = range(3) 

563 if a < b: 

564 a, b, i, j = b, a, j, i 

565 if a < c: 

566 a, c, i, k = c, a, k, i 

567 if b < c: 

568 b, c, j, k = c, b, k, j 

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

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

571 return (a, b, c), ijk 

572 

573 abc, T = self._abc3, self 

574 if not self.isOrdered: 

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

576 if ijk: 

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

578 T = _UnOrderedTriaxialBase(*abc) 

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

580 return abc + (T,) 

581 

582 def _order_ijk(self, reverse=False): 

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

584 ''' 

585 return self._kji if reverse else self._ijk 

586 

587 @deprecated_property_RO 

588 def perimeter4ab(self): 

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

590 a, b, _ = self._abc3 

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

592 

593 @deprecated_property_RO 

594 def perimeter4ac(self): 

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

596 a, _, c = self._abc3 

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

598 

599 @deprecated_property_RO 

600 def perimeter4bc(self): 

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

602 _, b, c = self._abc3 

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

604 

605 @Property_RO 

606 def R2(self): 

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

608 ''' 

609 r = self.isSpherical 

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

611 

612 Rauthalic = R2 

613 

614 @Property_RO 

615 def R3(self): 

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

617 ''' 

618 a, b, c = self._abc3 

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

620 

621 Rvolumetric = R3 

622 

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

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

625 ''' 

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

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

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

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

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

631 

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

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

634 

635 a, b, c = self._abc3 

636 if a != b: 

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

638 if a != c: 

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

640 t = c * ca 

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

642 

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

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

645 

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

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

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

649 ignored otherwise. 

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

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

652 

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

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

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

656 this triaxial. 

657 

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

659 ''' 

660 v = _otherV3d_(x_xyz, y, z) 

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

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

663 

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

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

666 if s: # PYCHOK no cover 

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

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

669 return s 

670 

671 def toEllipsoid(self, **name): 

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

673 

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

675 

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

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

678 

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

680 

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

682 ''' 

683 a, b, c = self._abc3 

684 if a == b: 

685 b = c # N = c-Z 

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

687 a, b = b, a 

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

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

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

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

692 

693 toBiaxial = toEllipsoid 

694 

695 def toStr(self, prec=9, terse=-3, **name): # PYCHOK signature 

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

697 

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

699 @kwarg terse: Limit the number of items (C{int}, 3..11), 

700 use C{B{terse}=0} or C{=None} for all. 

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

702 to exclude this triaxial's name. 

703 

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

705 ''' 

706 T = _UnOrderedTriaxialBase 

707 m = _MODS.triaxials 

708 C = m.Triaxial3B 

709 if isinstance(self, C): 

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

711 else: 

712 t = T.a, # props 

713 C = m.ConformalSphere 

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

715 C = _Triaxial3Base 

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

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

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

719 if isinstance(self, C): 

720 t += C.xyQ2, 

721 break 

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

723 if terse: 

724 t = t[:terse] 

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

726 

727 @Property_RO 

728 def unOrdered(self): 

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

730 ''' 

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

732 

733 @Property_RO 

734 def volume(self): 

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

736 ''' 

737 a, b, c = self._abc3 

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

739 

740 

741class _OrderedTriaxialBase(_UnOrderedTriaxialBase): 

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

743 ''' 

744 _unordered = False 

745 

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

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

748 

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

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

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

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

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

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

755 

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

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

758 

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

760 ''' 

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

762 

763 @Property_RO 

764 def _a2b2_a2c2(self): 

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

766 ''' 

767 s = self._a2c2 

768 if s: 

769 s = self._a2b2 / s 

770 return s or _0_0 

771 

772 @Property_RO 

773 def area(self): 

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

775 

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

777 ''' 

778 a = self._areax 

779 if a is None: 

780 a = _UnOrderedTriaxialBase(self).area # or self.area21k 

781 return a 

782 

783 @Property_RO 

784 def area21k(self): 

785 '''Get the surface area using incomplete elliptic integrals of the 

786 2nd and 1st kind (C{meter} I{squared}), see also C{Elliptic.fE} 

787 respectively C{Elliptic.fF}. 

788 ''' 

789 a = self._areax 

790 if a is None: 

791 k2, kp2 = t = self._k2_kp2E 

792 if self.e2ac < EPS or min(t) < EPS or max(t) > _1_0: 

793 a = self.areaKT() # "flat" or near-spherical 

794 else: 

795 aE = _MODS.elliptic.Elliptic(k2=kp2, kp2=k2) # swapped! 

796 s = sqrt(self.e2ac) # == sin(phi) 

797 t = self._c2_a2 / s # == cos(phi)**2 / sin(phi) 

798 r = asin1(s) # phi 

799 a, b, c = self._abc3 

800 t = (aE.fE(r) * s + aE.fF(r) * t) * a * b 

801 a = Meter2(area21k=(c**2 + t) * PI2) 

802 return a 

803 

804 @Property_RO 

805 def _areax(self): 

806 '''(INTERNAL) Get the area as ellipsoidal or C{None}. 

807 ''' 

808 a, b, c = self._abc3 

809 return None if a != b else \ 

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