1

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

3

4 u'''Triaxal ellipsoid classes L{JacobiConformal}, Jacobi's conformal projection, trancoded

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

6 doc/classGeographicLib_1_1JacobiConformal.html#details>} to pure Python, I{ordered} L{Triaxial}

7 and I{unordered} L{Triaxial_} and miscellaneous classes L{BetaOmega2Tuple}, L{BetaOmega3Tuple},

8 L{Jacobi2Tuple} and L{TriaxialError}.

9

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

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

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

13

14@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)

15@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)

16@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)

17@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)

18@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)

19@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)

20@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)

21@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)

22@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)

23@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)

24@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)

25@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)

26'''

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

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

29

30# from pygeodesy.basics import isscalar, _zip # from .fsums, .namedTuples, .streprs

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

32 _0_0, _0_5, _1_0, _N_2_0, isfinite, isnear1, \

33 _4_0 # PYCHOK used!

34 from pygeodesy.datums import Datum, Ellipsoid, _spherical_datum, _WGS84

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

36# from pygeodesy.elliptic import Elliptic # ._MODS

37# from pygeodesy.errors import _ValueError # from .streprs

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

39 from pygeodesy.fsums import Fsum, fsum_, isscalar, Property_RO

40 from pygeodesy.interns import NN, _a_, _b_, _beta_, _c_, _distant_, _height_, \

41 _inside_, _near_, _not_, _NL_, _NLATvar_, _NOTEQUAL_, \

42 _null_, _opposite_, _outside_, _SPACE_, _spherical_, \

43 _too_, _x_, _y_, _COMMA_ # PYCHOK used!

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

45 from pygeodesy.named import _NamedBase, _NamedEnum, _NamedEnumItem, \

46 _NamedTuple, _Pass, _lazyNamedEnumItem as _lazy

47 from pygeodesy.namedTuples import LatLon3Tuple, Vector3Tuple, Vector4Tuple

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

49 from pygeodesy.streprs import Fmt, _ValueError, _zip

50 from pygeodesy.units import Degrees, Float, Height_, Meter, Meter2, Meter3, Radians, Radius

51 from pygeodesy.utily import asin1, atan2d, km2m, m2km, sincos2, sincos2d, sincos2d_

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

53

54 from math import atan2, fabs, sqrt

55

56__all__ = _ALL_LAZY.triaxials

57__version__ = '23.03.20'

58

59_E = _WGS84.ellipsoid

60_not_ordered_ = _not_('ordered')

61_omega_ = 'omega'

62_TRIPS = 537 # max 55, Eberly 1074?

63_WGS84_35abc = _E.a + 35, _E.a - 35, _E.b

64 del _E

65

66

67 class _ToNamedBase(_NamedBase):

68 '''(INTERNAL) C{-.toDegrees}, C{-.toRadians} base.

69 '''

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

71 if toDMS_kwds:

72 toDMS = _MODS.dms.toDMS

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

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

75 elif isinstance(a, Degrees) and \

76 isinstance(b, Degrees):

77 return self

78 else:

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

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

81

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

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

84 isinstance(b, Radians) else \

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

86 *c, name=self.name)

87

88

89 class BetaOmega2Tuple(_NamedTuple, _ToNamedBase):

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

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

92 C{Degrees}).

93 '''

94 _Names_ = (_beta_, _omega_)

95 _Units_ = (_Pass, _Pass)

96

97 def toDegrees(self, **toDMS_kwds):

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

99

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

101 C{beta} and C{omega} both in C{Degrees}

102 or as an L{toDMS} string provided some

103 B{C{toDMS_kwds}} are supplied.

104 '''

105 return _ToNamedBase._toDegrees(self, *self, **toDMS_kwds)

106

107 def toRadians(self):

108 '''Convert this L{BetaOmega2Tuple} to C{Radians}.

109

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

111 C{beta} and C{omega} both in C{Radians}.

112 '''

113 return _ToNamedBase._toRadians(self, *self)

114

115

116 class BetaOmega3Tuple(_NamedTuple, _ToNamedBase):

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

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

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

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

121 in C{meter}, conventionally.

122 '''

123 _Names_ = BetaOmega2Tuple._Names_ + (_height_,)

124 _Units_ = BetaOmega2Tuple._Units_ + ( Meter,)

125

126 def toDegrees(self, **toDMS_kwds):

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

128

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

130 C{beta} and C{omega} both in C{Degrees} or as an

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

132 are supplied.

133 '''

134 return _ToNamedBase._toDegrees(self, *self, **toDMS_kwds)

135

136 def toRadians(self):

137 '''Convert this L{BetaOmega3Tuple} to C{Radians}.

138

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

140 C{beta} and C{omega} both in C{Radians}.

141 '''

142 return _ToNamedBase._toRadians(self, *self)

143

144 def to2Tuple(self):

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

146 '''

147 return BetaOmega2Tuple(*self[:2])

148

149

150 class Jacobi2Tuple(_NamedTuple, _ToNamedBase):

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

152 projection, both in C{Radians} (or C{Degrees}).

153 '''

154 _Names_ = (_x_, _y_)

155 _Units_ = (_Pass, _Pass)

156

157 def toDegrees(self, **toDMS_kwds):

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

159

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

161 both in C{Degrees} or as an L{toDMS} string

162 provided some B{C{toDMS_kwds}} are supplied.

163 '''

164 return _ToNamedBase._toDegrees(self, *self, **toDMS_kwds)

165

166 def toRadians(self):

167 '''Convert this L{Jacobi2Tuple} to C{Radians}.

168

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

170 and C{y} both in C{Radians}.

171 '''

172 return _ToNamedBase._toRadians(self, *self)

173

174

175 class Triaxial_(_NamedEnumItem):

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

177

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

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

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

181 I{omega}=0.

182

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

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

185

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

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

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

189 C{Degrees} if converted).

190 '''

191 _ijk = _kji = None

192 _unordered = True

193

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

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

196

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

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

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

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

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

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

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

204

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

206 '''

207 try:

208 a = a_triaxial

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

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

211 except (TypeError, ValueError) as x:

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

213 if name:

214 self.name = name

215

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

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

218 s, _, t = sorted(t)

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

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

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

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

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

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

225

226 def __str__(self):

227 return self.toStr()

228

229 @Property_RO

230 def a(self):

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

232 '''

233 a, _, _ = self._abc3

234 return a

235

236 @Property_RO

237 def _a2b2(self):

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

239 '''

240 a, b, _ = self._abc3

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

242

243 @Property_RO

244 def _a2_b2(self):

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

246 '''

247 a, b, _ = self._abc3

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

249

250 @Property_RO

251 def _a2c2(self):

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

253 '''

254 a, _, c = self._abc3

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

256

257 @Property_RO

258 def area(self):

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

260 '''

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

262 if a > c:

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

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

265 else: # a == c == b

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

267 return a

268

269 def area_p(self, p=1.6075):

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

271

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

273 for "near-flat" triaxials.

274

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

276 '''

277 a, b, c = self._abc3

278 if a == b == c:

279 a *= a

280 else:

281 _p = pow

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

283 return Meter2(area_p=a * PI4)

284

285 @Property_RO

286 def b(self):

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

288 '''

289 _, b, _ = self._abc3

290 return b

291

292 @Property_RO

293 def _b2c2(self):

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

295 '''

296 _, b, c = self._abc3

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

298

299 @Property_RO

300 def c(self):

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

302 '''

303 _, _, c = self._abc3

304 return c

305

306 @Property_RO

307 def _c2_b2(self):

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

309 '''

310 _, b, c = self._abc3

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

312

313 @Property_RO

314 def e2ab(self):

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

316 '''

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

318

319 @Property_RO

320 def _1e2ab(self):

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

322 '''

323 a, b, _ = self._abc3

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

325

326 @Property_RO

327 def e2bc(self):

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

329 '''

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

331

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

333

334 @Property_RO

335 def e2ac(self):

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

337 '''

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

339

340 @Property_RO

341 def _1e2ac(self):

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

343 '''

344 a, _, c = self._abc3

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

346

347 @Property_RO

348 def _Elliptic(self):

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

350 '''

351 return _MODS.elliptic.Elliptic

352

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

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

355 from a Point-Of-View in space.

356

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

358 '''

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

360

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

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

363 this triaxial's surface.

364

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

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

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

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

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

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

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

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

373 negative value to skip validation.

374

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

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

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

378 surface in C{meter}, conventionally.

379

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

381 convergence in root finding or validation failed.

382

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

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

385 DistancePointEllipseEllipsoid.pdf>}.

386 '''

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

388

389 i, h = None, v.length

390 if h < EPS0: # EPS

391 x = y = z = _0_0

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

393 elif r: # .isSpherical

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

395 h -= r

396 else:

397 x, y, z = v.xyz

398 try:

399 if normal: # perpendicular to triaxial

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

401 else: # radially to triaxial's center

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

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

404 except Exception as e:

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

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

407 h = -h # below the surface

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

409

410 @Property_RO

411 def isOrdered(self):

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

413 '''

414 a, b, c = self._abc3

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

416

417 @Property_RO

418 def isSpherical(self):

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

420 '''

421 a, b, c = self._abc3

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

423

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

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

426

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

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

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

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

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

432

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

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

435

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

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

438 '''

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

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

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

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

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

444 if n.length < EPS0:

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

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

447

448 @Property_RO

449 def _normal3d(self):

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

451 '''

452 d = max(self._abc3)

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

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

455

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

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

458 '''

459 if fabs(s) > _1_0 or fabs(c) > _1_0 or \

460 fabs(_hypot21_(s, c)) > EPS0:

461 s, c = norm2(s, c)

462 if a:

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

464 return (s or _0_0), (c or _0_0)

465

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

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

468

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

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

471 '''

472 ijk = self._order_ijk(**reverse)

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

474

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

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

477

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

479 '''

480 ijk = self._order_ijk(**reverse)

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

482

483 @Property_RO

484 def _ordered4(self):

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

486 '''

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

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

489

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

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

492 '''

493 i, j, k = 0, 1, 2

494 if a < b:

495 a, b, i, j = b, a, j, i

496 if a < c:

497 a, c, i, k = c, a, k, i

498 if b < c:

499 b, c, j, k = c, b, k, j

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

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

502 return (a, b, c), ijk

503

504 abc, T = self._abc3, self

505 if not self.isOrdered:

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

507 if ijk:

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

509 T = Triaxial_(*abc)

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

511 return abc + (T,)

512

513 def _order_ijk(self, reverse=False):

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

515 '''

516 return self._kji if reverse else self._ijk

517

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

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

520 '''

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

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

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

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

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

526

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

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

529

530 a, b, c = self._abc3

531 if a != b:

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

533 if a != c:

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

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

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

537 return x, y, z

538

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

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

541

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

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

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

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

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

547

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

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

550 if in- respectively outside this triaxial.

551

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

553 '''

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

555

556 def _sqrt(self, x):

557 '''(INTERNAL) Helper.

558 '''

559 if x < 0:

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

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

562

563 def toEllipsoid(self, name=NN):

564 '''Convert this triaxial to an L{Ellipsoid}, provided C{a == b} or C{b == c}.

565

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

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

568

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

570

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

572 '''

573 a, b, c = self._abc3

574 if a == b: # N = Z

575 b = c

576 elif b == c: # N = X

577 a, b = b, a

578 elif a != c:

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

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

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

582

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

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

585

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

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

588 this triaxial's name.

589

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

591 '''

592 T = Triaxial_

593 t = T.a, T.b, T.c, T.e2ab, T.e2bc, T.e2ac

594 if isinstance(self, JacobiConformal):

595 t += JacobiConformal.xyQ2,

596 t += T.volume, T.area

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

598

599 @Property_RO

600 def volume(self):

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

602 '''

603 a, b, c = self._abc3

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

605

606

607 class Triaxial(Triaxial_):

608 '''I{Ordered} triaxial ellipsoid.

609

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

611 '''

612 _unordered = False

613

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

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

616

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

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

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

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

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

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

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

624

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

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

627

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

629 '''

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

631

632 @Property_RO

633 def _a2b2_a2c2(self):

634 ''' @see: Method C{forwardBetaOmega}.

635 '''

636 return self._a2b2 / self._a2c2

637

638 @Property_RO

639 def area(self):

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

641

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

643 '''

644 a, b, c = self._abc3

645 if a != b:

646 kp2, k2 = self._k2_kp2 # swapped!

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

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

649 s2 = self. e2ac # sin(phi)**2 == 1 - c2

650 s = sqrt(s2)

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

652 b *= fsum_(aE.fE(r) * s, c / a * c / b,

653 aE.fF(r) * c2 / s, floats=True)

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

655 else: # a == b > c

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

657 return a

658

659 def _exyz3(self, u):

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

661 '''

662 if u > 0:

663 u2 = u**2

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

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

666 else:

667 x = y = u = _0_0

668 return x, y, u

669

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

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

672 and height to cartesian.

673

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

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

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

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

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

679

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

681

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

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

684 '''

685 if height:

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

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

688 else:

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

690 if z: # and x and y:

691 sa, ca = _SinCos2(beta)

692 sb, cb = _SinCos2(omega)

693

694 r = self._a2b2_a2c2

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

696 y *= ca * sb

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

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

699

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

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

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

703

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

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

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

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

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

709

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

711

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

713

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

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

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

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

718 '''

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

720 return Vector3Tuple(*t, name=name)

721

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

723 '''Project a cartesian on this triaxial.

724

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

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

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

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

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

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

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

732

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

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

735 '''

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

737 _ = t.rename(name)

738 return t

739

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

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

742

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

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

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

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

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

748

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

750

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

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

753 '''

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

755

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

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

758

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

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

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

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

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

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

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

766

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

768

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

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

771 '''

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

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

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

775

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

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

778 '''

779 ca_x_sb = ca * sb

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

781 t = fsum_(_1_0, -self.e2ac * sa**2, -self.e2ab * ca_x_sb**2, floats=True)

782 N = self.a / self._sqrt(t) # prime vertical

783 x = (h + N) * ca * cb

784 y = (h + N * self._1e2ab) * ca_x_sb

785 z = (h + N * self._1e2ac) * sa

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

787

788 @Property_RO

789 def _k2_kp2(self):

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

791 '''

792 # k2 = a2b2 / a2c2 * c2_b2

793 # kp2 = b2c2 / a2c2 * a2_b2

794 # b2 = b**2

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

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

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

798 return (self._a2b2 / self._a2c2 * self._c2_b2,

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

800

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

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

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

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

805 '''

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

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

808

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

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

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

812

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

814 z = Fsum(c2_a2, sb**2, b2_a2 * cb**2).fover(a2c2_a2)

815

816 x = self.a * cb * self._sqrt(x)

817 y = self.b * ca * sb

818 z = self.c * sa * self._sqrt(z)

819 return x, y, z

820

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

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

823 and height.

824

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

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

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

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

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

830

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

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

833 units as this triaxial's axes.

834

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

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

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

838 '''

839 v = _otherV3d_(x_xyz, y, z)

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

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

842

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

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

845

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

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

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

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

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

851 as the axes).

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

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

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

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

856

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

858

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

860

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

862 '''

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

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

865 if s: # PYCHOK no cover

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

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

868

869 if h:

870 if normal:

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

872 elif v.length > EPS0:

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

874 return v.xyz # Vector3Tuple

875

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

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

878

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

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

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

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

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

884

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

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

887 as this triaxial's axes.

888

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

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

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

892 '''

893 v = _otherV3d_(x_xyz, y, z)

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

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

896 _1_0)

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

898 return LatLon3Tuple(*t, name=name)

899

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

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

902 '''

903 x, y, z = s.xyz

904 d = hypot( x, y)

905 a = atan2_(z, d)

906 b = atan2_(y, x)

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

908 return a, b, h

909

910

911 class JacobiConformal(Triaxial):

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

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

914

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

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

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

918 returned in the case of an ellipsoid of revolution (or a sphere, I{currently

919 not supported}).

920

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

922 licensed under the MIT/X11 License.

923

924 @note: This constructor can not be used to specify a sphere.

925

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

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

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

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

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

931 '''

932# @Property_RO

933# def ab(self):

934# '''Get relative magnitude C{ab} (C{None} or C{meter}, same units as B{C{a}}).

935# '''

936# return self._ab

937

938# @Property_RO

939# def bc(self):

940# '''Get relative magnitude C{bc} (C{None} or C{meter}, same units as B{C{a}}).

941# '''

942# return self._bc

943

944 @Property_RO

945 def _xE(self):

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

947 '''

948 k2, kp2 = self._k2_kp2

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

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

951

952 def xR(self, omega):

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

954

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

956

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

958 '''

959 return self.xR_(*_SinCos2(omega))

960

961 def xR_(self, somega, comega):

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

963

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

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

966

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

968 '''

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

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

971

972 @Property_RO

973 def xyQ2(self):

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

975 '''

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

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

978 name=JacobiConformal.xyQ2.name)

979

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

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

982

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

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

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

986

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

988 '''

989 return self.xyR2_(*(_SinCos2(beta) + _SinCos2(omega)),

990 name=name or self.xyR2.__name__)

991

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

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

994

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

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

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

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

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

1000

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

1002 '''

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

1004 self.yR_(sbeta, cbeta),

1005 name=name or self.xyR2_.__name__)

1006

1007 @Property_RO

1008 def _yE(self):

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

1010 '''

1011 kp2, k2 = self._k2_kp2 # swapped!

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

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

1014

1015 def yR(self, beta):

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

1017

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

1019

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

1021 '''

1022 return self.yR_(*_SinCos2(beta))

1023

1024 def yR_(self, sbeta, cbeta):

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

1026

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

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

1029

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

1031 '''

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

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

1034

1035

1036 class TriaxialError(_ValueError):

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

1038 '''

1039 pass # ...

1040

1041

1042 class Triaxials(_NamedEnum):

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

1044 to accommodate the L{_LazyNamedEnumItem} properties.

1045 '''

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

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

1048 '''

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

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

1051

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

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

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

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

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

1057 Triaxials._assert( # a (km) b (km) c (km) planet

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

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

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

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

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

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

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

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

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

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

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

1069 WGS84_35 = _lazy('WGS84_35', *map(m2km, _WGS84_35abc)))

1070

1071 del _WGS84_35abc

1072

1073

1074 def _getitems(items, *indices):

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

1076

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

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

1079 the special method C{__getitem__}.

1080 '''

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

1082

1083

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

1085 '''(INTERNAL) Hartzell's "Satellite Lin-of-Sight Intersection ...", I{un- or ordered}.

1086 '''

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

1088

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

1090 b2, p2 = (b**2, (b / a)**2) if b != a else (a2, _1_0)

1091 c2, q2 = c**2, (c / a)**2

1092

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

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

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

1096

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

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

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

1100

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

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

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

1104 raise _ValueError(_near_(_null_))

1105

1106 r = fsum_(b2 * w2, c2 * v2, -v2 * z2, vy * wz * 2,

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

1108 -w2 * x2 * p2, -u2 * y2 * q2, -v2 * x2 * q2, ux * vy * 2 * q2, floats=True)

1109 if r > 0: # a2 factored out

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

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

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

1113

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

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

1116 s = fsum_(_1_0, x2 / a2, y2 / b2, z2 / c2, _N_2_0, floats=True) # like _sideOf

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

1118 elif fsum_(x2, y2, z2, floats=True) < d**2: # d past triaxial's center

1119 raise _ValueError(_too_(_distant_))

1120

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

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

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

1124

1125

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

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

1128 from a Point-Of-View outside.

1129

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

1131 or L{Vector3d}).

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

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

1134 @kwarg tri_biax: A triaxial (L{Triaxial}, L{Triaxial_}, L{JacobiConformal})

1135 or biaxial ellipsoid (L{Datum}, L{Ellipsoid}, L{Ellipsoid2},

1136 L{a_f2Tuple} or C{scalar} radius in C{meter}).

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

1138

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

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

1141

1142 @raise TriaxialError: Null B{C{pov}} or B{C{los}} vector, B{C{pov}} is inside

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

1144 tri-/biaxial or points in an opposite direction.

1145

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

1147

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

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

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

1151 '''

1152 if isinstance(tri_biax, Triaxial_):

1153 T = tri_biax

1154 else:

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

1156 _spherical_datum(tri_biax, name=hartzell4.__name__)

1157 E = D.ellipsoid

1158 T = Triaxial_(E.a, E.a, E.b, name=E.name)

1159

1160 try:

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

1162 except Exception as x:

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

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

1165

1166

1167 def _normalTo4(x, y, a, b, eps=EPS): # MCCABE 14

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

1169

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

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

1172 Documentation/DistancePointEllipseEllipsoid.pdf>}.

1173 '''

1174 def _root2d(r, u, v, g, eps):

1175 # robust root finder

1176 _1, __2 = _1_0, _0_5

1177 _a, _h2 = fabs, _hypot21_

1178 u *= r

1179 t0 = v - _1

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

1181 for i in range(1, _TRIPS):

1182 t = (t0 + t1) * __2

1183 if t in (t0, t1) or _a(t0 - t1) < eps:

1184 break

1185 g = _h2(u / (t + r), v / (t + _1))

1186 if g > 0:

1187 t0 = t

1188 elif g < 0:

1189 t1 = t

1190 else:

1191 break

1192 else: # PYCHOK no cover

1193 e = _a(t0 - t1)

1194 t = _root2d.__name__

1195 raise _ValueError(Fmt.no_convergence(e, eps), txt=t)

1196 return t, i

1197

1198 if a < b:

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

1200 return a, b, d, i

1201

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

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

1204

1205 i = None

1206 if y:

1207 if x:

1208 u = fabs(x / a)

1209 v = fabs(y / b)

1210 g = _hypot21_(u, v)

1211 if g:

1212 r = (a / b)**2

1213 t, i = _root2d(r, u, v, g, eps)

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

1215 b = y / (t + _1_0)

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

1217 else: # on the ellipse

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

1219 else: # x == 0

1220 if y < 0:

1221 b = -b

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

1223

1224 else: # y == 0

1225 n = a * x

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

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

1228 r = n / d

1229 a *= r

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

1231 d = hypot(x - a, b)

1232 else:

1233 if x < 0:

1234 a = -a

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

1236 return a, b, d, i

1237

1238

1239 def _normalTo5(x, y, z, Tun, eps=EPS): # MCCABE 24

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

1241

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

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

1244 '''

1245 def _root3d(r, s, u, v, w, g, eps):

1246 # robust root finder

1247 _1, __2 = _1_0, _0_5

1248 _a, _h2 = fabs, _hypot21_

1249 u *= r

1250 v *= s

1251 t0 = w - _1

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

1253 for i in range(1, _TRIPS):

1254 t = (t0 + t1) * __2

1255 if t in (t0, t1) or _a(t0 - t1) < eps:

1256 break

1257 g = _h2(u / (t + r), v / (t + s), w / (t + _1))

1258 if g > 0:

1259 t0 = t

1260 elif g < 0:

1261 t1 = t

1262 else:

1263 break

1264 else: # PYCHOK no cover

1265 e = _a(t0 - t1)

1266 t = _root3d.__name__

1267 raise _ValueError(Fmt.no_convergence(e, eps), txt=t)

1268 return t, i

1269

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

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

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

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

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

1275

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

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

1278

1279 if eps > 0:

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

1281 else: # no validation

1282 val, eps = 0, -eps

1283

1284 i = None

1285 if z:

1286 if y:

1287 if x:

1288 u = fabs(x / a)

1289 v = fabs(y / b)

1290 w = fabs(z / c)

1291 g = _hypot21_(u, v, w)

1292 if g:

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

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

1295 t, i = _root3d(_1_0 / r, _1_0 / s, u, v, w, g, eps)

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

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

1298 c = z / (t + _1_0)

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

1300 else: # on the ellipsoid

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

1302 else: # x == 0

1303 a = x # 0

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

1305 elif x: # y == 0

1306 b = y # 0

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

1308 else: # x == y == 0

1309 if z < 0:

1310 c = -c

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

1312

1313 else: # z == 0

1314 t = False

1315 n = a * x

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

1317 if d > fabs(n):

1318 u = n / d

1319 n = b * y

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

1321 if d > fabs(n):

1322 v = n / d

1323 n = _hypot21_(u, v)

1324 if n < 0:

1325 a *= u

1326 b *= v

1327 c *= sqrt(-n)

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

1329 t = True

1330 if not t:

1331 c = z # 0

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

1333

1334 if val > 0: # validate

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

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

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

1338 sideOf=e, eps=val)

1339 if d: # angle of delta and normal vector

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

1341 if m.euclid > val:

1342 m = m.unit()

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

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

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

1346 raise _ValueError(n=n, m=m,

1347 dot=e, eps=val)

1348 return a, b, c, d, i

1349

1350

1351 def _otherV3d_(x_xyz, y, z, name=NN):

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

1353 '''

1354 return Vector3d(x_xyz, y, z, name=name) if isscalar(x_xyz) else \

1355 _otherV3d(x_xyz=x_xyz)

1356

1357

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

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

1360

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

1362 '''

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

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

1365

1366

1367 def _SinCos2(x):

1368 '''Get C{sin} and C{cos} of C{x} in C{Degrees}, C{Radians} or {radians}.

1369 '''

1370 return sincos2d(x) if isinstance(x, Degrees) else (

1371 sincos2(x) if isinstance(x, Radians) else

1372 sincos2(float(x))) # assume C{radians}

1373

1374

1375 if __name__ == '__main__':

1376

1377 from pygeodesy import printf

1378

1379 # __doc__ of this file, force all into registery

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

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

1382

1383# **) MIT License

1384#

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

1386#

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

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

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

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

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

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

1393#

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

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

1396#

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

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

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

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

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

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

1403# OTHER DEALINGS IN THE SOFTWARE.