Coverage for pygeodesy/triaxials.py: 96%

582 statements  

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

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:Charles@Karney.com>} (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 isscalar, map1, _zip, _ValueError 

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

35 _0_0, _0_5, _1_0, _N_2_0, float0_, isfinite, isnear1, \ 

36 _4_0 # PYCHOK used! 

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

38# from pygeodesy.dms import toDMS # _MODS 

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

40# from pygeodesy.elliptic import Elliptic # _MODS 

41# from pygeodesy.errors import _ValueError # from .basics 

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

43from pygeodesy.fsums import Fsum, fsumf_, fsum1f_, Property_RO 

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

45 _height_, _inside_, _near_, _not_, _NOTEQUAL_, _null_, \ 

46 _opposite_, _outside_, _SPACE_, _spherical_, _too_, \ 

47 _x_, _y_ 

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

49from pygeodesy.named import _NamedEnum, _NamedEnumItem, _NamedTuple, _Pass, \ 

50 _lazyNamedEnumItem as _lazy 

51from pygeodesy.namedTuples import LatLon3Tuple, Vector3Tuple, Vector4Tuple 

52# from pygeodesy.props import Property_RO # from .fsums 

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__ = '23.05.22' 

63 

64_not_ordered_ = _not_('ordered') 

65_omega_ = 'omega' 

66_TRIPS = 537 # 52..58, Eberly 1074? 

67 

68 

69class _NamedTupleTo(_NamedTuple): # in .testNamedTuples 

70 '''(INTERNAL) Base for C{-.toDegrees}, C{-.toRadians}. 

71 ''' 

72 def _toDegrees(self, a, b, *c, **toDMS_kwds): 

73 if toDMS_kwds: 

74 toDMS = _MODS.dms.toDMS 

75 a = toDMS(a.toDegrees(), **toDMS_kwds) 

76 b = toDMS(b.toDegrees(), **toDMS_kwds) 

77 elif isinstance(a, Degrees) and \ 

78 isinstance(b, Degrees): 

79 return self 

80 else: 

81 a, b = a.toDegrees(), b.toDegrees() 

82 return self.classof(a, b, *c, name=self.name) 

83 

84 def _toRadians(self, a, b, *c): 

85 return self if isinstance(a, Radians) and \ 

86 isinstance(b, Radians) else \ 

87 self.classof(a.toRadians(), b.toRadians(), 

88 *c, name=self.name) 

89 

90 

91class BetaOmega2Tuple(_NamedTupleTo): 

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, **toDMS_kwds): 

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

101 

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

103 C{beta} and C{omega} both in L{Degrees} 

104 or as a L{toDMS} string provided some 

105 B{C{toDMS_kwds}} keyword arguments are 

106 specified. 

107 ''' 

108 return _NamedTupleTo._toDegrees(self, *self, **toDMS_kwds) 

109 

110 def toRadians(self): 

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

112 

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

114 C{beta} and C{omega} both in L{Radians}. 

115 ''' 

116 return _NamedTupleTo._toRadians(self, *self) 

117 

118 

119class BetaOmega3Tuple(_NamedTupleTo): 

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

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

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

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

124 in C{meter}, conventionally. 

125 ''' 

126 _Names_ = BetaOmega2Tuple._Names_ + (_height_,) 

127 _Units_ = BetaOmega2Tuple._Units_ + ( Meter,) 

128 

129 def toDegrees(self, **toDMS_kwds): 

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

131 

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

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

134 L{toDMS} string provided some B{C{toDMS_kwds}} 

135 keyword arguments are specified. 

136 ''' 

137 return _NamedTupleTo._toDegrees(self, *self, **toDMS_kwds) 

138 

139 def toRadians(self): 

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

141 

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

143 C{beta} and C{omega} both in L{Radians}. 

144 ''' 

145 return _NamedTupleTo._toRadians(self, *self) 

146 

147 def to2Tuple(self): 

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

149 ''' 

150 return BetaOmega2Tuple(*self[:2]) 

151 

152 

153class Jacobi2Tuple(_NamedTupleTo): 

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

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

156 ''' 

157 _Names_ = (_x_, _y_) 

158 _Units_ = (_Pass, _Pass) 

159 

160 def toDegrees(self, **toDMS_kwds): 

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

162 

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

164 both in L{Degrees} or as a L{toDMS} string 

165 provided some B{C{toDMS_kwds}} keyword 

166 arguments are specified. 

167 ''' 

168 return _NamedTupleTo._toDegrees(self, *self, **toDMS_kwds) 

169 

170 def toRadians(self): 

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

172 

173 @return: L{Jacobi2Tuple}C{(x, y)} with C{x} 

174 and C{y} both in L{Radians}. 

175 ''' 

176 return _NamedTupleTo._toRadians(self, *self) 

177 

178 

179class Triaxial_(_NamedEnumItem): 

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

181 

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

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

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

185 I{omega}=0. 

186 

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

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

189 

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

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

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

193 L{Degrees} if converted). 

194 ''' 

195 _ijk = _kji = None 

196 _unordered = True 

197 

198 def __init__(self, a_triaxial, b=None, c=None, name=NN): 

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

200 

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

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

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

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

205 @kwarg c: Small, C{Z} semi-axis (C{meter}, same units as B{C{a}}), required 

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

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

208 

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

210 ''' 

211 try: 

212 a = a_triaxial 

213 t = a._abc3 if isinstance(a, Triaxial_) else ( 

214 Radius(a=a), Radius(b=b), Radius(c=c)) 

215 except (TypeError, ValueError) as x: 

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

217 if name: 

218 self.name = name 

219 

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

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

222 s, _, t = sorted(t) 

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

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

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

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

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

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

229 

230 def __str__(self): 

231 return self.toStr() 

232 

233 @Property_RO 

234 def a(self): 

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

236 ''' 

237 a, _, _ = self._abc3 

238 return a 

239 

240 @Property_RO 

241 def _a2b2(self): 

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

243 ''' 

244 a, b, _ = self._abc3 

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

246 

247 @Property_RO 

248 def _a2_b2(self): 

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

250 ''' 

251 a, b, _ = self._abc3 

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

253 

254 @Property_RO 

255 def _a2c2(self): 

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

257 ''' 

258 a, _, c = self._abc3 

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

260 

261 @Property_RO 

262 def area(self): 

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

264 ''' 

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

266 if a > c: 

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

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

269 else: # a == c == b 

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

271 return a 

272 

273 def area_p(self, p=1.6075): 

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

275 

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

277 for "near-flat" triaxials. 

278 

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

280 ''' 

281 a, b, c = self._abc3 

282 if a == b == c: 

283 a *= a 

284 else: 

285 _p = pow 

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

287 return Meter2(area_p=a * PI4) 

288 

289 @Property_RO 

290 def b(self): 

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

292 ''' 

293 _, b, _ = self._abc3 

294 return b 

295 

296 @Property_RO 

297 def _b2c2(self): 

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

299 ''' 

300 _, b, c = self._abc3 

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

302 

303 @Property_RO 

304 def c(self): 

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

306 ''' 

307 _, _, c = self._abc3 

308 return c 

309 

310 @Property_RO 

311 def _c2_b2(self): 

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

313 ''' 

314 _, b, c = self._abc3 

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

316 

317 @Property_RO 

318 def e2ab(self): 

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

320 ''' 

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

322 

323 @Property_RO 

324 def _1e2ab(self): 

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

326 ''' 

327 a, b, _ = self._abc3 

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

329 

330 @Property_RO 

331 def e2ac(self): 

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

333 ''' 

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

335 

336 @Property_RO 

337 def _1e2ac(self): 

338 '''(INTERNAL) Get C{1 - e2ac} == 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 e2bc(self): 

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

346 ''' 

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

348 

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

350 

351 @Property_RO 

352 def _Elliptic(self): 

353 '''(INTERNAL) Get class L{Elliptic} once. 

354 ''' 

355 return _MODS.elliptic.Elliptic 

356 

357 def hartzell4(self, pov, los=None, name=NN): 

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

359 from a Point-Of-View in space. 

360 

361 @see: Function L{pygeodesy.hartzell4} for further details. 

362 ''' 

363 return hartzell4(pov, los=los, tri_biax=self, name=name) 

364 

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

366 '''Compute the projection on and the height of a cartesian above or below 

367 this triaxial's surface. 

368 

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

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

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

372 @kwarg z: Z component (C{scalar}), required if B{C{x_xyz}} if C{scalar}. 

373 @kwarg normal: If C{True} the projection is perpendicular to (the nearest 

374 point on) this triaxial's surface, otherwise the C{radial} 

375 line to this triaxial's center (C{bool}). 

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

377 negative value to skip validation. 

378 

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

380 C{x}, C{y} and C{z} of the projection on or the intersection 

381 with and with the height C{h} above or below the triaxial's 

382 surface in C{meter}, conventionally. 

383 

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

385 convergence in root finding or validation failed. 

386 

387 @see: Method L{Ellipsoid.height4} and I{Eberly}'s U{Distance from a Point 

388 to ... an Ellipsoid ...<https://www.GeometricTools.com/Documentation/ 

389 DistancePointEllipseEllipsoid.pdf>}. 

390 ''' 

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

392 

393 i, h = None, v.length 

394 if h < EPS0: # EPS 

395 x = y = z = _0_0 

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

397 elif r: # .isSpherical 

398 x, y, z = v.times(r / h).xyz 

399 h -= r 

400 else: 

401 x, y, z = v.xyz 

402 try: 

403 if normal: # perpendicular to triaxial 

404 x, y, z, h, i = _normalTo5(x, y, z, self, eps=eps) 

405 else: # radially to triaxial's center 

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

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

408 except Exception as e: 

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

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

411 h = -h # below the surface 

412 return Vector4Tuple(x, y, z, h, iteration=i, name=self.height4.__name__) 

413 

414 @Property_RO 

415 def isOrdered(self): 

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

417 ''' 

418 a, b, c = self._abc3 

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

420 

421 @Property_RO 

422 def isSpherical(self): 

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

424 ''' 

425 a, b, c = self._abc3 

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

427 

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

429 '''Get a 3-D vector perpendicular to at a cartesian on this triaxial's surface. 

430 

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

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

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

434 @kwarg z: Z component (C{scalar}), required if B{C{x_xyz}} if C{scalar}. 

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

436 

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

438 in- or outward for neg- respectively positive B{C{length}}. 

439 

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

441 surface, use method L{Triaxial.sideOf} to validate. 

442 ''' 

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

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

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

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

447 n = self._normal3d.times_(*_otherV3d_(x_xyz, y, z).xyz) 

448 if n.length < EPS0: 

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

450 return n.times(length / n.length) 

451 

452 @Property_RO 

453 def _normal3d(self): 

454 '''(INTERNAL) Get M{Vector3d((d/a)**2, (d/b)**2, (d/c)**2)}, M{d = max(a, b, c)}. 

455 ''' 

456 d = max(self._abc3) 

457 t = tuple(((d / x)**2 if x != d else _1_0) for x in self._abc3) 

458 return Vector3d(*t, name=self.normal3d.__name__) 

459 

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

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

462 ''' 

463 if fabs(_hypot21(s, c)) > EPS02: 

464 s, c = norm2(s, c) 

465 if a: 

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

467 return float0_(s, c) 

468 

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

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

471 

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

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

474 ''' 

475 ijk = self._order_ijk(**reverse) 

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

477 

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

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

480 

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

482 ''' 

483 ijk = self._order_ijk(**reverse) 

484 return v.classof(*_getitems(v.xyz, *ijk)) if ijk else v 

485 

486 @Property_RO 

487 def _ordered4(self): 

488 '''(INTERNAL) Helper for C{_hartzell3d2} and C{_normalTo5}. 

489 ''' 

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

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

492 

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

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

495 ''' 

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

497 if a < b: 

498 a, b, i, j = b, a, j, i 

499 if a < c: 

500 a, c, i, k = c, a, k, i 

501 if b < c: 

502 b, c, j, k = c, b, k, j 

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

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

505 return (a, b, c), ijk 

506 

507 abc, T = self._abc3, self 

508 if not self.isOrdered: 

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

510 if ijk: 

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

512 T = Triaxial_(*abc) 

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

514 return abc + (T,) 

515 

516 def _order_ijk(self, reverse=False): 

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

518 ''' 

519 return self._kji if reverse else self._ijk 

520 

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

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

523 ''' 

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

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

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

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

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

529 

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

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

532 

533 a, b, c = self._abc3 

534 if a != b: 

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

536 if a != c: 

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

538 z, r = c * sa, c * ca 

539 x, y = r * cb, r * sb 

540 return x, y, z 

541 

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

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

544 

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

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

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

548 @kwarg z: Z component (C{scalar}), required if B{C{x_xyz}} if C{scalar}. 

549 @kwarg eps: Near surface tolerance(C{scalar}). 

550 

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

552 within tolerance B{C{eps}}, otherwise a neg- or positive C{float} 

553 if in- respectively outside this triaxial. 

554 

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

556 ''' 

557 return _sideOf(_otherV3d_(x_xyz, y, z).xyz, self._abc3, eps=eps) 

558 

559 def _sqrt(self, x): 

560 '''(INTERNAL) Helper, see L{pygeodesy.sqrt0}. 

561 ''' 

562 if x < 0: 

563 raise TriaxialError(Fmt.PAREN(sqrt=x)) 

564 return _0_0 if x < EPS02 else sqrt(x) 

565 

566 def toEllipsoid(self, name=NN): 

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

568 

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

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

571 

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

573 

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

575 ''' 

576 a, b, c = self._abc3 

577 if a == b: 

578 b = c # N = c-Z 

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

580 a, b = b, a 

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

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

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

584 return Ellipsoid(a, b=b, name=name or self.name) 

585 

586 def toStr(self, prec=9, name=NN, **unused): # PYCHOK signature 

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

588 

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

590 @kwarg name: Override name (C{str}) or C{None} to exclude 

591 this triaxial's name. 

592 

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

594 ''' 

595 T = Triaxial_ 

596 t = T.a, 

597 J = JacobiConformalSpherical 

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

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

600 J = JacobiConformal 

601 if isinstance(self, J): 

602 t += J.xyQ2, 

603 t += T.volume, T.area 

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

605 

606 @Property_RO 

607 def volume(self): 

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

609 ''' 

610 a, b, c = self._abc3 

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

612 

613 

614class Triaxial(Triaxial_): 

615 '''I{Ordered} triaxial ellipsoid. 

616 

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

618 ''' 

619 _unordered = False 

620 

621 def __init__(self, a_triaxial, b=None, c=None, name=NN): 

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

623 

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

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

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

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

628 @kwarg c: Smallest semi-axis (C{meter}, same units as B{C{a}}), required 

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

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

631 

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

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

634 

635 @raise TriaxialError: Semi-axes not ordered, spherical or invalid. 

636 ''' 

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

638 

639 @Property_RO 

640 def _a2b2_a2c2(self): 

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

642 ''' 

643 return self._a2b2 / self._a2c2 

644 

645 @Property_RO 

646 def area(self): 

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

648 

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

650 ''' 

651 a, b, c = self._abc3 

652 if a != b: 

653 kp2, k2 = self._k2_kp2 # swapped! 

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

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

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

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

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

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

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

661 else: # a == b > c 

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

663 return a 

664 

665 def _exyz3(self, u): 

666 '''(INTERNAL) Helper for C{.forwardBetOmg}. 

667 ''' 

668 if u > 0: 

669 u2 = u**2 

670 x = u * self._sqrt(_1_0 + self._a2c2 / u2) 

671 y = u * self._sqrt(_1_0 + self._b2c2 / u2) 

672 else: 

673 x = y = u = _0_0 

674 return x, y, u 

675 

676 def forwardBetaOmega(self, beta, omega, height=0, name=NN): 

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

678 and height to cartesian. 

679 

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

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

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

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

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

685 

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

687 

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

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

690 ''' 

691 if height: 

692 h = Height_(height=height, low=-self.c, Error=TriaxialError) 

693 x, y, z = self._exyz3(h + self.c) 

694 else: 

695 x, y, z = self._abc3 # == self._exyz3(self.c) 

696 if z: # and x and y: 

697 sa, ca = SinCos2(beta) 

698 sb, cb = SinCos2(omega) 

699 

700 r = self._a2b2_a2c2 

701 x *= cb * self._sqrt(ca**2 + r * sa**2) 

702 y *= ca * sb 

703 z *= sa * self._sqrt(_1_0 - r * cb**2) 

704 return Vector3Tuple(x, y, z, name=name) 

705 

706 def forwardBetaOmega_(self, sbeta, cbeta, somega, comega, name=NN): 

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

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

709 

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

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

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

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

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

715 

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

717 

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

719 

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

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

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

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

724 ''' 

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

726 return Vector3Tuple(*t, name=name) 

727 

728 def forwardCartesian(self, x_xyz, y=None, z=None, name=NN, **normal_eps): 

729 '''Project a cartesian on this triaxial. 

730 

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

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

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

734 @kwarg z: Z component (C{scalar}), required if B{C{x_xyz}} if C{scalar}. 

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

736 @kwarg normal_eps: Optional keyword arguments C{B{normal}=True} and 

737 C{B{eps}=EPS}, see method L{Triaxial.height4}. 

738 

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

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

741 ''' 

742 t = self.height4(x_xyz, y, z, **normal_eps) 

743 _ = t.rename(name) 

744 return t 

745 

746 def forwardLatLon(self, lat, lon, height=0, name=NN): 

747 '''Convert I{geodetic} lat-, longitude and heigth to cartesian. 

748 

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

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

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

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

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

754 

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

756 

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

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

759 ''' 

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

761 

762 def forwardLatLon_(self, slat, clat, slon, clon, height=0, name=NN): 

763 '''Convert I{geodetic} lat-, longitude and heigth to cartesian. 

764 

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

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

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

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

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

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

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

772 

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

774 

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

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

777 ''' 

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

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

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

781 

782 def _forwardLatLon3(self, h, name, sa, ca, sb, cb): 

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

784 ''' 

785 ca_x_sb = ca * sb 

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

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

788 n = self.a / self._sqrt(t) # prime vertical 

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

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

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

792 return Vector3Tuple(x, y, z, name=name) 

793 

794 @Property_RO 

795 def _k2_kp2(self): 

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

797 ''' 

798 # k2 = a2b2 / a2c2 * c2_b2 

799 # kp2 = b2c2 / a2c2 * a2_b2 

800 # b2 = b**2 

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

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

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

804 return (self._a2b2_a2c2 * self._c2_b2, 

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

806 

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

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

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

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

811 ''' 

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

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

814 

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

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

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

818 

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

820 z2 = Fsum(c2_a2, sb**2, b2_a2 * cb**2).fover(a2c2_a2) 

821 

822 x, y, z = self._abc3 

823 x *= cb * self._sqrt(x2) 

824 y *= ca * sb 

825 z *= sa * self._sqrt(z2) 

826 return x, y, z 

827 

828 def reverseBetaOmega(self, x_xyz, y=None, z=None, name=NN): 

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

830 and height. 

831 

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

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

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

835 @kwarg z: Z component (C{scalar}), required if B{C{x_xyz}} if C{scalar}. 

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

837 

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

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

840 units as this triaxial's axes. 

841 

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

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

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

845 ''' 

846 v = _otherV3d_(x_xyz, y, z) 

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

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

849 

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

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

852 

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

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

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

856 @kwarg z: Z component (C{scalar}), required if B{C{x_xyz}} if C{scalar}. 

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

858 as the axes). 

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

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

861 @kwarg eps: Tolerance for surface test (C{scalar}). 

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

863 

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

865 

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

867 

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

869 ''' 

870 v = _otherV3d_(x_xyz, y, z, name=name) 

871 s = _sideOf(v.xyz, self._abc3, eps=eps) 

872 if s: # PYCHOK no cover 

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

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

875 

876 if h: 

877 if normal: 

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

879 elif v.length > EPS0: 

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

881 return v.xyz # Vector3Tuple 

882 

883 def reverseLatLon(self, x_xyz, y=None, z=None, name=NN): 

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

885 

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

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

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

889 @kwarg z: Z component (C{scalar}), required if B{C{x_xyz}} if C{scalar}. 

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

891 

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

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

894 as this triaxial's axes. 

895 

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

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

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

899 ''' 

900 v = _otherV3d_(x_xyz, y, z) 

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

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

903 _1_0) 

904 t = self._reverseLatLon3(s, atan2d, v, self.forwardLatLon_) 

905 return LatLon3Tuple(*t, name=name) 

906 

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

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

909 ''' 

910 x, y, z = s.xyz 

911 d = hypot( x, y) 

912 a = atan2_(z, d) 

913 b = atan2_(y, x) 

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

915 return a, b, h 

916 

917 

918class JacobiConformal(Triaxial): 

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

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

921 

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

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

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

925 returned in the case of an ellipsoid of revolution. 

926 

927 Copyright (C) U{Charles Karney<mailto:Charles@Karney.com>} (2014-2023) and 

928 licensed under the MIT/X11 License. 

929 

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

931 L{JacobiConformalSpherical}. 

932 

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

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

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

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

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

938 ''' 

939 

940 @Property_RO 

941 def _xE(self): 

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

943 ''' 

944 k2, kp2 = self._k2_kp2 

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

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

947 

948 def xR(self, omega): 

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

950 

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

952 

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

954 ''' 

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

956 

957 def xR_(self, somega, comega): 

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

959 

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

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

962 

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

964 ''' 

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

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

967 

968 @Property_RO 

969 def xyQ2(self): 

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

971 ''' 

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

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

974 name=JacobiConformal.xyQ2.name) 

975 

976 def xyR2(self, beta, omega, name=NN): 

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

978 

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

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

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

982 

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

984 ''' 

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

986 name=name or self.xyR2.__name__) 

987 

988 def xyR2_(self, sbeta, cbeta, somega, comega, name=NN): 

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

990 

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

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

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

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

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

996 

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

998 ''' 

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

1000 self.yR_(sbeta, cbeta), 

1001 name=name or self.xyR2_.__name__) 

1002 

1003 @Property_RO 

1004 def _yE(self): 

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

1006 ''' 

1007 kp2, k2 = self._k2_kp2 # swapped! 

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

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

1010 

1011 def yR(self, beta): 

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

1013 

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

1015 

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

1017 ''' 

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

1019 

1020 def yR_(self, sbeta, cbeta): 

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

1022 

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

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

1025 

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

1027 ''' 

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

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

1030 

1031 

1032class JacobiConformalSpherical(JacobiConformal): 

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

1034 

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

1036 ''' 

1037 _ab = _bc = 0 

1038 

1039 def __init__(self, radius_triaxial, ab=0, bc=0, name=NN): 

1040 '''New L{JacobiConformalSpherical}. 

1041 

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

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

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

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

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

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

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

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

1050 

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

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

1053 + B{bc})} not positive. 

1054 

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

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

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

1058 and C{c} are copied. 

1059 ''' 

1060 try: 

1061 r, j = radius_triaxial, False 

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

1063 if (not (ab or bc)) and isinstance(r, JacobiConformalSpherical): 

1064 j = True 

1065 t = r._abc3 

1066 else: 

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

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

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

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

1071 raise ValueError('(ab + bc)') 

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

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

1074 and isfinite(self._a2c2)): 

1075 raise ValueError(_not_(_finite_)) 

1076 except (TypeError, ValueError) as x: 

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

1078 if name: 

1079 self.name = name 

1080 

1081 @Property_RO 

1082 def ab(self): 

1083 '''Get relative magnitude C{ab} (C{meter}, same units as B{C{a}}). 

1084 ''' 

1085 return self._ab 

1086 

1087 @Property_RO 

1088 def _a2b2(self): 

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

1090 ''' 

1091 a, b, _ = self._abc3 

1092 return self.ab * (a + b) 

1093 

1094 @Property_RO 

1095 def _a2c2(self): 

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

1097 ''' 

1098 return self._a2b2 + self._b2c2 

1099 

1100 @Property_RO 

1101 def bc(self): 

1102 '''Get relative magnitude C{bc} (C{meter}, same units as B{C{a}}). 

1103 ''' 

1104 return self._bc 

1105 

1106 @Property_RO 

1107 def _b2c2(self): 

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

1109 ''' 

1110 _, b, c = self._abc3 

1111 return self.bc * (b + c) 

1112 

1113 @Property_RO 

1114 def radius(self): 

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

1116 ''' 

1117 return self.a 

1118 

1119 

1120class TriaxialError(_ValueError): 

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

1122 ''' 

1123 pass # ... 

1124 

1125 

1126class Triaxials(_NamedEnum): 

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

1128 to accommodate the L{_LazyNamedEnumItem} properties. 

1129 ''' 

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

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

1132 ''' 

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

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

1135 

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

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

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

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

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

1141_E = _WGS84.ellipsoid 

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

1143 Amalthea = _lazy('Amalthea', 125.0, 73.0, 64), # Jupiter 

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

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

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

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

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

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

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

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

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

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

1154 WGS84_35 = _lazy('WGS84_35', *map1(m2km, _E.a + 35, _E.a - 35, _E.b))) 

1155del _E 

1156 

1157 

1158def _getitems(items, *indices): 

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

1160 

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

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

1163 the special method C{__getitem__}. 

1164 ''' 

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

1166 

1167 

1168def _hartzell3d2(pov, los, Tun): # MCCABE 13 in .ellipsoidal.hartzell4, .formy.hartzell 

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

1170 formula for I{un-/ordered} triaxials. 

1171 ''' 

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

1173 

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

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

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

1177 

1178 p3 = T._order3d(_otherV3d(pov=pov)) 

1179 u3 = T._order3d(_otherV3d(los=los)) if los else p3.negate() 

1180 u3 = u3.unit() # unit vector, opposing signs 

1181 

1182 x2, y2, z2 = p3.x2y2z2 # p3.times_(p3).xyz 

1183 ux, vy, wz = u3.times_(p3).xyz 

1184 u2, v2, w2 = u3.x2y2z2 # u3.times_(u3).xyz 

1185 

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

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

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

1189 raise _ValueError(_near_(_null_)) 

1190 

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

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

1193 -w2 * x2 * p2, -u2 * y2 * q2, -v2 * x2 * q2, ux * vy * 2 * q2) 

1194 if r > 0: # a2 factored out 

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

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

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

1198 

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

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

1201 s = fsumf_(_1_0, x2 / a2, y2 / b2, z2 / c2, _N_2_0) # like _sideOf 

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

1203 elif fsum1f_(x2, y2, z2) < d**2: # d past triaxial's center 

1204 raise _ValueError(_too_(_distant_)) 

1205 

1206 v = p3.minus(u3.times(d)) # Vector3d 

1207 h = p3.minus(v).length # distance to triaxial 

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

1209 

1210 

1211def hartzell4(pov, los=None, tri_biax=_WGS84, name=NN): 

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

1213 from a Point-Of-View outside. 

1214 

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

1216 or L{Vector3d}). 

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

1218 C{None} to point to the tri-/biaxial's center. 

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

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

1221 L{Ellipsoid}, L{Ellipsoid2}, L{a_f2Tuple} or C{scalar} radius, 

1222 conventionally in C{meter}). 

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

1224 

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

1226 C{h} the distance from B{C{pov}} to C{(x, y, z)} along the B{C{los}}, 

1227 all in C{meter}, conventionally. 

1228 

1229 @raise TriaxialError: Null B{C{pov}} or B{C{los}}, or B{C{pov}} is inside the 

1230 tri-/biaxial or B{C{los}} points outside the tri-/biaxial 

1231 or points in an opposite direction. 

1232 

1233 @raise TypeError: Invalid B{C{pov}} or B{C{los}}. 

1234 

1235 @see: Function L{pygeodesy.hartzell}, L{pygeodesy.tyr3d} for B{C{los}} and 

1236 U{I{Satellite Line-of-Sight Intersection with Earth}<https://StephenHartzell. 

1237 Medium.com/satellite-line-of-sight-intersection-with-earth-d786b4a6a9b6>}. 

1238 ''' 

1239 if isinstance(tri_biax, Triaxial_): 

1240 T = tri_biax 

1241 else: 

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

1243 _spherical_datum(tri_biax, name=hartzell4.__name__) 

1244 T = D.ellipsoid._triaxial 

1245 

1246 try: 

1247 v, h = _hartzell3d2(pov, los, T) 

1248 except Exception as x: 

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

1250 return Vector4Tuple(v.x, v.y, v.z, h, name=name or hartzell4.__name__) 

1251 

1252 

1253def _hypot21(x, y, z=0): 

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

1255 fabs(y), fabs(z))} rarely greater than 1.0. 

1256 ''' 

1257 return fsumf_(_1_0, x**2, y**2, z**2, _N_2_0) if z else \ 

1258 fsumf_(_1_0, x**2, y**2, _N_2_0) 

1259 

1260 

1261def _normalTo4(x, y, a, b, eps=EPS): 

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

1263 

1264 @see: Function C{pygeodesy.ellipsoids._normalTo3} and I{Eberly}'s U{Distance 

1265 from a Point to ... an Ellipsoid ...<https://www.GeometricTools.com/ 

1266 Documentation/DistancePointEllipseEllipsoid.pdf>}. 

1267 ''' 

1268 if a < b: 

1269 b, a, d, i = _normalTo4(y, x, b, a, eps=eps) 

1270 return a, b, d, i 

1271 

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

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

1274 

1275 i = None 

1276 if y: 

1277 if x: 

1278 u = fabs(x / a) 

1279 v = fabs(y / b) 

1280 g = _hypot21(u, v) 

1281 if g: 

1282 r = (a / b)**2 

1283 t, i = _rootXd(r, 0, u, 0, v, g, eps) 

1284 a = x / (t / r + _1_0) 

1285 b = y / (t + _1_0) 

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

1287 else: # on the ellipse 

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

1289 else: # x == 0 

1290 if y < 0: 

1291 b = -b 

1292 a, d = x, fabs(y - b) 

1293 

1294 else: # y == 0 

1295 n = a * x 

1296 d = (a + b) * (a - b) 

1297 if d > fabs(n): # PYCHOK no cover 

1298 r = n / d 

1299 a *= r 

1300 b *= sqrt(_1_0 - r**2) 

1301 d = hypot(x - a, b) 

1302 else: 

1303 if x < 0: 

1304 a = -a 

1305 b, d = y, fabs(x - a) 

1306 return a, b, d, i 

1307 

1308 

1309def _normalTo5(x, y, z, Tun, eps=EPS): # MCCABE 19 

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

1311 

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

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

1314 ''' 

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

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

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

1318 a, b, c, d, i = _normalTo5(*t, eps=eps) 

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

1320 

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

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

1323 

1324 if eps > 0: 

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

1326 else: # no validation 

1327 val, eps = 0, -eps 

1328 

1329 i = None 

1330 if z: 

1331 if y: 

1332 if x: 

1333 u = fabs(x / a) 

1334 v = fabs(y / b) 

1335 w = fabs(z / c) 

1336 g = _hypot21(u, v, w) 

1337 if g: 

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

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

1340 t, i = _rootXd(_1_0 / r, _1_0 / s, u, v, w, g, eps) 

1341 a = x / (t * r + _1_0) 

1342 b = y / (t * s + _1_0) 

1343 c = z / (t + _1_0) 

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

1345 else: # on the ellipsoid 

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

1347 else: # x == 0 

1348 a = x # 0 

1349 b, c, d, i = _normalTo4(y, z, b, c, eps=eps) 

1350 elif x: # y == 0 

1351 b = y # 0 

1352 a, c, d, i = _normalTo4(x, z, a, c, eps=eps) 

1353 else: # x == y == 0 

1354 if z < 0: 

1355 c = -c 

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

1357 

1358 else: # z == 0 

1359 t = False 

1360 n = a * x 

1361 d = T._a2c2 # (a + c) * (a - c) 

1362 if d > fabs(n): 

1363 u = n / d 

1364 n = b * y 

1365 d = T._b2c2 # (b + c) * (b - c) 

1366 if d > fabs(n): 

1367 v = n / d 

1368 n = _hypot21(u, v) 

1369 if n < 0: 

1370 a *= u 

1371 b *= v 

1372 c *= sqrt(-n) 

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

1374 t = True 

1375 if not t: 

1376 c = z # 0 

1377 a, b, d, i = _normalTo4(x, y, a, b, eps=eps) 

1378 

1379 if val > 0: # validate 

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

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

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

1383 sideOf=e, eps=val) 

1384 if d: # angle of delta and normal vector 

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

1386 if m.euclid > val: 

1387 m = m.unit() 

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

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

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

1391 raise _ValueError(n=n, m=m, 

1392 dot=e, eps=val) 

1393 return a, b, c, d, i 

1394 

1395 

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

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

1398 ''' 

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

1400 _otherV3d(x_xyz=x_xyz) 

1401 

1402 

1403def _rootXd(r, s, u, v, w, g, eps): 

1404 '''(INTERNAL) Robust 2d- or 3d-root finder: 

1405 2d- if C{s == v == 0} otherwise 3d-root. 

1406 ''' 

1407 _1, __2 = _1_0, _0_5 

1408 _a, _h2 = fabs, _hypot21 

1409 

1410 u *= r 

1411 v *= s # 0 for 2d-root 

1412 t0 = w - _1 

1413 t1 = _0_0 if g < 0 else _h2(u, w, v) 

1414 for i in range(1, _TRIPS): 

1415 e = _a(t0 - t1) 

1416 if e < eps: 

1417 break 

1418 t = (t0 + t1) * __2 

1419 if t in (t0, t1): 

1420 break 

1421 g = _h2(u / (t + r), w / (t + _1), 

1422 (v / (t + s)) if v else 0) 

1423 if g > 0: 

1424 t0 = t 

1425 elif g < 0: 

1426 t1 = t 

1427 else: 

1428 break 

1429 else: # PYCHOK no cover 

1430 t = Fmt.no_convergence(e, eps) 

1431 raise _ValueError(t, txt=_rootXd.__name__) 

1432 return t, i 

1433 

1434 

1435def _sideOf(xyz, abc, eps=EPS): # in .formy 

1436 '''(INTERNAL) Helper for C{_hartzell3d2}, M{.sideOf} and M{.reverseCartesian}. 

1437 

1438 @return: M{sum((x / a)**2 for x, a in zip(xyz, abc)) - 1} or C{INT0}, 

1439 ''' 

1440 s = _hypot21(*((x / a) for x, a in _zip(xyz, abc) if a)) # strict=True 

1441 return s if fabs(s) > eps else INT0 

1442 

1443 

1444if __name__ == '__main__': 

1445 

1446 from pygeodesy import printf 

1447 from pygeodesy.interns import _COMMA_, _NL_, _NLATvar_ 

1448 

1449 # __doc__ of this file, force all into registery 

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

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

1452 

1453# **) MIT License 

1454# 

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

1456# 

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

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

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

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

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

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

1463# 

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

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

1466# 

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

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

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

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

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

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

1473# OTHER DEALINGS IN THE SOFTWARE.