Coverage for pygeodesy / triaxials / triaxial3.py: 90%

1 

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

3 

4u'''I{Ordered} triaxial ellipsoid classes L{Triaxial3} and L{Triaxial3B} for conversion between 

5variuos lat-/longitudal and cartesian coordinates on a triaxial ellipsoid using L{Ang}, L{Deg}, 

6L{Rad} lat-, longitude and heading angles. 

7 

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

9GeographicLib.SourceForge.io/C++/doc/classGeographicLib_1_1Triaxial_1_1Ellipsoidal3.html>} and U{Cartesian3 

10<https://GeographicLib.SourceForge.io/C++/doc/classGeographicLib_1_1Triaxial_1_1Cartesian3.html>}. 

11 

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

13License. For more information, see the U{GeographicLib 2.7 <https://GeographicLib.SourceForge.io/>} 

14documentation. 

15 

16@var Triaxial3s.Amalthea: Triaxial3(name='Amalthea', a=125000, b=73000, c=64000, k2=0.106947697, kp2=0.893052303, volume=2446253479595252, area=93239507787.490356445, R2=86138.05359954) 

17@var Triaxial3s.Ariel: Triaxial3(name='Ariel', a=581100, b=577900, c=577700, k2=0.05866109, kp2=0.94133891, volume=812633172614203904, area=4211301462766.580078125, R2=578899.578791275) 

18@var Triaxial3s.Earth: Triaxial3(name='Earth', a=6378173.435, b=6378103.9, c=6356754.399999999, k2=0.996748146, kp2=0.003251854, volume=1083208241574987694080, area=510065911057440.9375, R2=6371008.987886564) 

19@var Triaxial3s.Enceladus: Triaxial3(name='Enceladus', a=256600, b=251400, c=248300, k2=0.369647336, kp2=0.630352664, volume=67094551514082248, area=798618496278.596679688, R2=252095.300756832) 

20@var Triaxial3s.Europa: Triaxial3(name='Europa', a=1564130, b=1561230, c=1560930, k2=0.093663002, kp2=0.906336998, volume=15966575194402123776, area=30663773697323.51953125, R2=1562096.533153486) 

21@var Triaxial3s.Io: Triaxial3(name='Io', a=1829400, b=1819300, c=1815700, k2=0.262045618, kp2=0.737954382, volume=25313121117889765376, area=41691875849096.734375, R2=1821464.812747882) 

22@var Triaxial3s.Mars: Triaxial3(name='Mars', a=3394600, b=3393300, c=3376300, k2=0.92878339, kp2=0.07121661, volume=162907283585817247744, area=144249140795107.4375, R2=3388064.624110653) 

23@var Triaxial3s.Mimas: Triaxial3(name='Mimas', a=207400, b=196800, c=190600, k2=0.359218713, kp2=0.640781287, volume=32587072869017956, area=493855762247.691833496, R2=198241.75359411) 

24@var Triaxial3s.Miranda: Triaxial3(name='Miranda', a=240400, b=234200, c=232900, k2=0.171062751, kp2=0.828937249, volume=54926187094835456, area=698880863325.757080078, R2=235828.692095158) 

25@var Triaxial3s.Moon: Triaxial3(name='Moon', a=1735550, b=1735324, c=1734898, k2=0.653331685, kp2=0.346668315, volume=21886698675223740416, area=37838824729886.09375, R2=1735257.329122863) 

26@var Triaxial3s.Tethys: Triaxial3(name='Tethys', a=535600, b=528200, c=525800, k2=0.243190549, kp2=0.756809451, volume=623086233855821440, area=3528073490771.393554688, R2=529863.348254881) 

27@var Triaxial3s.WGS84_3: Triaxial3(name='WGS84_3', a=6378171.36, b=6378101.609999999, c=6356751.84, k2=0.996738165, kp2=0.003261835, volume=1083207064030173855744, area=510065541435967.5, R2=6371006.679496506) 

28@var Triaxial3s.WGS84_35: Triaxial3(name='WGS84_35', a=6378172, b=6378102, c=6356752.314245179, k2=0.996726499, kp2=0.003273501, volume=1083207319768789942272, area=510065621722018.25, R2=6371007.180905545) 

29@var Triaxial3s.WGS84_3r: Triaxial3(name='WGS84_3r', a=6378172, b=6378102, c=6356752, k2=0.996726547, kp2=0.003273453, volume=1083207266220584468480, area=510065604942135.875, R2=6371007.076110449) 

30''' 

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

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

33 

34from pygeodesy.angles import Ang, Ang_, _Ang3Tuple, atan2, sincos2, _SinCos2 

35from pygeodesy.basics import _copysign, map1 

36from pygeodesy.constants import EPS, EPS_2, EPS02, EPS8, INT0, NAN, \ 

37 _EPSqrt, _SQRT3, _copysign_0_0, _copysign_1_0, \ 

38 _flipsign, _isfinite, _over, _1_over, _0_0, \ 

39 _0_5, _N_1_0, _1_0, _2_0, _3_0, _4_0, _9_0 

40from pygeodesy.errors import _xattr, _xkwds, _xkwds_get, _xkwds_pop2 

41from pygeodesy.fmath import cbrt2, fdot, hypot, hypot2, norm2, fabs, sqrt 

42from pygeodesy.fsums import Fsum, fsumf_, Fmt 

43from pygeodesy.interns import NN, _azimuth_, _DMAIN_, _h_, _lam_, _lat_, \ 

44 _lon_, _phi_ 

45# from pygeodesy.lazily import _ALL_LAZY # from .vector3d 

46# from pygeodesy.named import _Pass # from .namedTuples 

47from pygeodesy.namedTuples import _NamedTuple, _xyzh_Tuple, _Pass, Property_RO 

48# from pygeodesy.props import Property_RO # from .namedTuples 

49# from pygeodesy.streprs import Fmt # from .fsums 

50from pygeodesy.triaxials.bases import _bet_, _HeightINT0, LLK, _llk_, \ 

51 _MAXIT, _omg_, _otherV3d_, _sqrt0, \ 

52 _Triaxial3Base, TriaxialError, \ 

53 _TriaxialsBase 

54from pygeodesy.units import Degrees, Lat, Lon, Meter, Radians, Radius_ 

55# from pygeodesy.utily import atan2, sincos2 # from .triaxials.angles 

56from pygeodesy.vector3d import Vector3d, _ALL_LAZY 

57 

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

59from random import random 

60 

61__all__ = _ALL_LAZY.triaxials_triaxial3 

62__version__ = '26.08.06' 

63 

64_alp_ = 'alp' 

65_NAN3d = Vector3d(NAN, NAN, NAN) 

66_TOL = cbrt2(EPS) 

67_TOL2 = _TOL**2 # cbrt(EPS)**4 

68_zet_ = 'zet' 

69_27_0 = 27.0 

70 

71 

72class BetOmgAlp5Tuple(_Ang3Tuple): 

73 '''5-Tuple C{(bet, omg, alp, h, llk)} with I{ellipsoidal} 

74 lat- C{bet}, longitude C{omg} and azimuth C{alp}, all 

75 in L{Ang}les on and height C{h} off the triaxial's 

76 surface and kind C{llk} set to C{LLK.ELLIPSOIDAL}. 

77 ''' 

78 _Names_ = (_bet_, _omg_, _alp_, _h_, _llk_) 

79 _Units_ = ( Ang, Ang, _Pass, _HeightINT0, _Pass) 

80 

81 

82class Cartesian5Tuple(_xyzh_Tuple): 

83 '''5-Tuple C{(x, y, z, h, llk)} with I{cartesian} C{x}, C{y} 

84 and C{z} coordinates on and height C{h} above or below the 

85 triaxial's surface and kind C{llk} set to the original 

86 C{LLK} or C{None}. 

87 ''' 

88 _Names_ = _xyzh_Tuple._Names_ + (_llk_,) 

89 _Units_ = (Meter, Meter, Meter, _HeightINT0, _Pass) 

90 

91 def __new__(cls, x, y, z, h=0, llk=None, **kwds): # **iteration_name 

92 args = x, y, z, (h or INT0), llk 

93 return _xyzh_Tuple.__new__(cls, args, **kwds) 

94 

95 

96class _Fp2(object): 

97 '''(INTERNAL) Function and derivate evaluation. 

98 ''' 

99 def __init__(self, rs, ls, n=1): 

100 # assert 0 < n <= 2 

101 self._2 = n == 2 

102 self._rls = tuple((p, q) for p, q in zip(rs, ls) if p) 

103 

104 def __call__(self, p): 

105 # Evaluate C{f(p) = sum((rs[k] / (p + ls[k]))**n, 

106 # k=0..2) - 1} and its derivative C{fp}. 

107 f = _N_1_0 

108 fc = fp = _0_0 

109 _D = EPS_2 

110 _2 = self._2 

111 for g, q in self._rls: 

112 q = _1_over(p + q) 

113 g *= q 

114 if _2: 

115 g *= g 

116 q += q 

117 r = round(g / _D) * _D 

118 f += r 

119 fc += g - r 

120 fp -= g * q 

121 return (f + fc), fp 

122 

123 

124class LatLonAzi5Tuple(_NamedTuple): 

125 '''5-Tuple C{(lat, lon, azimuth, h, llk)} with triaxial 

126 C{lat}-, C{lon}gitude and C{azimuth}, all in C{degrees}, 

127 height C{h} off the triaxial's surface and kind C{llk} 

128 either C{LLK.GEODETIC} or C{LLK.GEODESIC_LON0}. 

129 ''' 

130 _Names_ = (_lat_, _lon_, _azimuth_, _h_, _llk_) 

131 _Units_ = ( Lat, Lon, _Pass, _HeightINT0, _Pass) 

132 

133 

134class PhiLamZet5Tuple(_Ang3Tuple): 

135 '''5-Tuple C{(phi, lam, zet, h, llk)} with triaxial lat- 

136 lat- C{phi}, longitude C{lam} and azimuth C{zet}, all 

137 in L{Ang}les on and height C{h} off the triaxial's 

138 surface and kind C{llk} set to an C{LLK}. 

139 ''' 

140 _Names_ = (_phi_, _lam_, _zet_, _h_, _llk_) 

141 _Units_ = ( Ang, Ang, _Pass, _HeightINT0, _Pass) 

142 

143 def toLatLonAzi5(self, **name): 

144 '''Return this tuple as L{LatLonAzi5Tuple} with 

145 C{lat}, C{lon} and C{azimuth} in C{degrees} or 

146 C{azimuth} C{INT0}. 

147 ''' 

148 u = self # .toUnits(Error=TriaxialError) 

149 n = _xkwds_get(name, name=u.name) 

150 z = u.zet # None or isAng with turns 

151 z = INT0 if z is None else z.degrees 

152 t = u.phi.degrees0, u.lam.degrees0, z, u.h, u.llk 

153 return LatLonAzi5Tuple(t, name=n or NN) 

154 

155 

156class Triaxial3(_Triaxial3Base): 

157 '''I{Ordered} triaxial ellipsoid convering between cartesian and 

158 lat-/longitudes using using class L{Ang}. 

159 

160 @see: L{Triaxial<triaxials.triaxial5.Triaxial>} for details. 

161 ''' 

162 def _cardinal2(self, v, mer, llk): # cardinaldir 

163 '''(INTERNAL) Get 2-tuple C{(n, e)} at C{mer}idian. 

164 ''' 

165 # assert isinstance(v, Vector3d) and isinstance(mer, Ang) \ 

166 # and isinstance(llk, LLK.__class__) 

167 a2, b2, c2 = self._a2b2c23 

168 if llk._X: 

169 a2, c2 = c2, a2 

170 v = v._roty(True) # +1 

171 x, y, z = v.xyz3 

172 if x or y: 

173 s = (-z) / c2 

174 z = x**2 / a2 + y**2 / b2 

175 else: 

176 y, x, _ = mer.scn3 

177 s = _copysign_1_0(-z) 

178 z = _0_0 

179 n = Vector3d(x * s, y * s, z).unit() 

180 e = v.dividedBy_(a2, b2, c2).unit() # normvec 

181 e = n.cross(e).unit() 

182 if llk._X: 

183 e = e._roty(False) # -1 

184 n = n._roty(False) # -1 

185 return n, e 

186 

187 def forward(self, lat, lon, **height_unit_name): 

188 '''Convert a I{geodetic} lat- and longitude to a cartesian 

189 on this triaxial's surface. 

190 

191 @arg lat: Geodetic latitude (C{degrees} or B{C{unit}}). 

192 @arg lon: Geodetic longitude (C{degrees} or B{C{unit}}). 

193 @kwarg height_unit_name: Optional C{B{height}=0} (C{meter}), 

194 scalar C{B{unit}=}L{Radians} (or L{Degrees}) 

195 and C{B{name}=NN} (C{str}). 

196 

197 @return: A L{Cartesian5Tuple}C{(x, y, z, h, llk)} with 

198 C{h=B{height}} and kind C{llk=LLK.GEODETIC} or 

199 C{llk=LLK.GEODETIC_LON0}. 

200 

201 @see: Methods L{Triaxial3.forwardPhiLam} and L{Triaxial3.reverse}. 

202 ''' 

203 llk = LLK.GEODETIC_LON0 if self.Lon0 else LLK.GEODETIC 

204 height_unit_name = _xkwds(height_unit_name, unit=Degrees) 

205 return self.forwardPhiLam(lat, lon, llk=llk, **height_unit_name) 

206 

207 def forwardBetOmg(self, bet, omg, height=0, **unit_name): # elliptocart2 

208 '''Convert an I{ellipsoidal} lat- and longitude to a cartesian 

209 on this triaxial's surface. 

210 

211 @arg bet: Ellipsoidal latitude (C{Ang} or B{C{unit}}). 

212 @arg omg: Ellipsoidal longitude (C{Ang} or B{C{unit}}). 

213 @kwarg height: Height above or below this triaxial's surface (C{meter}, 

214 same units as this triaxial's semi-axes). 

215 @kwarg unit_name: Optional C{B{name}=NN} (C{str}) and scalar 

216 C{B{unit}=}L{Radians} (or L{Degrees}). 

217 

218 @return: A L{Cartesian5Tuple}C{(x, y, z, h, llk)} with C{h=B{height}} 

219 and kind C{llk=LLK.ELLIPSOIDAL}. 

220 

221 @see: Method L{Triaxial3.reverseBetOmg}. 

222 ''' 

223 ct, _ = self.forwardBetOmgAlp2(bet, omg, None, height, **unit_name) 

224 return ct 

225 

226 forwardBetaOmega = forwardBetOmg # for backward compatibility 

227 

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

229 '''DEPRECATED on 2025.11.15, like C{Triaxial.forwardBetaOmega_}.''' 

230 return self.forwardBetaOmega(Ang_(sbeta, cbeta), 

231 Ang_(somega, comega), **name) 

232 

233 def forwardBetOmgAlp2(self, bet, omg, alp, height=0, **unit_name): # elliptocart2 

234 '''Convert an I{ellipsoidal} lat-, longitude and heading to a 

235 cartesian and a direction on this triaxial's surface. 

236 

237 @arg bet: Ellipsoidal latitude (C{Ang} or B{C{unit}}). 

238 @arg omg: Ellipsoidal longitude (C{Ang} or B{C{unit}}). 

239 @arg alp: Azimuth of the heading (C{Ang}, B{C{unit}} or C{None}). 

240 @kwarg height: Height above or below this triaxial's surface (C{meter}, 

241 same units as this triaxial's semi-axes). 

242 @kwarg unit_name: Optional C{B{name}=NN} (C{str}), scalar 

243 C{B{unit}=}L{Radians} (or L{Degrees}). 

244 

245 @return: 2-Tuple C{(cartesian, direction)} with C{cartesian} a 

246 L{Cartesian5Tuple}C{(x, y, z, h, llk)} with C{h=B{height}}, 

247 kind C{llk=LLK.ELLIPSOIDAL} and C{direction} a C{Vector3d} 

248 tangent to this triaxial's surface or C{None}. 

249 

250 @see: Method L{Triaxial3.reverseBetOmgAlp}. 

251 ''' 

252 h, llk, unit, name = _h_llk_unit_name(height, **unit_name) 

253 a, b, c = self._abc3 

254 if h: # Cartesian.elliptocart 

255 a, b, _ = self._a2b2c23 

256 h = _HeightINT0(h) 

257 s = (c * _2_0 + h) * h 

258 if s < 0: 

259 s = -min(a, b, -s) 

260 a = sqrt(a + s) 

261 b = sqrt(b + s) 

262 c += h 

263 sb, cb = _SinCos2(bet, unit) 

264 so, co = _SinCos2(omg, unit) 

265 k, kp = self._k_kp 

266 tx, tz = _txtz2(cb, so, k, kp) 

267 ct = Cartesian5Tuple(a * co * tx, 

268 b * cb * so, 

269 c * sb * tz, 

270 h, llk, **name) 

271 

272 if alp is None: # or h? 

273 dir3d = None # _NAN3d? 

274 else: 

275 try: 

276 sa, ca = _SinCos2(alp, unit) 

277 except Exception as X: 

278 raise TriaxialError(alp=alp, cause=X) 

279 a, b, c = self._abc3 

280 if k and kp and not (cb or so): 

281 c2s2_b = (ca - sa) * (ca + sa) / b 

282 dir3d = Vector3d(a * k * co * c2s2_b 

283 -co * sb * ca * sa * _2_0, 

284 c * kp * sb * c2s2_b) 

285 else: 

286 if not tx: # at oblate pole tx -> |cos(bet)| 

287 c = _flipsign(co, cb) 

288 n = Vector3d(-c * sb, 

289 -so * sb, _0_0) 

290 e = Vector3d(-so, c, _0_0) 

291 elif not tz: # at prolate pole tz -> |sin(omg)| 

292 s = _flipsign(sb, so) 

293 n = Vector3d(_0_0, -s, cb) 

294 e = Vector3d(_0_0, cb * co, co * s) 

295 else: 

296 k2, kp2 = self._k2_kp2 

297 n = Vector3d(-a * k2 * sb * cb * co / tx, 

298 -b * sb * so, c * cb * tz) 

299 e = Vector3d(-a * tx * so, b * cb * co, 

300 c * kp2 * sb * so * co / tz) 

301 dir3d = n.unit().times(ca) # NAN 

302 dir3d += e.unit().times(sa) # NAN 

303 dir3d.name = ct.name 

304 return ct, dir3d 

305 

306 def forwardCartesian(self, x_ct, y=None, z=None, normal=True, **eps_llk_name): 

307 '''Project any cartesian I{onto} this triaxial's surface. 

308 

309 @arg x_ct: X component (C{scalar}) or a cartesian (L{Cartesian5Tuple} or 

310 any C{Cartesian}, L{Ecef9Tuple}, L{Vector3d}, L{Vector3Tuple} 

311 or L{Vector4Tuple}). 

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

313 ignored otherwise. 

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

315 @kwarg normal: If C{True}, the projection is C{perpendicular} to the surface, 

316 otherwise C{radial} to the center of this triaxial (C{bool}). 

317 @kwarg eps_llk_name: Root finder tolerance C{B{eps}=EPS}, kind C{B{llk}=None} 

318 overriding C{B{x_ct}.llk} and optional C{B{name}="height4"} (C{str}). 

319 

320 @return: A L{Cartesian5Tuple}C{(x, y, z, h, llk)}. 

321 

322 @see: Method L{Triaxial3.reverseCartesian} to reverse the projection and 

323 function L{height4<triaxials.triaxial5.height4>} for more details. 

324 ''' 

325 llk, kwds = _xkwds_pop2(eps_llk_name, llk=_xattr(x_ct, llk=None)) 

326 h = self.sideOf(x_ct, y, z) 

327 if h: # signed, square 

328 v = self.height4(x_ct, y, z, normal=normal, **kwds) 

329 h = v.h 

330 else: # on the surface 

331 v = _otherV3d_(x_ct, y, z) 

332 n = _xkwds_get(kwds, name=NN) 

333 return Cartesian5Tuple(v.x, v.y, v.z, h, llk, iteration=v.iteration, name=n) 

334 

335 def forwardLatLon(self, lat, lon, height=0, llk=LLK.ELLIPSOIDAL, **unit_name): # anytocart2 

336 '''Convert any lat-/longitude kind to a cartesian on this triaxial's surface. 

337 

338 @arg lat: Latitude (C{Ang} or B{C{unit}}). 

339 @arg lon: Longitude (C{Ang} or B{C{unit}}). 

340 @kwarg height: Height above or below this triaxial's surface (C{meter}, same 

341 units as this triaxial's semi-axes). 

342 @kwarg llk: The kind (an L{LLK}). 

343 @kwarg unit_name: Optional C{B{name}=NN} (C{str}) and scalar C{B{unit}=}L{Degrees} 

344 (or L{Radians}). 

345 

346 @return: A L{Cartesian5Tuple}C{(x, y, z, h, llk)} with height C{h=B{height}} and 

347 kind C{llk=B{llk}}. 

348 

349 @see: Method L{Triaxial3.reverseLatLon}. 

350 ''' 

351 _fwd = self.forwardBetOmg if llk in LLK._NOIDAL else \ 

352 self.forwardPhiLam # PYCHOK OK 

353 return _fwd(lat, lon, height=height, llk=llk, **_xkwds(unit_name, unit=Degrees)) 

354 

355 def forwardPhiLam(self, phi, lam, height=0, llk=LLK.GEODETIC, **unit_name): # generictocart2 

356 '''Convert any lat-/longitude kind to a cartesian on this triaxial's surface. 

357 

358 @arg phi: Latitude (C{Ang} or B{C{unit}}). 

359 @arg lam: Longitude (C{Ang} or B{C{unit}}). 

360 @kwarg height: Height above or below this triaxial's surface (C{meter}, same 

361 units as this triaxial's semi-axes). 

362 @kwarg llk: The kind (an L{LLK}). 

363 @kwarg unit_name: Optional C{B{name}=NN} (C{str}) and scalar C{B{unit}=}L{Radians} 

364 (or L{Degrees}). 

365 

366 @return: A L{Cartesian5Tuple}C{(x, y, z, h, llk)} with height C{h=B{height}} and 

367 kind C{llk=B{llk}}. 

368 

369 @note: Longitude C{B{lam} -= Lon0} if C{B{llk} is LLK.GEODETIC_LON0}. 

370 

371 @see: Method L{Triaxial3.reverseLatLon}. 

372 ''' 

373 ct, _ = self.forwardPhiLamZet2(phi, lam, None, height=height, llk=llk, **unit_name) 

374 return ct 

375 

376 def forwardPhiLamZet2(self, phi, lam, zet, height=0, llk=LLK.GEODETIC, **unit_name): # generictocart2 

377 '''Convert a lat-, longitude and heading to a cartesian and a direction 

378 on this triaxial's surface. 

379 

380 @arg phi: Latitude (C{Ang} or B{C{unit}}). 

381 @arg lam: Longitude (C{Ang} or B{C{unit}}). 

382 @arg zet: Azimuth of the heading (C{Ang}, B{C{unit}} or C{None}). 

383 @kwarg height: Height above or below this triaxial's surface (C{meter}, 

384 same units as this triaxial's semi-axes). 

385 @kwarg llk: The kind (an L{LLK}). 

386 @kwarg unit_name: Optional C{B{name}=NN} (C{str}) and scalar 

387 C{B{unit}=}L{Radians} (or L{Degrees}). 

388 

389 @return: 2-Tuple C{(cartesian, direction)} with the C{cartesian} a 

390 L{Cartesian5Tuple}C{(x, y, z, h, llk)} with height C{h=0}, 

391 kind C{llk=B{llk}} and C{direction}, a C{Vector3d} on and 

392 tangent to this triaxial's surface. 

393 

394 @note: Longitude C{B{lam} -= Lon0} if C{B{llk} is LLK.GEODETIC_LON0}. 

395 

396 @see: Method L{Triaxial3.reversePhiLamZet}. 

397 ''' 

398 unit, name = _xkwds_pop2(unit_name, unit=Radians) 

399 try: 

400 sa, ca = _SinCos2(phi, unit) 

401 if llk is LLK.GEODETIC_LON0 and self.Lon0: 

402 lam = Ang.fromScalar(lam, unit=unit) 

403 lam -= self.Lon0 

404 sb, cb = _SinCos2(lam, unit) 

405 except Exception as X: 

406 raise TriaxialError(phi=phi, lam=lam, llk=llk, cause=X) 

407 v, _, llk, name = _v_h_llk_name(ca * cb, ca * sb, sa, llk=llk, **name) 

408 if llk and llk._X: 

409 v = v._roty(False) # -1 

410 d, t = _d_t(self, llk) 

411 if t: 

412 v = v.times_(*t) 

413 if d: 

414 d = v.dividedBy_(*self._abc3).length 

415 v = v.dividedBy(d) 

416 

417 h = _HeightINT0(height) 

418 if h: # cart2cart 

419 v, h = self._toHeight2(v, h) 

420 ct = Cartesian5Tuple(v.x, v.y, v.z, h, llk, **name) 

421 

422 if zet is None: 

423 dir3d = None 

424 else: 

425 try: 

426 s, c = _SinCos2(zet, unit) 

427 except Exception as X: 

428 raise TriaxialError(zet=zet, cause=X) 

429 n, e = self._meridian2(v, lam, llk) 

430 dir3d = n.times(c) 

431 dir3d += e.times(s) 

432 dir3d.name = ct.name 

433 return ct, dir3d 

434 

435 def _meridian(self, lam, llk): 

436 '''(INTERNAL) Get the meridian plane's at C{lam}. 

437 ''' 

438 _, t = _d_t(self, llk) 

439 if t: 

440 a, b, c = t 

441 lam = lam.mod((c if llk._X else a) / b) 

442 return lam 

443 

444 def _meridian2(self, v, lam, llk): 

445 '''(INTERNAL) Get 2-tuple C{(n, e)} at C{lam} meridian. 

446 ''' 

447 mer = self._meridian(lam, llk) 

448 return self._cardinal2(v, mer, llk) 

449 

450 def normed2(self, x_ct, y=None, z=None, dir3d=None, **llk_name): # Ellipsoid3.Norm 

451 '''Scale a cartesian and direction to this triaxial's surface. 

452 

453 @arg x_ct: X component (C{scalar}) or a cartesian (L{Cartesian5Tuple} or 

454 any C{Cartesian}, L{Ecef9Tuple}, L{Vector3d}, L{Vector3Tuple} 

455 or L{Vector4Tuple}). 

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

457 ignored otherwise. 

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

459 @kwarg dir3d: The direction (C{Vector3d} or C{None}). 

460 @kwarg llk_name: Optional C{B{name}=NN} (C{str}) and kind C{B{llk}=None} 

461 overriding C{B{x_ct}.llk}. 

462 

463 @return: 2-Tuple C{(cartesian, direction)} with the C{cartesian} a 

464 L{Cartesian5Tuple}C{(x, y, z, h, llk)} and C{direction}, a 

465 C{Vector3d} tangent to this triaxial's surface or C{None} 

466 iff C{B{dir3d} is None}. 

467 ''' 

468 v, h, llk, name = _v_h_llk_name(x_ct, y, z, **llk_name) 

469 

470 u = v.dividedBy_(*self._abc3).length 

471 r = v.dividedBy(u) if u else _NAN3d 

472 ct = Cartesian5Tuple(r.x, r.y, r.z, h, llk, **name) 

473 

474 if isinstance(dir3d, Vector3d): 

475 if u: # and r is not _NAN3d 

476 u = r.dividedBy_(*self._a2b2c23) 

477 d = dir3d.dot(u) 

478 if _isfinite(d) and u.length2: 

479 u = u.times(d / u.length2) 

480 dir3d = dir3d.minus(u).unit() # NAN 

481 else: 

482 dir3d = _NAN3d 

483 else: 

484 dir3d = _NAN3d 

485 dir3d.name = ct.name 

486 return ct, dir3d 

487 

488 def reverse(self, x_ct, y=None, z=None, **name): 

489 '''Convert a cartesian I{on this triaxial's surface} to I{geodetic} lat- 

490 and longitude. 

491 

492 @arg x_ct: X component (C{scalar}) or a cartesian (L{Cartesian5Tuple} or 

493 any C{Cartesian}, L{Ecef9Tuple}, L{Vector3d}, L{Vector3Tuple} 

494 or L{Vector4Tuple}). 

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

496 ignored otherwise. 

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

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

499 

500 @return: A L{LatLonAzi5Tuple}C{(lat, lon, azimuth, h, llk)} with C{azimuth} 

501 usually C{INT0}, height C{h} off this triaxial's surface in C{meter} 

502 or C{INT0} and C{llk} either C{LLK.GEODETIC} or C{LLK.GEODETIC_LON0}. 

503 ''' 

504 llk = LLK.GEODETIC_LON0 if self.Lon0 else LLK.GEODETIC 

505 r = self.reversePhiLamZet(x_ct, y, z, llk=llk, **name) 

506 return r.toLatLonAzi5(name=r.name) 

507 

508 def reverseBetOmg(self, x_ct, y=None, z=None, **llk_name): # Cartesian3.carttoellip 

509 '''Convert a cartesian I{on this triaxial's surface} to an I{ellipsoidal} 

510 lat-/longitude. 

511 

512 @arg x_ct: X component (C{scalar}) or a cartesian (L{Cartesian5Tuple} or 

513 any C{Cartesian}, L{Ecef9Tuple}, L{Vector3d}, L{Vector3Tuple} 

514 or L{Vector4Tuple}). 

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

516 ignored otherwise. 

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

518 @kwarg llk_name: Optional C{B{name}=NN} (C{str}) and kind C{B{llk}=None} 

519 overriding C{B{x_ct}.llk}. 

520 

521 @return: A L{BetOmgAlp5Tuple}C{(bet, omg, alp, h, llk)} with C{alp=None} 

522 and C{llk=LLK.ELLIPSOIDAL}. 

523 ''' 

524 v, _, llk, name = _v_h_llk_name_NOIDAL(x_ct, y, z, **llk_name) 

525 

526 _, y2, z2 = rs = v.x2y2z23 

527 l0, l1, _ = ls = self._lcc23 

528 qmax = fsumf_(*rs) 

529 qmin = q = max(z2, y2 + z2 - l1, qmax - l0) 

530 _fp2 = _Fp2(rs, ls, n=1) 

531 f, _ = _fp2(q) 

532 if f > _TOL2: # neg means convergence 

533 q = max(qmin, min(qmax, _cubic(rs, qmax, l0, l1))) 

534 f, fp = _fp2(q) 

535 if fabs(f) > _TOL2: 

536 q = max(qmin, q - _over(f, fp)) 

537 q = _solve(_fp2, q, self.b2) 

538 

539 a, b, c = map1(_sqrt0, l0 + q, l1 + q, q) # axes (a, b, c) 

540 h = (c - self.c) or INT0 

541 bet, omg, _ = self._reverseBetOmgAlp3(v, None, a, b, c, **name) 

542 return BetOmgAlp5Tuple(bet, omg, None, h, llk, **name) 

543 

544 reverseBetaOmega = reverseBetOmg # for backward compatibility 

545 

546 def reverseBetOmgAlp(self, x_ct, y=None, z=None, dir3d=None, **llk_name): # Ellipsoid3.cart2toellip[int] 

547 '''Convert a cartesian and direction I{on this triaxial's surface} to an 

548 I{ellipsoidal} lat-, longitude and heading. 

549 

550 @arg x_ct: X component (C{scalar}) or a cartesian (L{Cartesian5Tuple} or 

551 any C{Cartesian}, L{Ecef9Tuple}, L{Vector3d}, L{Vector3Tuple} 

552 or L{Vector4Tuple}). 

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

554 ignored otherwise. 

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

556 @kwarg dir3d: The direction (C{Vector3d} or C{None}). 

557 @kwarg llk_name: Optional C{B{name}=NN} (C{str}) and kind C{B{llk}=None} 

558 overriding C{B{x_ct}.llk}. 

559 

560 @return: A L{BetOmgAlp5Tuple}C{(bet, omg, alp, h, llk)} with C{alp=None} 

561 if C{B{dir3d} is None} and C{llk=LLK.ELLIPSOIDAL}. 

562 ''' 

563 v, h, llk, name = _v_h_llk_name_NOIDAL(x_ct, y, z, **llk_name) 

564 bet, omg, alp = self._reverseBetOmgAlp3(v, dir3d, **name) 

565 return BetOmgAlp5Tuple(bet, omg, alp, h, llk, **name) 

566 

567 def _reverseBetOmgAlp3(self, v, dir3d, *a_b_c, **name): # cart2toellipint 

568 '''(INTERNAL) Helper for methods C{reverseBetOmg/-Alp}. 

569 ''' 

570 k, kp = self._k_kp 

571 k2, kp2 = self._k2_kp2 

572 a, b, c = a_b_c or self._abc3 

573 V = v.dividedBy_(a, b, c) 

574 X, E, Z = V.xyz3 # Xi, Eta, Zeta 

575 h = fabs(E * k * kp * _2_0) 

576 if v.y or fabs(v.x) != a * kp2 or \ 

577 fabs(v.z) != c * k2: 

578 g = fdot(V.x2y2z23, k2, (k2 - kp2), -kp2) 

579 h = hypot(g, h) 

580 else: 

581 g = _0_0 

582 if h < EPS02: 

583 so = cb = _0_0 

584 elif g < 0: 

585 h = _over(sqrt((h - g) * _0_5), kp) 

586 so = _copysign(h, E) 

587 cb = fabs(_over(E, so)) 

588 else: 

589 cb = _over(sqrt((h + g) * _0_5), k) 

590 so = _over(E, cb) 

591 tx, tz = _txtz2(cb, so, k, kp) 

592 sb = (Z / tz) if tz else _N_1_0 

593 co = (X / tx) if tx else _1_0 

594 bet = Ang_(sb, cb, **name) 

595 omg = Ang_(so, co, **name) 

596 

597 if isinstance(dir3d, Vector3d): # cart2toellip(bet, omg, V) -> alp 

598 if cb or so or not (tx and tz): # not umbilical 

599 if not tx: 

600 n = Vector3d(-co, -so, tx) * sb 

601 e = Vector3d(-so, co, tx) 

602 elif not tz: 

603 n = Vector3d(tz, -sb, cb) 

604 e = Vector3d(tz, cb, sb) * co 

605 else: 

606 n = Vector3d(-a * sb * k2 * cb * co / tx, 

607 -b * sb * so, c * cb * tz) 

608 e = Vector3d(-a * so * tx, b * cb * co, 

609 c * so * kp2 * sb * co / tz) 

610 sa = dir3d.dot(e.unit()) # NAN 

611 ca = dir3d.dot(n.unit()) # NAN 

612 else: # at umbilicial PYCHOK no cover 

613 x, z = norm2(co * tx / a, sb * tz / c) # _MODS.karney._norm2 

614 v = dir3d * (sb * co) # dir3d.times(sb * co) 

615 s2a = -v.y 

616 c2a = fdot(v, z, 0, -x) # v.x * z - v.z * x 

617 sa = ca = -sb 

618 sa *= _copysign(_1_0 - c2a, s2a) if c2a < 0 else s2a 

619 ca *= fabs(s2a) if c2a < 0 else (c2a + _1_0) 

620 alp = Ang_(sa, ca, **name) 

621 elif dir3d is None: 

622 alp = None # Ang.NAN(**name) 

623 else: 

624 raise TriaxialError(dir3d=dir3d) 

625 return bet, omg, alp 

626 

627 def reverseCartesian(self, x_ct, y=None, z=None, height=0, normal=True, **llk_name): # cart2tocart 

628 '''"Unproject" a cartesian I{off} this triaxial's surface. 

629 

630 @arg x_ct: X component (C{scalar}) or a cartesian (L{Cartesian5Tuple} or 

631 any C{Cartesian}, L{Ecef9Tuple}, L{Vector3d}, L{Vector3Tuple} 

632 or L{Vector4Tuple}). 

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

634 ignored otherwise. 

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

636 @kwarg height: Height above or below this triaxial's surface (C{meter}, 

637 same units as this triaxial's semi-axes). 

638 @kwarg normal: If C{True}, B{C{height}} is C{perpendicular} to the surface, 

639 otherwise C{radial} to the center of this triaxial (C{bool}). 

640 @kwarg llk_name: Optional C{B{name}=NN} (C{str}) and kind C{B{llk}} 

641 overriding C{B{x_ct}.llk}. 

642 

643 @return: L{Cartesian5Tuple}C{(x, y, z, h, llk)}. 

644 

645 @raise TrialError: Cartesian B{C{x_ct}} or C{(x, y, z)} not on this 

646 triaxial's surface. 

647 

648 @see: Methods L{Triaxial3.forwardCartesian}. 

649 ''' 

650 kwds = _xkwds(llk_name, llk=_xattr(x_ct, llk=None)) 

651 v, _, llk, name = _v_h_llk_name(x_ct, y, z, **kwds) 

652 _ = self._sideOn(v) 

653 h = _HeightINT0(height) 

654 if h: 

655 v, h = self._toHeight2(v, h, normal) 

656 return Cartesian5Tuple(v.x, v.y, v.z, h, llk, **name) 

657 

658 def reverseLatLon(self, x_ct, y=None, z=None, **llk_name): # cart2toany 

659 '''Convert a cartesian I{on this triaxial's surface} to a lat-/longitude. 

660 

661 @arg x_ct: X component (C{scalar}) or a cartesian (L{Cartesian5Tuple} or 

662 any C{Cartesian}, L{Ecef9Tuple}, L{Vector3d}, L{Vector3Tuple} 

663 or L{Vector4Tuple}). 

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

665 ignored otherwise. 

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

667 @kwarg llk_name: Optional C{B{name}=NN} (C{str}) and kind C{B{llk}=None} 

668 overriding C{B{x_ct}.llk}. 

669 

670 @return: A L{BetOmgAlp5Tuple}C{(bet, omg, alp, h, llk)} with C{alp=None} or 

671 a L{PhiLamZet5Tuple}C{(phi, lam, zet, h, llk)} with C{zet=None}. 

672 

673 @note: Longitude C{B{lam} += Lon0} if C{B{llk} is LLK.GEODETIC_LON0}. 

674 ''' 

675 llk, name = _xkwds_pop2(llk_name, llk=_xattr(x_ct, 

676 llk=LLK.ELLIPSOIDAL)) 

677 _rev = self.reverseBetOmg if llk in LLK._NOIDAL else \ 

678 self.reversePhiLam # PYCHOK OK 

679 return _rev(x_ct, y, z, llk=llk, **name) 

680 

681 def reversePhiLam(self, x_ct, y=None, z=None, **llk_name): # cart2togeneric 

682 '''Convert a cartesian I{on this triaxial's surface} to lat-/longitude. 

683 

684 @arg x_ct: X component (C{scalar}) or a cartesian (L{Cartesian5Tuple} or 

685 any C{Cartesian}, L{Ecef9Tuple}, L{Vector3d}, L{Vector3Tuple} 

686 or L{Vector4Tuple}). 

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

688 ignored otherwise. 

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

690 @kwarg llk_name: Optional C{B{name}=NN} (C{str}) and kind C{B{llk}=None} 

691 overriding C{B{x_ct}.llk}. 

692 

693 @return: A L{PhiLamZet5Tuple}C{(phi, lam, zet, h, llk)} with C{zet=None}. 

694 

695 @note: Longitude C{B{lam} += Lon0} if C{B{llk} is LLK.GEODETIC_LON0}. 

696 ''' 

697 return self.reversePhiLamZet(x_ct, y, z, **llk_name) 

698 

699 def reversePhiLamZet(self, x_ct, y=None, z=None, dir3d=None, **llk_name): # cart2togeneric(R, V, ... 

700 '''Convert a cartesian and direction to lat-, longitude and azimuth. 

701 

702 @arg x_ct: X component (C{scalar}) or a cartesian (L{Cartesian5Tuple} or 

703 any C{Cartesian}, L{Ecef9Tuple}, L{Vector3d}, L{Vector3Tuple} 

704 or L{Vector4Tuple}). 

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

706 ignored otherwise. 

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

708 @kwarg dir3d: Optional direction (C{Vector3d} or C{None}). 

709 @kwarg llk_name: Optional C{B{name}=NN} (C{str}) and kind C{B{llk}=None} 

710 overriding C{B{x_ct}.llk}. 

711 

712 @return: A L{PhiLamZet5Tuple}C{(phi, lam, zet, h, llk)} with C{zet=None} 

713 if C{B{dir3d} is None}. 

714 

715 @note: Longitude C{B{lam} += Lon0} if C{B{llk} is LLK.GEODETIC_LON0}. 

716 ''' 

717 ct = self.toTriaxial5(x_ct, y, z, h=NAN, **llk_name) 

718 v, h, llk, name = _v_h_llk_name(ct) 

719 _, t = _d_t(self, llk) 

720 if t: 

721 v = v.dividedBy_(*t) 

722 if llk._X: 

723 v = v._roty(True) # +1 

724 phi = Ang_(v.z, hypot(v.x, v.y), **name) 

725 lam = Ang_(v.y, v.x, **name) # Ang(0, 0) -> 0 

726 

727 if dir3d is None: 

728 zet = None 

729 elif isinstance(dir3d, Vector3d): 

730 n, e = self._meridian2(v, lam, llk) 

731 zet = Ang_(dir3d.dot(e), 

732 dir3d.dot(n), **name) 

733 else: 

734 raise TriaxialError(dir3d=dir3d) 

735 if llk is LLK.GEODETIC_LON0 and self.Lon0: 

736 lam += self.Lon0 

737 return PhiLamZet5Tuple(phi, lam, zet, h, llk, **name) 

738 

739 def random2(self, llk=LLK.ELLIPSOIDAL, both=False, _rand=random): 

740 '''Return a random cartesian with/out direction on this triaxial's surface. 

741 

742 @kwarg llk: The kind (an L{LLK}). 

743 @kwarg both: If C{True}, generate a random direction (C{bool}). 

744 

745 @return: 2-Tuple C{(cartesian, direction)} with the C{cartesian} a 

746 L{Cartesian5Tuple}C{(x, y, z, h, llk)} and C{direction}, a 

747 C{Vector3d} tangent to this triaxial's surface or C{None} 

748 iff C{B{both} is False}. 

749 ''' 

750 for _ in range(_MAXIT): 

751 for _ in range(_MAXIT): 

752 v = Vector3d(_rand(), _rand(), _rand()) 

753 u = v.length 

754 if u and _isfinite(u): 

755 break 

756 else: 

757 raise TriaxialError(Fmt.no_convergence(u)) 

758 v = v.dividedBy(u).times_(*self._abc3) 

759 q = v.dividedBy_(*self._a2b2c23).length * self.c 

760 if 0 < q <= _1_0: # _uni(q) < q: 

761 break 

762 else: 

763 raise TriaxialError(Fmt.no_convergence(q)) 

764 ct = Cartesian5Tuple(v.x, v.y, v.z, INT0, llk, name__=self.random2) 

765 v = None 

766 if both: 

767 for _ in range(_MAXIT): 

768 v = Vector3d(_rand(), _rand(), _rand()) 

769 u = v.length 

770 if u: