Coverage for pygeodesy/triaxials.py: 94%

615 statements  

« prev     ^ index     » next       coverage.py v7.6.1, created at 2024-11-12 16:17 -0500

1 

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

3 

4u'''Triaxal ellipsoid classes I{ordered} L{Triaxial} and I{unordered} L{Triaxial_} and Jacobi 

5conformal projections L{JacobiConformal} and L{JacobiConformalSpherical}, transcoded from 

6I{Charles Karney}'s C++ class U{JacobiConformal<https://GeographicLib.SourceForge.io/C++/doc/ 

7classGeographicLib_1_1JacobiConformal.html#details>} to pure Python and miscellaneous classes 

8L{BetaOmega2Tuple}, L{BetaOmega3Tuple}, L{Jacobi2Tuple} and L{TriaxialError}. 

9 

10Copyright (C) U{Charles Karney<mailto:Karney@Alum.MIT.edu>} (2008-2023). For more information, 

11see the U{GeographicLib<https://GeographicLib.SourceForge.io>} documentation. 

12 

13@see: U{Geodesics on a triaxial ellipsoid<https://WikiPedia.org/wiki/Geodesics_on_an_ellipsoid# 

14 Geodesics_on_a_triaxial_ellipsoid>} and U{Triaxial coordinate systems and their geometrical 

15 interpretation<https://www.Topo.Auth.GR/wp-content/uploads/sites/111/2021/12/09_Panou.pdf>}. 

16 

17@var Triaxials.Amalthea: Triaxial(name='Amalthea', a=125000, b=73000, c=64000, e2ab=0.658944, e2bc=0.231375493, e2ac=0.737856, volume=2446253479595252, area=93239507787.490371704, area_p=93212299402.670425415) 

18@var Triaxials.Ariel: Triaxial(name='Ariel', a=581100, b=577900, c=577700, e2ab=0.01098327, e2bc=0.000692042, e2ac=0.011667711, volume=812633172614203904, area=4211301462766.580078125, area_p=4211301574065.829589844) 

19@var Triaxials.Earth: Triaxial(name='Earth', a=6378173.435, b=6378103.9, c=6356754.399999999, e2ab=0.000021804, e2bc=0.006683418, e2ac=0.006705077, volume=1083208241574987694080, area=510065911057441.0625, area_p=510065915922713.6875) 

20@var Triaxials.Enceladus: Triaxial(name='Enceladus', a=256600, b=251400, c=248300, e2ab=0.040119337, e2bc=0.024509841, e2ac=0.06364586, volume=67094551514082248, area=798618496278.596679688, area_p=798619018175.109863281) 

21@var Triaxials.Europa: Triaxial(name='Europa', a=1564130, b=1561230, c=1560930, e2ab=0.003704694, e2bc=0.000384275, e2ac=0.004087546, volume=15966575194402123776, area=30663773697323.51953125, area_p=30663773794562.45703125) 

22@var Triaxials.Io: Triaxial(name='Io', a=1829400, b=1819300, c=1815700, e2ab=0.011011391, e2bc=0.003953651, e2ac=0.014921506, volume=25313121117889765376, area=41691875849096.7421875, area_p=41691877397441.2109375) 

23@var Triaxials.Mars: Triaxial(name='Mars', a=3394600, b=3393300, c=3376300, e2ab=0.000765776, e2bc=0.009994646, e2ac=0.010752768, volume=162907283585817247744, area=144249140795107.4375, area_p=144249144150662.15625) 

24@var Triaxials.Mimas: Triaxial(name='Mimas', a=207400, b=196800, c=190600, e2ab=0.09960581, e2bc=0.062015624, e2ac=0.155444317, volume=32587072869017956, area=493855762247.691894531, area_p=493857714107.9375) 

25@var Triaxials.Miranda: Triaxial(name='Miranda', a=240400, b=234200, c=232900, e2ab=0.050915557, e2bc=0.011070811, e2ac=0.061422691, volume=54926187094835456, area=698880863325.756958008, area_p=698881306767.950317383) 

26@var Triaxials.Moon: Triaxial(name='Moon', a=1735550, b=1735324, c=1734898, e2ab=0.000260419, e2bc=0.000490914, e2ac=0.000751206, volume=21886698675223740416, area=37838824729886.09375, area_p=37838824733332.2265625) 

27@var Triaxials.Tethys: Triaxial(name='Tethys', a=535600, b=528200, c=525800, e2ab=0.027441672, e2bc=0.009066821, e2ac=0.036259685, volume=623086233855821440, area=3528073490771.394042969, area_p=3528074261832.738769531) 

28@var Triaxials.WGS84_35: Triaxial(name='WGS84_35', a=6378172, b=6378102, c=6356752.314245179, e2ab=0.00002195, e2bc=0.006683478, e2ac=0.006705281, volume=1083207319768789942272, area=510065621722018.125, area_p=510065626587483.3125) 

29''' 

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

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

32 

33from pygeodesy.basics import isLatLon, isscalar 

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

35 _EPS2e4, _SQRT2_2, float0_, isfinite, isnear1, _over, \ 

36 _0_0, _0_5, _1_0, _N_1_0, _64_0, _4_0 # PYCHOK used! 

37from pygeodesy.datums import Datum, _spherical_datum, _WGS84, Ellipsoid, _EWGS84, Fmt 

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

39# from pygeodesy.elliptic import Elliptic # _MODS 

40from pygeodesy.errors import _AssertionError, _ValueError 

41from pygeodesy.fmath import Fdot, fdot, fmean_, hypot, hypot_, norm2, sqrt0 

42from pygeodesy.fsums import _Fsumf_, fsumf_, fsum1f_ 

43from pygeodesy.interns import NN, _a_, _b_, _beta_, _c_, _distant_, _finite_, \ 

44 _height_, _inside_, _near_, _negative_, _not_, \ 

45 _NOTEQUAL_, _null_, _opposite_, _outside_, _SPACE_, \ 

46 _spherical_, _too_, _x_, _y_ 

47# from pygeodesy.lazily import _ALL_LAZY, _ALL_MODS as _MODS # from .vector3d 

48from pygeodesy.named import _lazyNamedEnumItem as _lazy, _name__, _NamedEnum, \ 

49 _NamedEnumItem, _Pass 

50from pygeodesy.namedTuples import LatLon3Tuple, _NamedTupleTo, Vector3Tuple, \ 

51 Vector4Tuple 

52from pygeodesy.props import Property_RO, property_ROver 

53# from pygeodesy.streprs import Fmt # from .datums 

54from pygeodesy.units import Degrees, Float, Height_, Meter, Meter2, Meter3, \ 

55 Radians, Radius, Scalar_ 

56from pygeodesy.utily import asin1, atan2d, km2m, m2km, SinCos2, sincos2d_ 

57from pygeodesy.vector3d import _otherV3d, Vector3d, _ALL_LAZY, _MODS 

58 

59from math import atan2, fabs, sqrt 

60 

61__all__ = _ALL_LAZY.triaxials 

62__version__ = '24.10.15' 

63 

64_not_ordered_ = _not_('ordered') 

65_omega_ = 'omega' 

66_TRIPS = 359 # Eberly 1074? 

67 

68 

69class _NamedTupleToX(_NamedTupleTo): # in .testNamedTuples 

70 '''(INTERNAL) Base class for L{BetaOmega2Tuple}, 

71 L{BetaOmega3Tuple} and L{Jacobi2Tuple}. 

72 ''' 

73 def _toDegrees(self, name, **toDMS_kwds): 

74 '''(INTERNAL) Convert C{self[0:2]} to L{Degrees} or C{toDMS}. 

75 ''' 

76 return self._toX3U(_NamedTupleTo._Degrees3, Degrees, name, *self, **toDMS_kwds) 

77 

78 def _toRadians(self, name): 

79 '''(INTERNAL) Convert C{self[0:2]} to L{Radians}. 

80 ''' 

81 return self._toX3U(_NamedTupleTo._Radians3, Radians, name, *self) 

82 

83 def _toX3U(self, _X3, U, name, a, b, *c, **kwds): 

84 a, b, s = _X3(self, a, b, **kwds) 

85 if s is None or name: 

86 n = self._name__(name) 

87 s = self.classof(a, b, *c, name=n).reUnit(U, U).toUnits() 

88 return s 

89 

90 

91class BetaOmega2Tuple(_NamedTupleToX): 

92 '''2-Tuple C{(beta, omega)} with I{ellipsoidal} lat- and 

93 longitude C{beta} and C{omega} both in L{Radians} (or 

94 L{Degrees}). 

95 ''' 

96 _Names_ = (_beta_, _omega_) 

97 _Units_ = (_Pass, _Pass) 

98 

99 def toDegrees(self, name=NN, **toDMS_kwds): 

100 '''Convert this L{BetaOmega2Tuple} to L{Degrees} or C{toDMS}. 

101 

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

103 

104 @return: L{BetaOmega2Tuple}C{(beta, omega)} with C{beta} and 

105 C{omega} both in L{Degrees} or as L{toDMS} strings 

106 provided some B{C{toDMS_kwds}} keyword arguments are 

107 specified. 

108 ''' 

109 return self._toDegrees(name, **toDMS_kwds) 

110 

111 def toRadians(self, **name): 

112 '''Convert this L{BetaOmega2Tuple} to L{Radians}. 

113 

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

115 

116 @return: L{BetaOmega2Tuple}C{(beta, omega)} with C{beta} and C{omega} 

117 both in L{Radians}. 

118 ''' 

119 return self._toRadians(name) 

120 

121 

122class BetaOmega3Tuple(_NamedTupleToX): 

123 '''3-Tuple C{(beta, omega, height)} with I{ellipsoidal} lat- and 

124 longitude C{beta} and C{omega} both in L{Radians} (or L{Degrees}) 

125 and the C{height}, rather the (signed) I{distance} to the triaxial's 

126 surface (measured along the radial line to the triaxial's center) 

127 in C{meter}, conventionally. 

128 ''' 

129 _Names_ = BetaOmega2Tuple._Names_ + (_height_,) 

130 _Units_ = BetaOmega2Tuple._Units_ + ( Meter,) 

131 

132 def toDegrees(self, name=NN, **toDMS_kwds): 

133 '''Convert this L{BetaOmega3Tuple} to L{Degrees} or C{toDMS}. 

134 

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

136 

137 @return: L{BetaOmega3Tuple}C{(beta, omega, height)} with 

138 C{beta} and C{omega} both in L{Degrees} or as 

139 L{toDMS} strings provided some B{C{toDMS_kwds}} 

140 keyword arguments are specified. 

141 ''' 

142 return self._toDegrees(name, **toDMS_kwds) 

143 

144 def toRadians(self, **name): 

145 '''Convert this L{BetaOmega3Tuple} to L{Radians}. 

146 

147 @kwarg name: Optional C{B{name}=NN} (C{str}), overriding this name. 

148 

149 @return: L{BetaOmega3Tuple}C{(beta, omega, height)} with C{beta} 

150 and C{omega} both in L{Radians}. 

151 ''' 

152 return self._toRadians(name) 

153 

154 def to2Tuple(self, **name): 

155 '''Reduce this L{BetaOmega3Tuple} to a L{BetaOmega2Tuple}. 

156 

157 @kwarg name: Optional C{B{name}=NN} (C{str}), overriding this name. 

158 ''' 

159 return BetaOmega2Tuple(*self[:2], name=self._name__(name)) 

160 

161 

162class Jacobi2Tuple(_NamedTupleToX): 

163 '''2-Tuple C{(x, y)} with a Jacobi Conformal C{x} and C{y} 

164 projection, both in L{Radians} (or L{Degrees}). 

165 ''' 

166 _Names_ = (_x_, _y_) 

167 _Units_ = (_Pass, _Pass) 

168 

169 def toDegrees(self, name=NN, **toDMS_kwds): 

170 '''Convert this L{Jacobi2Tuple} to L{Degrees} or C{toDMS}. 

171 

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

173 

174 @return: L{Jacobi2Tuple}C{(x, y)} with C{x} and C{y} both 

175 in L{Degrees} or as L{toDMS} strings provided some 

176 B{C{toDMS_kwds}} keyword arguments are specified. 

177 ''' 

178 return self._toDegrees(name, **toDMS_kwds) 

179 

180 def toRadians(self, **name): 

181 '''Convert this L{Jacobi2Tuple} to L{Radians}. 

182 

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

184 

185 @return: L{Jacobi2Tuple}C{(x, y)} with C{x} and C{y} both in L{Radians}. 

186 ''' 

187 return self._toRadians(name) 

188 

189 

190class Triaxial_(_NamedEnumItem): 

191 '''I{Unordered} triaxial ellipsoid and base class. 

192 

193 Triaxial ellipsoids with right-handed semi-axes C{a}, C{b} and C{c}, oriented 

194 such that the large principal ellipse C{ab} is the equator I{Z}=0, I{beta}=0, 

195 while the small principal ellipse C{ac} is the prime meridian, plane I{Y}=0, 

196 I{omega}=0. 

197 

198 The four umbilic points, C{abs}(I{omega}) = C{abs}(I{beta}) = C{PI/2}, lie on 

199 the middle principal ellipse C{bc} in plane I{X}=0, I{omega}=C{PI/2}. 

200 

201 @note: I{Geodetic} C{lat}- and C{lon}gitudes are in C{degrees}, I{geodetic} 

202 C{phi} and C{lam}bda are in C{radians}, but I{ellipsoidal} lat- and 

203 longitude C{beta} and C{omega} are in L{Radians} by default (or in 

204 L{Degrees} if converted). 

205 ''' 

206 _ijk = _kji = None 

207 _unordered = True 

208 

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

210 '''New I{unordered} L{Triaxial_}. 

211 

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

213 C{meter}) or an other L{Triaxial} or L{Triaxial_} instance. 

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

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

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

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

218 

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

220 ''' 

221 try: 

222 try: 

223 a = a_triaxial 

224 t = a._abc3 

225 except AttributeError: 

226 t = Radius(a=a), Radius(b=b), Radius(c=c) 

227 except (TypeError, ValueError) as x: 

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

229 if name: 

230 self.name = name 

231 

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

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

234 s, _, t = sorted(t) 

235 if not (isfinite(t) and s > 0): 

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

237 elif not (isfinite(a) and a >= b >= c > 0): 

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

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

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

241 

242 def __str__(self): 

243 return self.toStr() 

244 

245 @Property_RO 

246 def a(self): 

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

248 ''' 

249 a, _, _ = self._abc3 

250 return a 

251 

252 @Property_RO 

253 def _a2b2(self): 

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

255 ''' 

256 a, b, _ = self._abc3 

257 return ((a - b) * (a + b)) if a != b else _0_0 

258 

259 @Property_RO 

260 def _a2_b2(self): 

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

262 ''' 

263 a, b, _ = self._abc3 

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

265 

266 @Property_RO 

267 def _a2c2(self): 

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

269 ''' 

270 a, _, c = self._abc3 

271 return ((a - c) * (a + c)) if a != c else _0_0 

272 

273 @Property_RO 

274 def area(self): 

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

276 ''' 

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

278 if a > c: 

279 a = Triaxial(a, b, c).area if a > b else \ 

280 Ellipsoid(a, b=c).areax # a == b 

281 else: # a == c == b 

282 a = Meter2(area=a**2 * PI4) 

283 return a 

284 

285 def area_p(self, p=1.6075): 

286 '''I{Approximate} the surface area (C{meter} I{squared}). 

287 

288 @kwarg p: Exponent (C{scalar} > 0), 1.6 for near-spherical or 1.5849625007 

289 for "near-flat" triaxials. 

290 

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

292 ''' 

293 a, b, c = self._abc3 

294 if a == b == c: 

295 a *= a 

296 else: 

297 _p = pow 

298 a = _p(fmean_(_p(a * b, p), _p(a * c, p), _p(b * c, p)), _1_0 / p) 

299 return Meter2(area_p=a * PI4) 

300 

301 @Property_RO 

302 def b(self): 

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

304 ''' 

305 _, b, _ = self._abc3 

306 return b 

307 

308 @Property_RO 

309 def _b2c2(self): 

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

311 ''' 

312 _, b, c = self._abc3 

313 return ((b - c) * (b + c)) if b != c else _0_0 

314 

315 @Property_RO 

316 def c(self): 

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

318 ''' 

319 _, _, c = self._abc3 

320 return c 

321 

322 @Property_RO 

323 def _c2_b2(self): 

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

325 ''' 

326 _, b, c = self._abc3 

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

328 

329 @Property_RO 

330 def e2ab(self): 

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

332 ''' 

333 return Float(e2ab=(_1_0 - self._1e2ab) or _0_0) 

334 

335 @Property_RO 

336 def _1e2ab(self): 

337 '''(INTERNAL) Get C{1 - e2ab} == C{(b/a)**2}. 

338 ''' 

339 a, b, _ = self._abc3 

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

341 

342 @Property_RO 

343 def e2ac(self): 

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

345 ''' 

346 return Float(e2ac=(_1_0 - self._1e2ac) or _0_0) 

347 

348 @Property_RO 

349 def _1e2ac(self): 

350 '''(INTERNAL) Get C{1 - e2ac} == C{(c/a)**2}. 

351 ''' 

352 a, _, c = self._abc3 

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

354 

355 @Property_RO 

356 def e2bc(self): 

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

358 ''' 

359 return Float(e2bc=(_1_0 - self._1e2bc) or _0_0) 

360 

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

362 

363 @property_ROver 

364 def _Elliptic(self): 

365 '''(INTERNAL) Get class L{Elliptic}, I{once}. 

366 ''' 

367 return _MODS.elliptic.Elliptic # overwrite property_ROver 

368 

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

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

371 from a Point-Of-View in space. 

372 

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

374 ''' 

375 return hartzell4(pov, los=los, tri_biax=self, **name) 

376 

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

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

379 

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

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

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

383 otherwise. 

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

385 @kwarg normal: If C{True}, the projection is the I{perpendicular, plumb} to the 

386 triaxial's surface, otherwise the C{radial} line to the center of 

387 this triaxial (C{bool}). 

388 @kwarg eps: Tolerance for root finding and validation (C{scalar}), use a negative 

389 value to skip validation. 

390 @kwarg name: Optional C{B{name}="heigh4"} (C{str}). 

391 

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

393 and C{z} of the projection on or the intersection with and with the height 

394 C{h} above or below the triaxial's surface in C{meter}, conventionally. 

395 

396 @raise TriaxialError: Non-cartesian B{C{xyz}}, invalid B{C{eps}}, no convergence in 

397 root finding or validation failed. 

398 

399 @see: Methods L{Triaxial.normal3d} and L{Ellipsoid.height4} and I{Eberly}'s U{Distance from a 

400 Point to ...<https://www.GeometricTools.com/Documentation/DistancePointEllipseEllipsoid.pdf>}. 

401 ''' 

402 v, r = _otherV3d_(x_xyz, y, z), self.isSpherical 

403 

404 i, h = None, v.length 

405 if h < EPS0: # EPS 

406 x = y = z = _0_0 

407 h -= min(self._abc3) # nearest 

408 elif r: # .isSpherical 

409 x, y, z = v.times(r / h).xyz3 

410 h -= r 

411 else: 

412 x, y, z = v.xyz3 

413 try: 

414 if normal: # plumb to surface 

415 x, y, z, h, i = _plumbTo5(x, y, z, self, eps=eps) 

416 else: # radial to center 

417 x, y, z = self._radialTo3(z, hypot(x, y), y, x) 

418 h = v.minus_(x, y, z).length 

419 except Exception as e: 

420 raise TriaxialError(x=x, y=y, z=z, cause=e) 

421 if h > 0 and self.sideOf(v, eps=EPS0) < 0: 

422 h = -h # inside 

423 n = _name__(name, name__=self.height4) 

424 return Vector4Tuple(x, y, z, h, iteration=i, name=n) 

425 

426 @Property_RO 

427 def isOrdered(self): 

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

429 ''' 

430 a, b, c = self._abc3 

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

432 

433 @Property_RO 

434 def isSpherical(self): 

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

436 ''' 

437 a, b, c = self._abc3 

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

439 

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

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

442 ''' 

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

444 s, c = norm2(s, c) 

445 if a: 

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

447 return float0_(s, c) 

448 

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

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

451 

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

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

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

455 otherwise. 

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

457 @kwarg length: Optional, signed length and in-/outward direction (C{scalar}). 

458 

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

460 outward for neg- respectively positive B{C{length}}. 

461 

462 @raise TriaxialError: Zero length cartesian or vector. 

463 

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

465 method L{Triaxial.sideOf} to validate. 

466 

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

468 ''' 

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

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

471 # == 2 * (x, y / _1e2ab, z / _1e2ac) / a2 

472 # == unit(x, y / _1e2ab, z / _1e2ac).times(length) 

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

474 n = Vector3d(x, y / self._1e2ab, 

475 z / self._1e2ac, name__=self.normal3d) 

476 u = n.length 

477 if u < EPS0: 

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

479 return n.times(length / u) 

480 

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

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

483 

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

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

486 ''' 

487 ijk = self._order_ijk(**reverse) 

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

489 

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

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

492 

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

494 ''' 

495 ijk = self._order_ijk(**reverse) 

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

497 

498 @Property_RO 

499 def _ordered4(self): 

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

501 ''' 

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

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

504 

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

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

507 ''' 

508 i, j, k = 0, 1, 2 # range(3) 

509 if a < b: 

510 a, b, i, j = b, a, j, i 

511 if a < c: 

512 a, c, i, k = c, a, k, i 

513 if b < c: 

514 b, c, j, k = c, b, k, j 

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

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

517 return (a, b, c), ijk 

518 

519 abc, T = self._abc3, self 

520 if not self.isOrdered: 

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

522 if ijk: 

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

524 T = Triaxial_(*abc) 

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

526 return abc + (T,) 

527 

528 def _order_ijk(self, reverse=False): 

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

530 ''' 

531 return self._kji if reverse else self._ijk 

532 

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

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

535 ''' 

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

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

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

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

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

541 

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

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

544 

545 a, b, c = self._abc3 

546 if a != b: 

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

548 if a != c: 

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

550 t = c * ca 

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

552 

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

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

555 

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

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

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

559 ignored otherwise. 

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

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

562 

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

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

565 I{squared} (C{float}), negative for in- or positive for outside 

566 this triaxial. 

567 

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

569 ''' 

570 v = _otherV3d_(x_xyz, y, z) 

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

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

573 

574 def toEllipsoid(self, **name): 

575 '''Convert this triaxial to an L{Ellipsoid}, provided 2 axes match. 

576 

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

578 

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

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

581 

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

583 

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

585 ''' 

586 a, b, c = self._abc3 

587 if a == b: 

588 b = c # N = c-Z 

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

590 a, b = b, a 

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

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

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

594 return Ellipsoid(a, b=b, name=self._name__(name)) 

595 

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

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

598 

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

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

601 to exclude this triaxial's name. 

602 

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

604 ''' 

605 T = Triaxial_ 

606 t = T.a, 

607 J = JacobiConformalSpherical 

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

609 t += T.e2ab, T.e2bc, T.e2ac 

610 J = JacobiConformal 

611 if isinstance(self, J): 

612 t += J.xyQ2, 

613 t += T.volume, T.area 

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

615 

616 @Property_RO 

617 def unOrdered(self): 

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

619 ''' 

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

621 

622 @Property_RO 

623 def volume(self): 

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

625 ''' 

626 a, b, c = self._abc3 

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

628 

629 

630class Triaxial(Triaxial_): 

631 '''I{Ordered} triaxial ellipsoid. 

632 

633 @see: L{Triaxial_} for more information. 

634 ''' 

635 _unordered = False 

636 

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

638 '''New I{ordered} L{Triaxial}. 

639 

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

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

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

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

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

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

646 

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

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

649 

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

651 ''' 

652 Triaxial_.__init__(self, a_triaxial, b=b, c=c, **name) 

653 

654 @Property_RO 

655 def _a2b2_a2c2(self): 

656 '''@see: Methods C{.forwardBetaOmega} and C{._k2_kp2}. 

657 ''' 

658 return self._a2b2 / self._a2c2 

659 

660 @Property_RO 

661 def area(self): 

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

663 

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

665 ''' 

666 a, b, c = self._abc3 

667 if a != b: 

668 kp2, k2 = self._k2_kp2 # swapped! 

669 aE = self._Elliptic(k2, _0_0, kp2, _1_0) 

670 c2 = self._1e2ac # cos(phi)**2 = (c/a)**2 

671 s = sqrt(self.e2ac) # sin(phi)**2 = 1 - c2 

672 r = asin1(s) # phi = atan2(sqrt(c2), s) 

673 b *= fsum1f_(aE.fE(r) * s, (c / a) * (c / b), 

674 aE.fF(r) * c2 / s) 

675 a = Meter2(area=a * b * PI2) 

676 else: # a == b > c 

677 a = Ellipsoid(a, b=c).areax 

678 return a 

679 

680 def forwardBetaOmega(self, beta, omega, height=0, **name): 

681 '''Convert I{ellipsoidal} lat- C{beta}, longitude C{omega} and C{height} 

682 to cartesian. 

683 

684 @arg beta: Ellipsoidal latitude (C{radians} or L{Degrees}). 

685 @arg omega: Ellipsoidal longitude (C{radians} or L{Degrees}). 

686 @arg height: Height above or below the ellipsoid's surface (C{meter}, same 

687 units as this triaxial's C{a}, C{b} and C{c} semi-axes). 

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

689 

690 @return: A L{Vector3Tuple}C{(x, y, z)}. 

691 

692 @see: Method L{Triaxial.reverseBetaOmega} and U{expressions (23-25)<https:// 

693 www.Topo.Auth.GR/wp-content/uploads/sites/111/2021/12/09_Panou.pdf>}. 

694 ''' 

695 if height: 

696 z = self._Height(height) + self.c 

697 if z > 0: 

698 z2 = z**2 

699 x = z * _sqrt0(_1_0 + self._a2c2 / z2) 

700 y = z * _sqrt0(_1_0 + self._b2c2 / z2) 

701 else: 

702 x = y = z = _0_0 

703 else: 

704 x, y, z = self._abc3 

705 if z: # and x and y: 

706 sa, ca = SinCos2(beta) 

707 sb, cb = SinCos2(omega) 

708 

709 r = self._a2b2_a2c2 

710 x *= cb * _sqrt0(ca**2 + sa**2 * r) 

711 y *= ca * sb 

712 z *= sa * _sqrt0(_1_0 - cb**2 * r) 

713 return Vector3Tuple(x, y, z, **name) 

714 

715 def forwardBetaOmega_(self, sbeta, cbeta, somega, comega, **name): 

716 '''Convert I{ellipsoidal} lat- and longitude C{beta} and C{omega} 

717 to cartesian coordinates I{on the triaxial's surface}. 

718 

719 @arg sbeta: Ellipsoidal latitude C{sin(beta)} (C{scalar}). 

720 @arg cbeta: Ellipsoidal latitude C{cos(beta)} (C{scalar}). 

721 @arg somega: Ellipsoidal longitude C{sin(omega)} (C{scalar}). 

722 @arg comega: Ellipsoidal longitude C{cos(omega)} (C{scalar}). 

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

724 

725 @return: A L{Vector3Tuple}C{(x, y, z)} on the surface. 

726 

727 @raise TriaxialError: This triaxial is near-spherical. 

728 

729 @see: Method L{Triaxial.reverseBetaOmega}, U{Triaxial ellipsoid coordinate 

730 system<https://WikiPedia.org/wiki/Geodesics_on_an_ellipsoid# 

731 Triaxial_ellipsoid_coordinate_system>} and U{expressions (23-25)<https:// 

732 www.Topo.Auth.GR/wp-content/uploads/sites/111/2021/12/09_Panou.pdf>}. 

733 ''' 

734 t = self._radialTo3(sbeta, cbeta, somega, comega) 

735 return Vector3Tuple(*t, **name) 

736 

737 def forwardCartesian(self, x_xyz, y=None, z=None, **normal_eps_name): 

738 '''Project a cartesian on this triaxial. 

739 

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

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

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

743 ignored otherwise. 

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

745 @kwarg normal_eps_name: Optional keyword arguments C{B{normal}=True}, 

746 C{B{eps}=EPS} and overriding C{B{name}="height4"} (C{str}), 

747 see method L{Triaxial.height4}. 

748 

749 @see: Method L{Triaxial.height4} for further information and method 

750 L{Triaxial.reverseCartesian} to reverse the projection. 

751 ''' 

752 return self.height4(x_xyz, y, z, **normal_eps_name) 

753 

754 def forwardLatLon(self, lat, lon, height=0, **name): 

755 '''Convert I{geodetic} lat-, longitude and height to cartesian. 

756 

757 @arg lat: Geodetic latitude (C{degrees}). 

758 @arg lon: Geodetic longitude (C{degrees}). 

759 @arg height: Height above the ellipsoid (C{meter}, same units 

760 as this triaxial's C{a}, C{b} and C{c} axes). 

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

762 

763 @return: A L{Vector3Tuple}C{(x, y, z)}. 

764 

765 @see: Method L{Triaxial.reverseLatLon} and U{expressions (9-11)<https:// 

766 www.Topo.Auth.GR/wp-content/uploads/sites/111/2021/12/09_Panou.pdf>}. 

767 ''' 

768 return self._forwardLatLon3(height, name, *sincos2d_(lat, lon)) 

769 

770 def forwardLatLon_(self, slat, clat, slon, clon, height=0, **name): 

771 '''Convert I{geodetic} lat-, longitude and height to cartesian. 

772 

773 @arg slat: Geodetic latitude C{sin(lat)} (C{scalar}). 

774 @arg clat: Geodetic latitude C{cos(lat)} (C{scalar}). 

775 @arg slon: Geodetic longitude C{sin(lon)} (C{scalar}). 

776 @arg clon: Geodetic longitude C{cos(lon)} (C{scalar}). 

777 @arg height: Height above the ellipsoid (C{meter}, same units 

778 as this triaxial's axes C{a}, C{b} and C{c}). 

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

780 

781 @return: A L{Vector3Tuple}C{(x, y, z)}. 

782 

783 @see: Method L{Triaxial.reverseLatLon} and U{expressions (9-11)<https:// 

784 www.Topo.Auth.GR/wp-content/uploads/sites/111/2021/12/09_Panou.pdf>}. 

785 ''' 

786 sa, ca = self._norm2(slat, clat) 

787 sb, cb = self._norm2(slon, clon) 

788 return self._forwardLatLon3(height, name, sa, ca, sb, cb) 

789 

790 def _forwardLatLon3(self, height, name, sa, ca, sb, cb): # name always **name 

791 '''(INTERNAL) Helper for C{.forwardLatLon} and C{.forwardLatLon_}. 

792 ''' 

793 ca_x_sb = ca * sb 

794 h = self._Height(height) 

795 # 1 - (1 - (c/a)**2) * sa**2 - (1 - (b/a)**2) * ca**2 * sb**2 

796 t = fsumf_(_1_0, -self.e2ac * sa**2, -self.e2ab * ca_x_sb**2) 

797 n = self.a / _sqrt0(t) # prime vertical 

798 x = (h + n) * ca * cb 

799 y = (h + n * self._1e2ab) * ca_x_sb 

800 z = (h + n * self._1e2ac) * sa 

801 return Vector3Tuple(x, y, z, **name) 

802 

803 def _Height(self, height): 

804 '''(INTERNAL) Validate a C{height}. 

805 ''' 

806 return Height_(height=height, low=-self.c, Error=TriaxialError) 

807 

808 @Property_RO 

809 def _k2_kp2(self): 

810 '''(INTERNAL) Get C{k2} and C{kp2} for C{._xE}, C{._yE} and C{.area}. 

811 ''' 

812 # k2 = a2b2 / a2c2 * c2_b2 

813 # kp2 = b2c2 / a2c2 * a2_b2 

814 # b2 = b**2 

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

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

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

818 return (self._a2b2_a2c2 * self._c2_b2, 

819 self._b2c2 / self._a2c2 * self._a2_b2) 

820 

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

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

823 C{omega} to cartesian coordinates I{on the triaxial's surface}, 

824 also I{ordered} helper for C{.height4}. 

825 ''' 

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

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

828 

829 b2_a2 = self._1e2ab # == (b/a)**2 

830 c2_a2 = -self._1e2ac # == -(c/a)**2 

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

832 

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

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

835 

836 x, y, z = self._abc3 

837 x *= cb * _sqrt0(x2) 

838 y *= ca * sb 

839 z *= sa * _sqrt0(z2) 

840 return x, y, z 

841 

842 def reverseBetaOmega(self, x_xyz, y=None, z=None, **name): 

843 '''Convert cartesian to I{ellipsoidal} lat- and longitude, C{beta}, C{omega} 

844 and height. 

845 

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

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

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

849 ignored otherwise. 

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

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

852 

853 @return: A L{BetaOmega3Tuple}C{(beta, omega, height)} with C{beta} and 

854 C{omega} in L{Radians} and (radial) C{height} in C{meter}, same 

855 units as this triaxial's axes. 

856 

857 @see: Methods L{Triaxial.forwardBetaOmega} and L{Triaxial.forwardBetaOmega_} 

858 and U{expressions (21-22)<https://www.Topo.Auth.GR/wp-content/uploads/ 

859 sites/111/2021/12/09_Panou.pdf>}. 

860 ''' 

861 v = _otherV3d_(x_xyz, y, z) 

862 a, b, h = self._reverseLatLon3(v, atan2, v, self.forwardBetaOmega_) 

863 return BetaOmega3Tuple(Radians(beta=a), Radians(omega=b), h, **name) 

864 

865 def reverseCartesian(self, x_xyz, y=None, z=None, h=0, normal=True, eps=_EPS2e4, **name): 

866 '''"Unproject" a cartesian on to a cartesion I{off} this triaxial's surface. 

867 

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

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

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

871 ignored otherwise. 

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

873 @arg h: Height above or below this triaxial's surface (C{meter}, same units 

874 as the axes). 

875 @kwarg normal: If C{True}, the height is C{normal} to the surface, otherwise 

876 C{radially} to the center of this triaxial (C{bool}). 

877 @kwarg eps: Tolerance for on-surface test (C{scalar}), see method L{Triaxial.sideOf}. 

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

879 

880 @return: A L{Vector3Tuple}C{(x, y, z)}. 

881 

882 @raise TrialError: Cartesian C{(x, y, z)} not on this triaxial's surface. 

883 

884 @see: Methods L{Triaxial.forwardCartesian} and L{Triaxial.height4}. 

885 ''' 

886 v = _otherV3d_(x_xyz, y, z, **name) 

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

888 if s: # PYCHOK no cover 

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

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

891 

892 if h: 

893 if normal: 

894 v = v.plus(self.normal3d(*v.xyz, length=h)) 

895 elif v.length > EPS0: 

896 v = v.times(_1_0 + (h / v.length)) 

897 return v.xyz # Vector3Tuple 

898 

899 def reverseLatLon(self, x_xyz, y=None, z=None, **name): 

900 '''Convert cartesian to I{geodetic} lat-, longitude and height. 

901 

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

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

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

905 ignored otherwise. 

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

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

908 

909 @return: A L{LatLon3Tuple}C{(lat, lon, height)} with C{lat} and C{lon} 

910 in C{degrees} and (radial) C{height} in C{meter}, same units 

911 as this triaxial's axes. 

912 

913 @see: Methods L{Triaxial.forwardLatLon} and L{Triaxial.forwardLatLon_} 

914 and U{expressions (4-5)<https://www.Topo.Auth.GR/wp-content/uploads/ 

915 sites/111/2021/12/09_Panou.pdf>}. 

916 ''' 

917 v = _otherV3d_(x_xyz, y, z) 

918 s = v.times_(self._1e2ac, # == 1 - e_sub_x**2 

919 self._1e2bc, # == 1 - e_sub_y**2 

920 _1_0) 

921 a, b, h = self._reverseLatLon3(s, atan2d, v, self.forwardLatLon_) 

922 return LatLon3Tuple(Degrees(lat=a), Degrees(lon=b), h, **name) 

923 

924 def _reverseLatLon3(self, s, atan2_, v, forward_): 

925 '''(INTERNAL) Helper for C{.reverseBetOmg} and C{.reverseLatLon}. 

926 ''' 

927 x, y, z = s.xyz3 

928 d = hypot( x, y) 

929 a = atan2_(z, d) 

930 b = atan2_(y, x) 

931 h = v.minus_(*forward_(z, d, y, x)).length 

932 return a, b, h 

933 

934 

935class JacobiConformal(Triaxial): 

936 '''This is a conformal projection of a triaxial ellipsoid to a plane in which the 

937 C{X} and C{Y} grid lines are straight. 

938 

939 Ellipsoidal coordinates I{beta} and I{omega} are converted to Jacobi Conformal 

940 I{y} respectively I{x} separately. Jacobi's coordinates have been multiplied 

941 by C{sqrt(B{a}**2 - B{c}**2) / (2 * B{b})} so that the customary results are 

942 returned in the case of an ellipsoid of revolution. 

943 

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

945 licensed under the MIT/X11 License. 

946 

947 @note: This constructor can I{not be used to specify a sphere}, see alternate 

948 L{JacobiConformalSpherical}. 

949 

950 @see: L{Triaxial}, C++ class U{JacobiConformal<https://GeographicLib.SourceForge.io/ 

951 C++/doc/classGeographicLib_1_1JacobiConformal.html#details>}, U{Jacobi's conformal 

952 projection<https://GeographicLib.SourceForge.io/C++/doc/jacobi.html>} and Jacobi, 

953 C. G. J. I{U{Vorlesungen über Dynamik<https://Books.Google.com/books? 

954 id=ryEOAAAAQAAJ&pg=PA212>}}, page 212ff. 

955 ''' 

956 

957 @Property_RO 

958 def _xE(self): 

959 '''(INTERNAL) Get the x-elliptic function. 

960 ''' 

961 k2, kp2 = self._k2_kp2 

962 # -a2b2 / b2 == (b2 - a2) / b2 == 1 - a2 / b2 == 1 - a2_b2 

963 return self._Elliptic(k2, _1_0 - self._a2_b2, kp2, self._a2_b2) 

964 

965 def xR(self, omega): 

966 '''Compute a Jacobi Conformal C{x} projection. 

967 

968 @arg omega: Ellipsoidal longitude (C{radians} or L{Degrees}). 

969 

970 @return: The C{x} projection (L{Radians}). 

971 ''' 

972 return self.xR_(*SinCos2(omega)) 

973 

974 def xR_(self, somega, comega): 

975 '''Compute a Jacobi Conformal C{x} projection. 

976 

977 @arg somega: Ellipsoidal longitude C{sin(omega)} (C{scalar}). 

978 @arg comega: Ellipsoidal longitude C{cos(omega)} (C{scalar}). 

979 

980 @return: The C{x} projection (L{Radians}). 

981 ''' 

982 s, c = self._norm2(somega, comega, self.a) 

983 return Radians(x=self._xE.fPi(s, c) * self._a2_b2) 

984 

985 @Property_RO 

986 def xyQ2(self): 

987 '''Get the Jacobi Conformal quadrant size (L{Jacobi2Tuple}C{(x, y)}). 

988 ''' 

989 return Jacobi2Tuple(Radians(x=self._a2_b2 * self._xE.cPi), 

990 Radians(y=self._c2_b2 * self._yE.cPi), 

991 name=JacobiConformal.xyQ2.name) 

992 

993 def xyR2(self, beta, omega, **name): 

994 '''Compute a Jacobi Conformal C{x} and C{y} projection. 

995 

996 @arg beta: Ellipsoidal latitude (C{radians} or L{Degrees}). 

997 @arg omega: Ellipsoidal longitude (C{radians} or L{Degrees}). 

998 @kwarg name: Optional name (C{str}), overriding C{B{name}="xyR2"}. 

999 

1000 @return: A L{Jacobi2Tuple}C{(x, y)}. 

1001 ''' 

1002 return self.xyR2_(*(SinCos2(beta) + SinCos2(omega)), 

1003 name=_name__(name, name__=self.xyR2)) 

1004 

1005 def xyR2_(self, sbeta, cbeta, somega, comega, **name): 

1006 '''Compute a Jacobi Conformal C{x} and C{y} projection. 

1007 

1008 @arg sbeta: Ellipsoidal latitude C{sin(beta)} (C{scalar}). 

1009 @arg cbeta: Ellipsoidal latitude C{cos(beta)} (C{scalar}). 

1010 @arg somega: Ellipsoidal longitude C{sin(omega)} (C{scalar}). 

1011 @arg comega: Ellipsoidal longitude C{cos(omega)} (C{scalar}). 

1012 @kwarg name: Optional name (C{str}), overriding C{B{name}="xyR2_"}. 

1013 

1014 @return: A L{Jacobi2Tuple}C{(x, y)}. 

1015 ''' 

1016 return Jacobi2Tuple(self.xR_(somega, comega), 

1017 self.yR_(sbeta, cbeta), 

1018 name=_name__(name, name__=self.xyR2_)) 

1019 

1020 @Property_RO 

1021 def _yE(self): 

1022 '''(INTERNAL) Get the x-elliptic function. 

1023 ''' 

1024 kp2, k2 = self._k2_kp2 # swapped! 

1025 # b2c2 / b2 == (b2 - c2) / b2 == 1 - c2 / b2 == e2bc 

1026 return self._Elliptic(k2, self.e2bc, kp2, self._c2_b2) 

1027 

1028 def yR(self, beta): 

1029 '''Compute a Jacobi Conformal C{y} projection. 

1030 

1031 @arg beta: Ellipsoidal latitude (C{radians} or L{Degrees}). 

1032 

1033 @return: The C{y} projection (L{Radians}). 

1034 ''' 

1035 return self.yR_(*SinCos2(beta)) 

1036 

1037 def yR_(self, sbeta, cbeta): 

1038 '''Compute a Jacobi Conformal C{y} projection. 

1039 

1040 @arg sbeta: Ellipsoidal latitude C{sin(beta)} (C{scalar}). 

1041 @arg cbeta: Ellipsoidal latitude C{cos(beta)} (C{scalar}). 

1042 

1043 @return: The C{y} projection (L{Radians}). 

1044 ''' 

1045 s, c = self._norm2(sbeta, cbeta, self.c) 

1046 return Radians(y=self._yE.fPi(s, c) * self._c2_b2) 

1047 

1048 

1049class JacobiConformalSpherical(JacobiConformal): 

1050 '''An alternate, I{spherical} L{JacobiConformal} projection. 

1051 

1052 @see: L{JacobiConformal} for other and more details. 

1053 ''' 

1054 _ab = _bc = 0 

1055 

1056 def __init__(self, radius_triaxial, ab=0, bc=0, **name): 

1057 '''New L{JacobiConformalSpherical}. 

1058 

1059 @arg radius_triaxial: Radius (C{scalar}, conventionally in 

1060 C{meter}) or an other L{JacobiConformalSpherical}, 

1061 L{JacobiConformal} or ordered L{Triaxial}. 

1062 @kwarg ab: Relative magnitude of C{B{a} - B{b}} (C{meter}, 

1063 same units as C{scalar B{radius}}. 

1064 @kwarg bc: Relative magnitude of C{B{b} - B{c}} (C{meter}, 

1065 same units as C{scalar B{radius}}. 

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

1067 

1068 @raise TriaxialError: Invalid B{C{radius_triaxial}}, negative 

1069 B{C{ab}}, negative B{C{bc}} or C{(B{ab} 

1070 + B{bc})} not positive. 

1071 

1072 @note: If B{C{radius_triaxial}} is a L{JacobiConformalSpherical} 

1073 and if B{C{ab}} and B{C{bc}} are both zero or C{None}, 

1074 the B{C{radius_triaxial}}'s C{ab}, C{bc}, C{a}, C{b} 

1075 and C{c} are copied. 

1076 ''' 

1077 try: 

1078 r = radius_triaxial 

1079 if isinstance(r, Triaxial): # ordered only 

1080 t = r._abc3 

1081 j = isinstance(r, JacobiConformalSpherical) and not bool(ab or bc) 

1082 else: 

1083 t = (Radius(radius=r),) * 3 

1084 j = False 

1085 self._ab = r.ab if j else Scalar_(ab=ab) # low=0 

1086 self._bc = r.bc if j else Scalar_(bc=bc) # low=0 

1087 if (self.ab + self.bc) <= 0: 

1088 raise ValueError(_negative_) 

1089 a, _, c = self._abc3 = t 

1090 if not (a >= c and isfinite(self._a2b2) 

1091 and isfinite(self._a2c2)): 

1092 raise ValueError(_not_(_finite_)) 

1093 except (TypeError, ValueError) as x: 

1094 raise TriaxialError(radius_triaxial=r, ab=ab, bc=bc, cause=x) 

1095 if name: 

1096 self.name = name 

1097 

1098 @Property_RO 

1099 def ab(self): 

1100 '''Get relative magnitude C{a - b} (C{meter}, same units as B{C{a}}). 

1101 ''' 

1102 return self._ab 

1103 

1104 @Property_RO 

1105 def _a2b2(self): 

1106 '''(INTERNAL) Get C{a**2 - b**2} == ab * (a + b). 

1107 ''' 

1108 a, b, _ = self._abc3 

1109 return self.ab * (a + b) 

1110 

1111 @Property_RO 

1112 def _a2c2(self): 

1113 '''(INTERNAL) Get C{a**2 - c**2} == a2b2 + b2c2. 

1114 ''' 

1115 return self._a2b2 + self._b2c2 

1116 

1117 @Property_RO 

1118 def bc(self): 

1119 '''Get relative magnitude C{b - c} (C{meter}, same units as B{C{a}}). 

1120 ''' 

1121 return self._bc 

1122 

1123 @Property_RO 

1124 def _b2c2(self): 

1125 '''(INTERNAL) Get C{b**2 - c**2} == bc * (b + c). 

1126 ''' 

1127 _, b, c = self._abc3 

1128 return self.bc * (b + c) 

1129 

1130 @Property_RO 

1131 def radius(self): 

1132 '''Get radius (C{meter}, conventionally). 

1133 ''' 

1134 return self.a 

1135 

1136 

1137class TriaxialError(_ValueError): 

1138 '''Raised for L{Triaxial} issues. 

1139 ''' 

1140 pass # ... 

1141 

1142 

1143class Triaxials(_NamedEnum): 

1144 '''(INTERNAL) L{Triaxial} registry, I{must} be a sub-class 

1145 to accommodate the L{_LazyNamedEnumItem} properties. 

1146 ''' 

1147 def _Lazy(self, *abc, **name): 

1148 '''(INTERNAL) Instantiate the C{Triaxial}. 

1149 ''' 

1150 a, b, c = map(km2m, abc) 

1151 return Triaxial(a, b, c, **name) 

1152 

1153Triaxials = Triaxials(Triaxial, Triaxial_) # PYCHOK singleton 

1154'''Some pre-defined L{Triaxial}s, all I{lazily} instantiated.''' 

1155# <https://ArxIV.org/pdf/1909.06452.pdf> Table 1 Semi-axes in Km 

1156# <https://www.JPS.NASA.gov/education/images/pdf/ss-moons.pdf> 

1157# <https://link.Springer.com/article/10.1007/s00190-022-01650-9> 

1158_abc84_35 = (_EWGS84.a + 35), (_EWGS84.a - 35), _EWGS84.b 

1159Triaxials._assert( # a (Km) b (Km) c (Km) planet 

1160 Amalthea = _lazy('Amalthea', 125.0, 73.0, _64_0), # Jupiter 

1161 Ariel = _lazy('Ariel', 581.1, 577.9, 577.7), # Uranus 

1162 Earth = _lazy('Earth', 6378.173435, 6378.1039, 6356.7544), 

1163 Enceladus = _lazy('Enceladus', 256.6, 251.4, 248.3), # Saturn 

1164 Europa = _lazy('Europa', 1564.13, 1561.23, 1560.93), # Jupiter 

1165 Io = _lazy('Io', 1829.4, 1819.3, 1815.7), # Jupiter 

1166 Mars = _lazy('Mars', 3394.6, 3393.3, 3376.3), 

1167 Mimas = _lazy('Mimas', 207.4, 196.8, 190.6), # Saturn 

1168 Miranda = _lazy('Miranda', 240.4, 234.2, 232.9), # Uranus 

1169 Moon = _lazy('Moon', 1735.55, 1735.324, 1734.898), # Earth 

1170 Tethys = _lazy('Tethys', 535.6, 528.2, 525.8), # Saturn 

1171 WGS84_35 = _lazy('WGS84_35', *map(m2km, _abc84_35))) 

1172del _abc84_35, _EWGS84 

1173 

1174 

1175def _getitems(items, *indices): 

1176 '''(INTERNAL) Get the C{items} at the given I{indices}. 

1177 

1178 @return: C{Type(items[i] for i in indices)} with 

1179 C{Type = type(items)}, any C{type} having 

1180 the special method C{__getitem__}. 

1181 ''' 

1182 return type(items)(map(items.__getitem__, indices)) 

1183 

1184 

1185def _hartzell3(pov, los, Tun): # in .Ellipsoid.hartzell4, .formy.hartzell 

1186 '''(INTERNAL) Hartzell's "Satellite Line-of-Sight Intersection ...", 

1187 formula from a Point-Of-View to an I{un-/ordered} Triaxial. 

1188 ''' 

1189 def _toUvwV3d(los, pov): 

1190 try: # pov must be LatLon or Cartesian if los is a Los 

1191 v = los.toUvw(pov) 

1192 except (AttributeError, TypeError): 

1193 v = _otherV3d(los=los) 

1194 return v 

1195 

1196 p3 = _otherV3d(pov=pov.toCartesian() if isLatLon(pov) else pov) 

1197 if los is True: # normal 

1198 a, b, c, d, i = _plumbTo5(p3.x, p3.y, p3.z, Tun) 

1199 return type(p3)(a, b, c), d, i 

1200 

1201 u3 = p3.negate() if los is False or los is None else _toUvwV3d(los, pov) 

1202 

1203 a, b, c, T = Tun._ordered4 

1204 

1205 a2 = a**2 # largest, factored out 

1206 b2, p2 = (b**2, T._1e2ab) if b != a else (a2, _1_0) 

1207 c2, q2 = (c**2, T._1e2ac) if c != a else (a2, _1_0) 

1208 

1209 p3 = T._order3d(p3) 

1210 u3 = T._order3d(u3).unit() # unit vector, opposing signs 

1211 

1212 x2, y2, z2 = p3.x2y2z2 # p3.times_(p3).xyz3 

1213 ux, vy, wz = u3.times_(p3).xyz3 

1214 u2, v2, w2 = u3.x2y2z2 # u3.times_(u3).xyz3 

1215 

1216 t = (p2 * c2), c2, b2 

1217 m = fdot(t, u2, v2, w2) # a2 factored out 

1218 if m < EPS0: # zero or near-null LOS vector 

1219 raise _ValueError(_near_(_null_)) 

1220 

1221 r = fsumf_(b2 * w2, c2 * v2, -v2 * z2, vy * wz * 2, 

1222 -w2 * y2, -u2 * y2 * q2, -u2 * z2 * p2, ux * wz * 2 * p2, 

1223 -w2 * x2 * p2, b2 * u2 * q2, -v2 * x2 * q2, ux * vy * 2 * q2) 

1224 if r > 0: # a2 factored out 

1225 r = sqrt(r) * b * c # == a * a * b * c / a2 

1226 elif r < 0: # LOS pointing away from or missing the triaxial 

1227 raise _ValueError(_opposite_ if max(ux, vy, wz) > 0 else _outside_) 

1228 

1229 d = Fdot(t, ux, vy, wz).fadd_(r).fover(m) # -r for antipode, a2 factored out 

1230 if d > 0: # POV inside or LOS outside or missing the triaxial 

1231 s = fsumf_(_N_1_0, _over(x2, a2), _over(y2, b2), _over(z2, c2)) # like _sideOf 

1232 raise _ValueError(_outside_ if s > 0 else _inside_) 

1233 elif fsum1f_(x2, y2, z2, -d*d) < 0: # d past triaxial's center 

1234 raise _ValueError(_too_(_distant_)) 

1235 

1236 v = p3.minus(u3.times(d)) # cartesian type(pov) or Vector3d 

1237 h = p3.minus(v).length # distance to pov == -d 

1238 return T._order3d(v, reverse=True), h, None 

1239 

1240 

1241def hartzell4(pov, los=False, tri_biax=_WGS84, **name): 

1242 '''Compute the intersection of a tri-/biaxial ellipsoid and a Line-Of-Sight from 

1243 a Point-Of-View outside. 

1244 

1245 @arg pov: Point-Of-View outside the tri-/biaxial (C{Cartesian}, L{Ecef9Tuple} 

1246 C{LatLon} or L{Vector3d}). 

1247 @kwarg los: Line-Of-Sight, I{direction} to the tri-/biaxial (L{Los}, L{Vector3d}), 

1248 or C{True} for the I{normal, perpendicular, plumb} to the surface of 

1249 the tri-/biaxial or C{False} or C{None} to point to its center. 

1250 @kwarg tri_biax: A triaxial (L{Triaxial}, L{Triaxial_}, L{JacobiConformal} or 

1251 L{JacobiConformalSpherical}) or biaxial ellipsoid (L{Datum}, 

1252 L{Ellipsoid}, L{Ellipsoid2}, L{a_f2Tuple}) or spherical earth 

1253 radius (C{scalar}, conventionally in C{meter}). 

1254 @kwarg name: Optional name (C{str}), overriding C{B{name}="hartzell4"}. 

1255 

1256 @return: L{Vector4Tuple}C{(x, y, z, h)} on the tri-/biaxial's surface, with C{h} 

1257 the distance from B{C{pov}} to C{(x, y, z)} I{along the} B{C{los}}, all 

1258 in C{meter}, conventionally. 

1259 

1260 @raise TriaxialError: Invalid B{C{pov}} or B{C{pov}} inside the tri-/biaxial or 

1261 invalid B{C{los}} or B{C{los}} points outside or away from 

1262 the tri-/biaxial. 

1263 

1264 @raise TypeError: Invalid B{C{tri_biax}}, C{ellipsoid} or C{datum}. 

1265 

1266 @see: Class L{pygeodesy3.Los}, functions L{pygeodesy.tyr3d} and L{pygeodesy.hartzell} 

1267 and U{lookAtSpheroid<https://PyPI.org/project/pymap3d>} and U{"Satellite 

1268 Line-of-Sight Intersection with Earth"<https://StephenHartzell.Medium.com/ 

1269 satellite-line-of-sight-intersection-with-earth-d786b4a6a9b6>}. 

1270 ''' 

1271 if isinstance(tri_biax, Triaxial_): 

1272 T = tri_biax 

1273 else: 

1274 D = tri_biax if isinstance(tri_biax, Datum) else \ 

1275 _spherical_datum(tri_biax, name__=hartzell4) # _DUNDER_nameof 

1276 T = D.ellipsoid._triaxial 

1277 try: 

1278 v, h, i = _hartzell3(pov, los, T) 

1279 except Exception as x: 

1280 raise TriaxialError(pov=pov, los=los, tri_biax=tri_biax, cause=x) 

1281 return Vector4Tuple(v.x, v.y, v.z, h, iteration=i, # _DUNDER_nameof 

1282 name=_name__(name, name__=hartzell4)) 

1283 

1284 

1285def _hypot2_1(x, y, z=0): 

1286 '''(INTERNAL) Compute M{x**2 + y**2 + z**2 - 1} with C{max(fabs(x), fabs(y), 

1287 fabs(z))} rarely greater than 1.0. 

1288 ''' 

1289 return fsumf_(_N_1_0, x*x, y*y, z*z) 

1290 

1291 

1292def _otherV3d_(x_xyz, y, z, **name): 

1293 '''(INTERNAL) Get a Vector3d from C{x_xyz}, C{y} and C{z}. 

1294 ''' 

1295 return Vector3d(x_xyz, y, z, **name) if isscalar(x_xyz) else \ 

1296 _otherV3d(x_xyz=x_xyz, **name) 

1297 

1298 

1299def _over0(p, q): 

1300 '''(INTERNAL) Return C{p / q} or C{0}. 

1301 ''' 

1302 return (p / q) if q > fabs(p) else _0_0 

1303 

1304 

1305def _over02(p, q): 

1306 '''(INTERNAL) Return C{(p / q)**2} or C{0}. 

1307 ''' 

1308 return (p / q)**2 if p and q else _0_0 

1309 

1310 

1311def _plumbTo3(px, py, E, eps=EPS): # in .ellipsoids.Ellipsoid.height4 

1312 '''(INTERNAL) Nearest point on a 2-D ellipse in 1st quadrant. 

1313 ''' 

1314 a, b = E.a, E.b 

1315 if min(px, py, a, b) < EPS0: 

1316 raise _AssertionError(px=px, py=py, a=a, b=b, E=E) 

1317 

1318 a2 = a - b * E.b_a 

1319 b2 = b - a * E.a_b 

1320 tx = ty = _SQRT2_2 

1321 for i in range(16): # max 5 

1322 ex = a2 * tx**3 

1323 ey = b2 * ty**3 

1324 

1325 qx = px - ex 

1326 qy = py - ey 

1327 q = hypot(qx, qy) 

1328 if q < EPS0: # PYCHOK no cover 

1329 break 

1330 r = hypot(ex - tx * a, 

1331 ey - ty * b) / q 

1332 

1333 sx, tx = tx, min(_1_0, max(0, (ex + qx * r) / a)) 

1334 sy, ty = ty, min(_1_0, max(0, (ey + qy * r) / b)) 

1335 t = hypot(ty, tx) 

1336 if t < EPS0: # PYCHOK no cover 

1337 break 

1338 tx = tx / t # /= chokes PyChecker 

1339 ty = ty / t 

1340 if fabs(sx - tx) < eps and \ 

1341 fabs(sy - ty) < eps: 

1342 break 

1343 

1344 tx *= a / px 

1345 ty *= b / py 

1346 return tx, ty, i # x and y as fractions 

1347 

1348 

1349def _plumbTo4(x, y, a, b, eps=EPS): 

1350 '''(INTERNAL) Nearest point on and distance to a 2-D ellipse, I{unordered}. 

1351 

1352 @see: Function C{_plumbTo3} and I{Eberly}'s U{Distance from a Point to ... an Ellipse ... 

1353 <https://www.GeometricTools.com/Documentation/DistancePointEllipseEllipsoid.pdf>}. 

1354 ''' 

1355 if b > a: 

1356 b, a, d, i = _plumbTo4(y, x, b, a, eps=eps) 

1357 return a, b, d, i 

1358 

1359 if not (b > 0 and isfinite(a)): 

1360 raise _ValueError(a=a, b=b) 

1361 

1362 i = None 

1363 if y: 

1364 if x: 

1365 u = fabs(x / a) 

1366 w = fabs(y / b) 

1367 g = _hypot2_1(u, w) 

1368 if fabs(g) > EPS02: 

1369 r = (b / a)**2 

1370 t, i = _rootNd(_1_0 / r, 0, u, 0, w, g) # eps 

1371 a = _over(x, t * r + _1_0) 

1372 b = _over(y, t + _1_0) 

1373 d = hypot(x - a, y - b) 

1374 else: # on the ellipse 

1375 a, b, d = x, y, _0_0 

1376 else: # x == 0 

1377 if y < 0: 

1378 b = -b 

1379 a = x # signed-0 

1380 d = fabs(y - b) 

1381 

1382 elif x: # y == 0 

1383 d, r = None, _over0(a * x, (a + b) * (a - b)) 

1384 if r: 

1385 a *= r 

1386 r = _1_0 - r**2 

1387 if r > EPS02: 

1388 b *= sqrt(r) 

1389 d = hypot(x - a, y - b) 

1390 elif x < 0: 

1391 a = -a 

1392 if d is None: 

1393 b = y # signed-0 

1394 d = fabs(x - a) 

1395 

1396 else: # x == y == 0 

1397 a = x # signed-0 

1398 d = b 

1399 

1400 return a, b, d, i 

1401 

1402 

1403def _plumbTo5(x, y, z, Tun, eps=EPS): # in .testTriaxials 

1404 '''(INTERNAL) Nearest point on and distance to an I{un-/ordered} triaxial. 

1405 

1406 @see: I{Eberly}'s U{Distance from a Point to ... an Ellipsoid ...<https:// 

1407 www.GeometricTools.com/Documentation/DistancePointEllipseEllipsoid.pdf>}. 

1408 ''' 

1409 a, b, c, T = Tun._ordered4 

1410 if Tun is not T: # T is ordered, Tun isn't 

1411 t = T._order3(x, y, z) + (T,) 

1412 a, b, c, d, i = _plumbTo5(*t, eps=eps) 

1413 return T._order3(a, b, c, reverse=True) + (d, i) 

1414 

1415 if not (c > 0 and isfinite(a)): 

1416 raise _ValueError(a=a, b=b, c=c) 

1417 

1418 if eps > 0: 

1419 val = max(eps * 1e8, EPS) 

1420 else: # no validation 

1421 val, eps = 0, max(EPS0, -eps) 

1422 

1423 i = None 

1424 if z: 

1425 if y: 

1426 if x: 

1427 u = fabs(x / a) 

1428 v = fabs(y / b) 

1429 w = fabs(z / c) 

1430 g = _hypot2_1(u, v, w) 

1431 if fabs(g) > EPS02: 

1432 r = T._1e2ac # (c / a)**2 

1433 s = T._1e2bc # (c / b)**2 

1434 t, i = _rootNd(_1_0 / r, _1_0 / s, u, v, w, g) # eps 

1435 a = _over(x, t * r + _1_0) 

1436 b = _over(y, t * s + _1_0) 

1437 c = _over(z, t + _1_0) 

1438 d = hypot_(x - a, y - b, z - c) 

1439 else: # on the ellipsoid 

1440 a, b, c, d = x, y, z, _0_0 

1441 else: # x == 0 

1442 a = x # 0 

1443 b, c, d, i = _plumbTo4(y, z, b, c, eps=eps) 

1444 elif x: # y == 0 

1445 b = y # 0 

1446 a, c, d, i = _plumbTo4(x, z, a, c, eps=eps) 

1447 else: # x == y == 0 

1448 if z < 0: 

1449 c = -c 

1450 a, b, d = x, y, fabs(z - c) 

1451 

1452 else: # z == 0 

1453 u = _over0(a * x, T._a2c2) # (a + c) * (a - c) 

1454 v = _over0(b * y, T._b2c2) # (b + c) * (b - c) 

1455 s = _hypot2_1(u, v) 

1456 if u and v and s < 0: 

1457 a *= u 

1458 b *= v 

1459 c *= sqrt(-s) 

1460 d = hypot_(x - a, y - b, c) 

1461 else: 

1462 c = z # signed-0 

1463 a, b, d, i = _plumbTo4(x, y, a, b, eps=eps) 

1464 

1465 if val > 0: 

1466 _validate(a, b, c, d, T, x, y, z, val) 

1467 return a, b, c, d, i 

1468 

1469 

1470def _rootNd(r, s, u, v, w, g, eps=EPS0): 

1471 '''(INTERNAL) Robust 2-D or 3-D root finder: 2-D if C{s == v == 0} else 3-D root. 

1472 

1473 @see: I{Eberly}'s U{Robust Root Finders ... and Listing 4<https:// 

1474 www.GeometricTools.com/Documentation/DistancePointEllipseEllipsoid.pdf>}. 

1475 ''' 

1476 u *= r 

1477 v *= s # 0 for 2-D root 

1478 t0 = w + _N_1_0 

1479 t1 = _0_0 if g < 0 else (hypot_(u, w, v) + _N_1_0) 

1480 # assert t0 <= t1 

1481 _e = -eps 

1482 for i in range(1, _TRIPS): # 47..65 

1483 t = (t1 + t0) * _0_5 

1484 e = t1 - t0 

1485 if _e < e < eps or t in (t0, t1): 

1486 break 

1487 g = fsumf_(_N_1_0, # ~= _hypot2_1 

1488 _over02(u, t + r), 

1489 _over02(w, t - _N_1_0), 

1490 _over02(v, t + s) if v else _0_0) 

1491 if g > 0: 

1492 t0 = t 

1493 elif g < 0: 

1494 t1 = t 

1495 else: 

1496 break 

1497 else: # PYCHOK no cover 

1498 t = Fmt.no_convergence(e, eps) 

1499 raise _ValueError(t, txt__=_rootNd) 

1500 return t, i 

1501 

1502 

1503def _sqrt0(x): 

1504 '''(INTERNAL) C{sqrt0} with C{TriaxialError}. 

1505 ''' 

1506 return sqrt0(x, Error=TriaxialError) 

1507 

1508 

1509def _validate(a, b, c, d, T, x, y, z, val): 

1510 '''(INTERNAL) Validate an C{_plumTo5} result. 

1511 ''' 

1512 e = T.sideOf(a, b, c, eps=val) 

1513 if e: # not near the ellipsoid's surface 

1514 raise _ValueError(a=a, b=b, c=c, d=d, 

1515 sideOf=e, eps=val) 

1516 if d: # angle of delta and normal vector 

1517 m = Vector3d(x, y, z).minus_(a, b, c) 

1518 if m.euclid > val: 

1519 m = m.unit() 

1520 n = T.normal3d(a, b, c) 

1521 e = n.dot(m) # n.negate().dot(m) 

1522 if not isnear1(fabs(e), eps1=val): 

1523 raise _ValueError(n=n, m=m, 

1524 dot=e, eps=val) 

1525 

1526 

1527if __name__ == '__main__': 

1528 

1529 from pygeodesy import printf 

1530 from pygeodesy.interns import _COMMA_, _NL_, _NLATvar_ 

1531 

1532 # __doc__ of this file, force all into registery 

1533 t = [NN] + Triaxials.toRepr(all=True, asorted=True).split(_NL_) 

1534 printf(_NLATvar_.join(i.strip(_COMMA_) for i in t)) 

1535 

1536# **) MIT License 

1537# 

1538# Copyright (C) 2022-2024 -- mrJean1 at Gmail -- All Rights Reserved. 

1539# 

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

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

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

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

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

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

1546# 

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

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

1549# 

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

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

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

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

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

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

1556# OTHER DEALINGS IN THE SOFTWARE.