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, _DMAIN_, _h_, _lam_, _phi_ 

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

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

46from pygeodesy.namedTuples import Vector4Tuple, _Pass 

47# from pygeodesy.props import Property_RO # from .units 

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

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

50 _MAXIT, _omg_, _otherV3d_, _sqrt0, \ 

51 _Triaxial3Base, TriaxialError, \ 

52 _TriaxialsBase 

53from pygeodesy.units import Degrees, Radians, Radius_, Property_RO 

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

55from pygeodesy.vector3d import Vector3d, _ALL_LAZY 

56 

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

58from random import random 

59 

60__all__ = _ALL_LAZY.triaxials_triaxial3 

61__version__ = '26.02.20' 

62 

63_alp_ = 'alp' 

64_NAN3d = Vector3d(NAN, NAN, NAN) 

65_TOL = cbrt2(EPS) 

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

67_zet_ = 'zet' 

68_27_0 = 27.0 

69 

70 

71class BetOmgAlp5Tuple(_Ang3Tuple): 

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

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

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

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

76 ''' 

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

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

79 

80 

81class Cartesian5Tuple(Vector4Tuple): 

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

83 C{y} and C{z} coordinates on and height C{h} above 

84 or below the triaxial's surface and kind C{llk} set 

85 to the original C{LLK} or C{None}. 

86 ''' 

87 _Names_ = Vector4Tuple._Names_ + (_llk_,) 

88 _Units_ = Vector4Tuple._Units_ + (_Pass,) 

89 

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

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

92 return Vector4Tuple.__new__(cls, args, **kwds) 

93 

94 

95class _Fp2(object): 

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

97 ''' 

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

99 # assert 0 < n <= 2 

100 self._2 = n == 2 

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

102 

103 def __call__(self, p): 

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

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

106 f = _N_1_0 

107 fc = fp = _0_0 

108 _D = EPS_2 

109 _2 = self._2 

110 for g, q in self._rls: 

111 q = _1_over(p + q) 

112 g *= q 

113 if _2: 

114 g *= g 

115 q += q 

116 r = round(g / _D) * _D 

117 f += r 

118 fc += g - r 

119 fp -= g * q 

120 return (f + fc), fp 

121 

122 

123class PhiLamZet5Tuple(_Ang3Tuple): 

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

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

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

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

128 ''' 

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

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

131 

132 

133class Triaxial3(_Triaxial3Base): 

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

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

136 

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

138 ''' 

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

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

141 ''' 

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

143 # and isinstance(llk, LLK.__class__) 

144 a2, b2, c2 = self._a2b2c23 

145 if llk._X: 

146 a2, c2 = c2, a2 

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

148 x, y, z = v.xyz3 

149 if x or y: 

150 s = (-z) / c2 

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

152 else: 

153 y, x, _ = mer.scn3 

154 s = _copysign_1_0(-z) 

155 z = _0_0 

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

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

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

159 if llk._X: 

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

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

162 return n, e 

163 

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

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

166 on this triaxial's surface. 

167 

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

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

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

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

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

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

174 

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

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

177 

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

179 ''' 

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

181 return ct 

182 

183 forwardBetaOmega = forwardBetOmg # for backward compatibility 

184 

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

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

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

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

189 

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

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

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

193 

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

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

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

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

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

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

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

201 

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

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

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

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

206 

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

208 ''' 

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

210 a, b, c = self._abc3 

211 if h: # Cartesian.elliptocart 

212 a, b, _ = self._a2b2c23 

213 h = _HeightINT0(h) 

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

215 if s < 0: 

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

217 a = sqrt(a + s) 

218 b = sqrt(b + s) 

219 c += h 

220 sb, cb = _SinCos2(bet, unit) 

221 so, co = _SinCos2(omg, unit) 

222 k, kp = self._k_kp 

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

224 ct = Cartesian5Tuple(a * co * tx, 

225 b * cb * so, 

226 c * sb * tz, 

227 h, llk, **name) 

228 

229 if alp is None: # or h? 

230 dir3d = None # _NAN3d? 

231 else: 

232 try: 

233 sa, ca = _SinCos2(alp, unit) 

234 except Exception as X: 

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

236 a, b, c = self._abc3 

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

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

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

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

241 c * kp * sb * c2s2_b) 

242 else: 

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

244 c = _flipsign(co, cb) 

245 n = Vector3d(-c * sb, 

246 -so * sb, _0_0) 

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

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

249 s = _flipsign(sb, so) 

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

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

252 else: 

253 k2, kp2 = self._k2_kp2 

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

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

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

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

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

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

260 dir3d.name = ct.name 

261 return ct, dir3d 

262 

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

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

265 

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

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

268 or L{Vector4Tuple}). 

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

270 ignored otherwise. 

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

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

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

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

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

276 

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

278 

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

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

281 ''' 

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

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

284 if h: # signed, square 

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

286 h = v.h 

287 else: # on the surface 

288 v = _otherV3d_(x_ct, y, z) 

289 n = _xkwds_get(kwds, name=NN) 

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

291 

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

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

294 

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

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

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

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

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

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

301 (or L{Radians}). 

302 

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

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

305 

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

307 ''' 

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

309 self.forwardPhiLam # PYCHOK OK 

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

311 

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

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

314 

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

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

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

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

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

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

321 (or L{Degrees}). 

322 

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

324 C{llk=B{llk}}. 

325 

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

327 

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

329 ''' 

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

331 return ct 

332 

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

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

335 on this trixial's surface. 

336 

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

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

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

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

341 same 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 

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

345 

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

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

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

349 tangent to this triaxial's surface. 

350 

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

352 

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

354 ''' 

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

356 try: 

357 sa, ca = _SinCos2(phi, unit) 

358 if llk is LLK.GEODETIC_LON0: 

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

360 lam -= self.Lon0 

361 sb, cb = _SinCos2(lam, unit) 

362 except Exception as X: 

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

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

365 if llk and llk._X: 

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

367 d, t = _d_t(self, llk) 

368 if t: 

369 v = v.times_(*t) 

370 if d: 

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

372 v = v.dividedBy(d) 

373 

374 h = _HeightINT0(height) 

375 if h: # cart2cart 

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

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

378 

379 if zet is None: 

380 dir3d = None 

381 else: 

382 try: 

383 s, c = _SinCos2(zet, unit) 

384 except Exception as X: 

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

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

387 dir3d = n.times(c) 

388 dir3d += e.times(s) 

389 dir3d.name = ct.name 

390 return ct, dir3d 

391 

392 def _meridian(self, lam, llk): 

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

394 ''' 

395 _, t = _d_t(self, llk) 

396 if t: 

397 a, b, c = t 

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

399 return lam 

400 

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

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

403 ''' 

404 mer = self._meridian(lam, llk) 

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

406 

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

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

409 

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

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

412 or L{Vector4Tuple}). 

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

414 ignored otherwise. 

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

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

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

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

419 

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

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

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

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

424 ''' 

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

426 

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

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

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

430 

431 if isinstance(dir3d, Vector3d): 

432 if u: # and r is not _NAN3d 

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

434 d = dir3d.dot(u) 

435 if _isfinite(d) and u.length2: 

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

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

438 else: 

439 dir3d = _NAN3d 

440 else: 

441 dir3d = _NAN3d 

442 dir3d.name = ct.name 

443 return ct, dir3d 

444 

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

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

447 lat-/longitude. 

448 

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

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

451 or L{Vector4Tuple}). 

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

453 ignored otherwise. 

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

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

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

457 

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

459 and C{llk=LLK.ELLIPSOIDAL}. 

460 ''' 

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

462 

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

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

465 qmax = fsumf_(*rs) 

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

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

468 f, _ = _fp2(q) 

469 if f > _TOL2: # neg means convergence 

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

471 f, fp = _fp2(q) 

472 if fabs(f) > _TOL2: 

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

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

475 

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

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

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

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

480 

481 reverseBetaOmega = reverseBetOmg # for backward compatibility 

482 

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

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

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

486 

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

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

489 or L{Vector4Tuple}). 

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

491 ignored otherwise. 

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

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

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

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

496 

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

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

499 ''' 

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

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

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

503 

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

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

506 ''' 

507 k, kp = self._k_kp 

508 k2, kp2 = self._k2_kp2 

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

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

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

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

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

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

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

516 h = hypot(g, h) 

517 else: 

518 g = _0_0 

519 if h < EPS02: 

520 so = cb = _0_0 

521 elif g < 0: 

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

523 so = _copysign(h, E) 

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

525 else: 

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

527 so = _over(E, cb) 

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

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

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

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

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

533 

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

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

536 if not tx: 

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

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

539 elif not tz: 

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

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

542 else: 

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

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

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

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

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

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

549 else: # at umbilicial PYCHOK no cover 

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

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

552 s2a = -v.y 

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

554 sa = ca = -sb 

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

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

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

558 elif dir3d is None: 

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

560 else: 

561 raise TriaxialError(dir3d=dir3d) 

562 return bet, omg, alp 

563 

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

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

566 

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

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

569 or L{Vector4Tuple}). 

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

571 ignored otherwise. 

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

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

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

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

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

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

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

579 

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

581 

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

583 triaxial's surface. 

584 

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

586 ''' 

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

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

589 _ = self._sideOn(v) 

590 h = _HeightINT0(height) 

591 if h: 

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

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

594 

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

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

597 

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

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

600 or L{Vector4Tuple}). 

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

602 ignored otherwise. 

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

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

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

606 

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

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

609 

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

611 ''' 

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

613 llk=LLK.ELLIPSOIDAL)) 

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

615 self.reversePhiLam # PYCHOK OK 

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

617 

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

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

620 

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

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

623 or L{Vector4Tuple}). 

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

625 ignored otherwise. 

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

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

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

629 

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

631 

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

633 ''' 

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

635 

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

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

638 

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

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

641 or L{Vector4Tuple}). 

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

643 ignored otherwise. 

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

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

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

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

648 

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

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

651 

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

653 ''' 

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

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

656 _, t = _d_t(self, llk) 

657 if t: 

658 v = v.dividedBy_(*t) 

659 if llk._X: 

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

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

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

663 

664 if dir3d is None: 

665 zet = None 

666 elif isinstance(dir3d, Vector3d): 

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

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

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

670 else: 

671 raise TriaxialError(dir3d=dir3d) 

672 if llk is LLK.GEODETIC_LON0: 

673 lam += self.Lon0 

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

675 

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

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

678 

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

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

681 

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

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

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

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

686 ''' 

687 for _ in range(_MAXIT): 

688 for _ in range(_MAXIT): 

689 v = Vector3d(_rand(), _rand(), _rand()) 

690 u = v.length 

691 if u and _isfinite(u): 

692 break 

693 else: 

694 raise TriaxialError(Fmt.no_convergence(u)) 

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

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

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

698 break 

699 else: 

700 raise TriaxialError(Fmt.no_convergence(q)) 

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

702 v = None 

703 if both: 

704 for _ in range(_MAXIT): 

705 v = Vector3d(_rand(), _rand(), _rand()) 

706 u = v.length 

707 if u: 

708 u = v.dividedBy(u).dividedBy_(*self._a2b2c23) 

709 d = v.dot(u) / u.length2 

710 v = v.minus(u.times(d)) 

711 u = v.length # normvec 

712 if u and _isfinite(u): 

713 v = v.dividedBy(u) 

714 break 

715 else: 

716 raise TriaxialError(Fmt.no_convergence(u)) 

717 v.name = ct.name 

718 return ct, v 

719 

720 def _toHeight2(self, v, h, normal=True): 

721 '''(INTERNAL) Move cartesian C{Vector3d B{v}} to height C{h}. 

722 ''' 

723 n = v.dividedBy_(*self._a2b2c23) if normal else v 

724 if n.length > EPS02: 

725 h = max(h, -self.c) 

726 v = v.plus(n.times(h / n.length)) 

727 return v, h 

728 

729 def toOther(self, lat, lon, llk1=LLK.GEODETIC, llk2=LLK.GEODETIC, **unit_name): # anytoany 

730 '''Convert one lat-/longitude kind to an other. 

731 

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

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

734 @kwarg llk1: The given kind (an L{LLK}). 

735 @kwarg llk2: The result kind (an L{LLK}). 

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

737 

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

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

740 

741 @see: Methods L{Triaxial3.forwardLatLon} and -L{reverseLatLon}. 

742 ''' 

743 ct = self.forwardLatLon(lat, lon, llk=llk1, **unit_name) 

744 r = self.reverseLatLon(ct, llk=llk2, name=ct.name) 

745# a, b = r[:2] 

746# if not isAng(lat): 

747# a = float(a) 

748# if not isAng(lon): 

749# b = float(b) 

750# if (a, b) =! r[:2]: 

751# r = r._dup(a, b) 

752 return r 

753 

754 def toTriaxial5(self, x_ct, y=None, z=None, **triaxial_h_llk_name): # carttocart2 

755 '''Find the closest cartesian on this or on another triaxial's surface. 

756 

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

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

759 or L{Vector4Tuple}). 

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

761 ignored otherwise. 

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

763 @kwarg triaxial_llk_name: Optional C{B{triaxial}=self} (C{Triaxial3}), 

764 C{B{name}=NN} (C{str}), height C{B{h}} and kind C{B{llk}} 

765 overriding C{B{x_ct}.h} respectively C{B{x_ct}.llk}. 

766 

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

768 

769 @raise TriaxialError: If C{B{triaxial}} is not a L{Triaxial3}. 

770 

771 @see: Functions L{hartzell4<triaxials.triaxial5.hartzell4>} and 

772 L{height4<triaxials.triaxial5.height4>} and methods.