Coverage for pygeodesy/triaxials.py: 95%
553 statements
« prev ^ index » next coverage.py v7.2.2, created at 2023-04-23 16:38 -0400
« prev ^ index » next coverage.py v7.2.2, created at 2023-04-23 16:38 -0400
2# -*- coding: utf-8 -*-
4u'''Triaxal ellipsoid classes L{JacobiConformal}, Jacobi's conformal projection, trancoded
5from I{Charles Karney}'s C++ class U{JacobiConformal<https://GeographicLib.SourceForge.io/C++/
6doc/classGeographicLib_1_1JacobiConformal.html#details>} to pure Python, I{ordered} L{Triaxial}
7and I{unordered} L{Triaxial_} and miscellaneous classes L{BetaOmega2Tuple}, L{BetaOmega3Tuple},
8L{Jacobi2Tuple} and L{TriaxialError}.
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>}.
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
28from __future__ import division as _; del _ # PYCHOK semicolon
30from pygeodesy.basics import isscalar, map1, _ValueError, _zip
31from pygeodesy.constants import EPS, EPS0, EPS02, EPS4, _EPS2e4, INT0, PI2, PI_3, PI4, \
32 _0_0, _0_5, _1_0, _N_2_0, float0_, isfinite, isnear1, \
33 _4_0 # PYCHOK used!
34from pygeodesy.datums import Datum, Ellipsoid, Fmt, _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 .basics
38from pygeodesy.fmath import Fdot, fdot, fmean_, hypot, hypot_, _hypot21_, norm2
39from pygeodesy.fsums import Fsum, fsum_, Property_RO
40from 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
45from pygeodesy.named import _NamedBase, _NamedEnum, _NamedEnumItem, \
46 _NamedTuple, _Pass, _lazyNamedEnumItem as _lazy
47from pygeodesy.namedTuples import LatLon3Tuple, Vector3Tuple, Vector4Tuple
48# from pygeodesy.props import Property_RO # from .fsums
49# from pygeodesy.streprs import Fmt # from .datums
50from pygeodesy.units import Degrees, Float, Height_, Meter, Meter2, Meter3, \
51 Radians, Radius
52from pygeodesy.utily import asin1, atan2d, km2m, m2km, SinCos2, sincos2d_
53from pygeodesy.vector3d import _ALL_LAZY, _MODS, _otherV3d, Vector3d
55from math import atan2, fabs, sqrt
57__all__ = _ALL_LAZY.triaxials
58__version__ = '23.04.14'
60_not_ordered_ = _not_('ordered')
61_omega_ = 'omega'
62_TRIPS = 537 # max 55, Eberly 1074?
65class _ToNamedBase(_NamedBase):
66 '''(INTERNAL) C{-.toDegrees}, C{-.toRadians} base.
67 '''
68 def _toDegrees(self, a, b, *c, **toDMS_kwds):
69 if toDMS_kwds:
70 toDMS = _MODS.dms.toDMS
71 a = toDMS(a.toDegrees(), **toDMS_kwds)
72 b = toDMS(b.toDegrees(), **toDMS_kwds)
73 elif isinstance(a, Degrees) and \
74 isinstance(b, Degrees):
75 return self
76 else:
77 a, b = a.toDegrees(), b.toDegrees()
78 return self.classof(a, b, *c, name=self.name)
80 def _toRadians(self, a, b, *c):
81 return self if isinstance(a, Radians) and \
82 isinstance(b, Radians) else \
83 self.classof(a.toRadians(), b.toRadians(),
84 *c, name=self.name)
87class BetaOmega2Tuple(_NamedTuple, _ToNamedBase):
88 '''2-Tuple C{(beta, omega)} with I{ellipsoidal} lat- and
89 longitude C{beta} and C{omega} both in C{Radians} (or
90 C{Degrees}).
91 '''
92 _Names_ = (_beta_, _omega_)
93 _Units_ = (_Pass, _Pass)
95 def toDegrees(self, **toDMS_kwds):
96 '''Convert this L{BetaOmega2Tuple} to C{Degrees} or C{toDMS}.
98 @return: L{BetaOmega2Tuple}C{(beta, omega)} with
99 C{beta} and C{omega} both in C{Degrees}
100 or as an L{toDMS} string provided some
101 B{C{toDMS_kwds}} are supplied.
102 '''
103 return _ToNamedBase._toDegrees(self, *self, **toDMS_kwds)
105 def toRadians(self):
106 '''Convert this L{BetaOmega2Tuple} to C{Radians}.
108 @return: L{BetaOmega2Tuple}C{(beta, omega)} with
109 C{beta} and C{omega} both in C{Radians}.
110 '''
111 return _ToNamedBase._toRadians(self, *self)
114class BetaOmega3Tuple(_NamedTuple, _ToNamedBase):
115 '''3-Tuple C{(beta, omega, height)} with I{ellipsoidal} lat- and
116 longitude C{beta} and C{omega} both in C{Radians} (or C{Degrees})
117 and the C{height}, rather the (signed) I{distance} to the triaxial's
118 surface (measured along the radial line to the triaxial's center)
119 in C{meter}, conventionally.
120 '''
121 _Names_ = BetaOmega2Tuple._Names_ + (_height_,)
122 _Units_ = BetaOmega2Tuple._Units_ + ( Meter,)
124 def toDegrees(self, **toDMS_kwds):
125 '''Convert this L{BetaOmega3Tuple} to C{Degrees} or C{toDMS}.
127 @return: L{BetaOmega3Tuple}C{(beta, omega, height)} with
128 C{beta} and C{omega} both in C{Degrees} or as an
129 L{toDMS} string provided some B{C{toDMS_kwds}}
130 are supplied.
131 '''
132 return _ToNamedBase._toDegrees(self, *self, **toDMS_kwds)
134 def toRadians(self):
135 '''Convert this L{BetaOmega3Tuple} to C{Radians}.
137 @return: L{BetaOmega3Tuple}C{(beta, omega, height)} with
138 C{beta} and C{omega} both in C{Radians}.
139 '''
140 return _ToNamedBase._toRadians(self, *self)
142 def to2Tuple(self):
143 '''Reduce this L{BetaOmega3Tuple} to a L{BetaOmega2Tuple}.
144 '''
145 return BetaOmega2Tuple(*self[:2])
148class Jacobi2Tuple(_NamedTuple, _ToNamedBase):
149 '''2-Tuple C{(x, y)} with a Jacobi Conformal C{x} and C{y}
150 projection, both in C{Radians} (or C{Degrees}).
151 '''
152 _Names_ = (_x_, _y_)
153 _Units_ = (_Pass, _Pass)
155 def toDegrees(self, **toDMS_kwds):
156 '''Convert this L{Jacobi2Tuple} to C{Degrees} or C{toDMS}.
158 @return: L{Jacobi2Tuple}C{(x, y)} with C{x} and C{y}
159 both in C{Degrees} or as an L{toDMS} string
160 provided some B{C{toDMS_kwds}} are supplied.
161 '''
162 return _ToNamedBase._toDegrees(self, *self, **toDMS_kwds)
164 def toRadians(self):
165 '''Convert this L{Jacobi2Tuple} to C{Radians}.
167 @return: L{Jacobi2Tuple}C{(x, y)} with C{x}
168 and C{y} both in C{Radians}.
169 '''
170 return _ToNamedBase._toRadians(self, *self)
173class Triaxial_(_NamedEnumItem):
174 '''I{Unordered} triaxial ellipsoid and base class.
176 Triaxial ellipsoids with right-handed semi-axes C{a}, C{b} and C{c}, oriented
177 such that the large principal ellipse C{ab} is the equator I{Z}=0, I{beta}=0,
178 while the small principal ellipse C{ac} is the prime meridian, plane I{Y}=0,
179 I{omega}=0.
181 The four umbilic points, C{abs}(I{omega}) = C{abs}(I{beta}) = C{PI/2}, lie on
182 the middle principal ellipse C{bc} in plane I{X}=0, I{omega}=C{PI/2}.
184 @note: I{Geodetic} C{lat}- and C{lon}gitudes are in C{degrees}, I{geodetic}
185 C{phi} and C{lam}bda are in C{radians}, but I{ellipsoidal} lat- and
186 longitude C{beta} and C{omega} are in C{Radians} by default (or in
187 C{Degrees} if converted).
188 '''
189 _ijk = _kji = None
190 _unordered = True
192 def __init__(self, a_triaxial, b=None, c=None, name=NN):
193 '''New I{unordered} L{Triaxial_}.
195 @arg a_triaxial: C{X} semi-axis (C{scalar}, conventionally in C{meter})
196 or an other L{Triaxial} or L{Triaxial_} instance.
197 @kwarg b: C{Y} semi-axis (C{meter}, same units as B{C{a}}), required
198 if C{B{a_triaxial} is scalar}, ignored otherwise.
199 @kwarg c: C{Z} semi-axis (C{meter}, same units as B{C{a}}), required
200 if C{B{a_triaxial} is scalar}, ignored otherwise.
201 @kwarg name: Optional name (C{str}).
203 @raise TriaxialError: Invalid semi-axis or -axes.
204 '''
205 try:
206 a = a_triaxial
207 t = a._abc3 if isinstance(a, Triaxial_) else (
208 Radius(a=a), Radius(b=b), Radius(c=c))
209 except (TypeError, ValueError) as x:
210 raise TriaxialError(a=a, b=b, c=c, cause=x)
211 if name:
212 self.name = name
214 a, b, c = self._abc3 = t
215 if self._unordered: # == not isinstance(self, Triaxial)
216 s, _, t = sorted(t)
217 if not (isfinite(t) and s > 0):
218 raise TriaxialError(a=a, b=b, c=c) # txt=_invalid_
219 elif not (isfinite(a) and a >= b >= c > 0):
220 raise TriaxialError(a=a, b=b, c=c, txt=_not_ordered_)
221 elif not (a > c and self._a2c2 > 0 and self.e2ac > 0):
222 raise TriaxialError(a=a, c=c, e2ac=self.e2ac, txt=_spherical_)
224 def __str__(self):
225 return self.toStr()
227 @Property_RO
228 def a(self):
229 '''Get the (largest) C{x} semi-axis (C{meter}, conventionally).
230 '''
231 a, _, _ = self._abc3
232 return a
234 @Property_RO
235 def _a2b2(self):
236 '''(INTERNAL) Get C{a**2 - b**2} == E_sub_e**2.
237 '''
238 a, b, _ = self._abc3
239 return ((a - b) * (a + b)) if a != b else _0_0
241 @Property_RO
242 def _a2_b2(self):
243 '''(INTERNAL) Get C{(a/b)**2}.
244 '''
245 a, b, _ = self._abc3
246 return (a / b)**2 if a != b else _1_0
248 @Property_RO
249 def _a2c2(self):
250 '''(INTERNAL) Get C{a**2 - c**2} == E_sub_x**2.
251 '''
252 a, _, c = self._abc3
253 return ((a - c) * (a + c)) if a != c else _0_0
255 @Property_RO
256 def area(self):
257 '''Get the surface area (C{meter} I{squared}).
258 '''
259 c, b, a = sorted(self._abc3)
260 if a > c:
261 a = Triaxial(a, b, c).area if a > b else \
262 Ellipsoid(a, b=c).areax # a == b
263 else: # a == c == b
264 a = Meter2(area=a**2 * PI4)
265 return a
267 def area_p(self, p=1.6075):
268 '''I{Approximate} the surface area (C{meter} I{squared}).
270 @kwarg p: Exponent (C{scalar} > 0), 1.6 for near-spherical or 1.5849625007
271 for "near-flat" triaxials.
273 @see: U{Surface area<https://WikiPedia.org/wiki/Ellipsoid#Approximate_formula>}.
274 '''
275 a, b, c = self._abc3
276 if a == b == c:
277 a *= a
278 else:
279 _p = pow
280 a = _p(fmean_(_p(a * b, p), _p(a * c, p), _p(b * c, p)), _1_0 / p)
281 return Meter2(area_p=a * PI4)
283 @Property_RO
284 def b(self):
285 '''Get the (middle) C{y} semi-axis (C{meter}, same units as B{C{a}}).
286 '''
287 _, b, _ = self._abc3
288 return b
290 @Property_RO
291 def _b2c2(self):
292 '''(INTERNAL) Get C{b**2 - c**2} == E_sub_y**2.
293 '''
294 _, b, c = self._abc3
295 return ((b - c) * (b + c)) if b != c else _0_0
297 @Property_RO
298 def c(self):
299 '''Get the (smallest) C{z} semi-axis (C{meter}, same units as B{C{a}}).
300 '''
301 _, _, c = self._abc3
302 return c
304 @Property_RO
305 def _c2_b2(self):
306 '''(INTERNAL) Get C{(c/b)**2}.
307 '''
308 _, b, c = self._abc3
309 return (c / b)**2 if b != c else _1_0
311 @Property_RO
312 def e2ab(self):
313 '''Get the C{ab} ellipse' I{(1st) eccentricity squared} (C{scalar}), M{1 - (b/a)**2}.
314 '''
315 return Float(e2ab=(_1_0 - self._1e2ab) or _0_0)
317 @Property_RO
318 def _1e2ab(self):
319 '''(INTERNAL) Get C{1 - e2ab} == C{(b/a)**2}.
320 '''
321 a, b, _ = self._abc3
322 return (b / a)**2 if a != b else _1_0
324 @Property_RO
325 def e2bc(self):
326 '''Get the C{bc} ellipse' I{(1st) eccentricity squared} (C{scalar}), M{1 - (c/b)**2}.
327 '''
328 return Float(e2bc=(_1_0 - self._1e2bc) or _0_0)
330 _1e2bc = _c2_b2 # C{1 - e2bc} == C{(c/b)**2}
332 @Property_RO
333 def e2ac(self):
334 '''Get the C{ac} ellipse' I{(1st) eccentricity squared} (C{scalar}), M{1 - (c/a)**2}.
335 '''
336 return Float(e2ac=(_1_0 - self._1e2ac) or _0_0)
338 @Property_RO
339 def _1e2ac(self):
340 '''(INTERNAL) Get C{1 - e2ac} == C{(c/a)**2}.
341 '''
342 a, _, c = self._abc3
343 return (c / a)**2 if a != c else _1_0
345 @Property_RO
346 def _Elliptic(self):
347 '''(INTERNAL) Get class L{Elliptic} once.
348 '''
349 return _MODS.elliptic.Elliptic
351 def hartzell4(self, pov, los=None, name=NN):
352 '''Compute the intersection of this triaxial's surface with a Line-Of-Sight
353 from a Point-Of-View in space.
355 @see: Function L{pygeodesy.hartzell4} for further details.
356 '''
357 return hartzell4(pov, los=los, tri_biax=self, name=name)
359 def height4(self, x_xyz, y=None, z=None, normal=True, eps=EPS):
360 '''Compute the projection on and the height of a cartesian above or below
361 this triaxial's surface.
363 @arg x_xyz: X component (C{scalar}) or a cartesian (C{Cartesian},
364 L{Ecef9Tuple}, L{Vector3d}, L{Vector3Tuple} or L{Vector4Tuple}).
365 @kwarg y: Y component (C{scalar}), required if B{C{x_xyz}} if C{scalar}.
366 @kwarg z: Z component (C{scalar}), required if B{C{x_xyz}} if C{scalar}.
367 @kwarg normal: If C{True} the projection is perpendicular to (the nearest
368 point on) this triaxial's surface, otherwise the C{radial}
369 line to this triaxial's center (C{bool}).
370 @kwarg eps: Tolerance for root finding and validation (C{scalar}), use a
371 negative value to skip validation.
373 @return: L{Vector4Tuple}C{(x, y, z, h)} with the cartesian coordinates
374 C{x}, C{y} and C{z} of the projection on or the intersection
375 with and with the height C{h} above or below the triaxial's
376 surface in C{meter}, conventionally.
378 @raise TriaxialError: Non-cartesian B{C{xyz}}, invalid B{C{eps}}, no
379 convergence in root finding or validation failed.
381 @see: Method L{Ellipsoid.height4} and I{Eberly}'s U{Distance from a Point
382 to ... an Ellipsoid ...<https://www.GeometricTools.com/Documentation/
383 DistancePointEllipseEllipsoid.pdf>}.
384 '''
385 v, r = _otherV3d_(x_xyz, y, z), self.isSpherical
387 i, h = None, v.length
388 if h < EPS0: # EPS
389 x = y = z = _0_0
390 h -= min(self._abc3) # nearest
391 elif r: # .isSpherical
392 x, y, z = v.times(r / h).xyz
393 h -= r
394 else:
395 x, y, z = v.xyz
396 try:
397 if normal: # perpendicular to triaxial
398 x, y, z, h, i = _normalTo5(x, y, z, self, eps=eps)
399 else: # radially to triaxial's center
400 x, y, z = self._radialTo3(z, hypot(x, y), y, x)
401 h = v.minus_(x, y, z).length
402 except Exception as e:
403 raise TriaxialError(x=x, y=y, z=z, cause=e)
404 if h > 0 and self.sideOf(v, eps=EPS0) < 0:
405 h = -h # below the surface
406 return Vector4Tuple(x, y, z, h, iteration=i, name=self.height4.__name__)
408 @Property_RO
409 def isOrdered(self):
410 '''Is this triaxial I{ordered} and I{not spherical} (C{bool})?
411 '''
412 a, b, c = self._abc3
413 return bool(a >= b > c) # b > c!
415 @Property_RO
416 def isSpherical(self):
417 '''Is this triaxial I{spherical} (C{Radius} or INT0)?
418 '''
419 a, b, c = self._abc3
420 return a if a == b == c else INT0
422 def normal3d(self, x_xyz, y=None, z=None, length=_1_0):
423 '''Get a 3-D vector perpendicular to at a cartesian on this triaxial's surface.
425 @arg x_xyz: X component (C{scalar}) or a cartesian (C{Cartesian},
426 L{Ecef9Tuple}, L{Vector3d}, L{Vector3Tuple} or L{Vector4Tuple}).
427 @kwarg y: Y component (C{scalar}), required if B{C{x_xyz}} if C{scalar}.
428 @kwarg z: Z component (C{scalar}), required if B{C{x_xyz}} if C{scalar}.
429 @kwarg length: Optional length and in-/outward direction (C{scalar}).
431 @return: A C{Vector3d(x_, y_, z_)} normalized to B{C{length}}, pointing
432 in- or outward for neg- respectively positive B{C{length}}.
434 @note: Cartesian location C{(B{x}, B{y}, B{z})} must be on this triaxial's
435 surface, use method L{Triaxial.sideOf} to validate.
436 '''
437 # n = 2 * (x / a2, y / b2, z / c2)
438 # == 2 * (x, y * a2 / b2, z * a2 / c2) / a2 # iff ordered
439 # == 2 * (x, y / _1e2ab, z / _1e2ac) / a2
440 # == unit(x, y / _1e2ab, z / _1e2ac).times(length)
441 n = self._normal3d.times_(*_otherV3d_(x_xyz, y, z).xyz)
442 if n.length < EPS0:
443 raise TriaxialError(x=x_xyz, y=y, z=z, txt=_null_)
444 return n.times(length / n.length)
446 @Property_RO
447 def _normal3d(self):
448 '''(INTERNAL) Get M{Vector3d((d/a)**2, (d/b)**2, (d/c)**2)}, M{d = max(a, b, c)}.
449 '''
450 d = max(self._abc3)
451 t = tuple(((d / x)**2 if x != d else _1_0) for x in self._abc3)
452 return Vector3d(*t, name=self.normal3d.__name__)
454 def _norm2(self, s, c, *a):
455 '''(INTERNAL) Normalize C{s} and C{c} iff not already.
456 '''
457 if fabs(s) > _1_0 or fabs(c) > _1_0 or \
458 fabs(_hypot21_(s, c)) > EPS0:
459 s, c = norm2(s, c)
460 if a:
461 s, c = norm2(s * self.b, c * a[0])
462 return float0_(s, c)
464 def _order3(self, *abc, **reverse): # reverse=False
465 '''(INTERNAL) Un-/Order C{a}, C{b} and C{c}.
467 @return: 3-Tuple C{(a, b, c)} ordered by or un-ordered
468 (reverse-ordered) C{ijk} if C{B{reverse}=True}.
469 '''
470 ijk = self._order_ijk(**reverse)
471 return _getitems(abc, *ijk) if ijk else abc
473 def _order3d(self, v, **reverse): # reverse=False
474 '''(INTERNAL) Un-/Order a C{Vector3d}.
476 @return: Vector3d(x, y, z) un-/ordered.
477 '''
478 ijk = self._order_ijk(**reverse)
479 return v.classof(*_getitems(v.xyz, *ijk)) if ijk else v
481 @Property_RO
482 def _ordered4(self):
483 '''(INTERNAL) Helper for C{_hartzell3d2} and C{_normalTo5}.
484 '''
485 def _order2(reverse, a, b, c):
486 '''(INTERNAL) Un-Order C{a}, C{b} and C{c}.
488 @return: 2-Tuple C{((a, b, c), ijk)} with C{a} >= C{b} >= C{c}
489 and C{ijk} a 3-tuple with the initial indices.
490 '''
491 i, j, k = 0, 1, 2 # range(3)
492 if a < b:
493 a, b, i, j = b, a, j, i
494 if a < c:
495 a, c, i, k = c, a, k, i
496 if b < c:
497 b, c, j, k = c, b, k, j
498 # reverse (k, j, i) since (a, b, c) is reversed-sorted
499 ijk = (k, j, i) if reverse else (None if i < j < k else (i, j, k))
500 return (a, b, c), ijk
502 abc, T = self._abc3, self
503 if not self.isOrdered:
504 abc, ijk = _order2(False, *abc)
505 if ijk:
506 _, kji = _order2(True, *ijk)
507 T = Triaxial_(*abc)
508 T._ijk, T._kji = ijk, kji
509 return abc + (T,)
511 def _order_ijk(self, reverse=False):
512 '''(INTERNAL) Get the un-/order indices.
513 '''
514 return self._kji if reverse else self._ijk
516 def _radialTo3(self, sbeta, cbeta, somega, comega):
517 '''(INTERNAL) I{Unordered} helper for C{.height4}.
518 '''
519 def _rphi(a, b, sphi, cphi):
520 # <https://WikiPedia.org/wiki/Ellipse#Polar_form_relative_to_focus>
521 # polar form: radius(phi) = a * b / hypot(a * sphi, b * cphi)
522 return (b / hypot(sphi, b / a * cphi)) if a > b else (
523 (a / hypot(cphi, a / b * sphi)) if a < b else a)
525 sa, ca = self._norm2(sbeta, cbeta)
526 sb, cb = self._norm2(somega, comega)
528 a, b, c = self._abc3
529 if a != b:
530 a = _rphi(a, b, sb, cb)
531 if a != c:
532 c = _rphi(a, c, sa, ca)
533 z, r = c * sa, c * ca
534 x, y = r * cb, r * sb
535 return x, y, z
537 def sideOf(self, x_xyz, y=None, z=None, eps=EPS4):
538 '''Is a cartesian above, below or on the surface of this triaxial?
540 @arg x_xyz: X component (C{scalar}) or a cartesian (C{Cartesian},
541 L{Ecef9Tuple}, L{Vector3d}, L{Vector3Tuple} or L{Vector4Tuple}).
542 @kwarg y: Y component (C{scalar}), required if B{C{x_xyz}} if C{scalar}.
543 @kwarg z: Z component (C{scalar}), required if B{C{x_xyz}} if C{scalar}.
544 @kwarg eps: Near surface tolerance(C{scalar}).
546 @return: C{INT0} if C{(B{x}, B{y}, B{z})} is near this triaxial's surface
547 within tolerance B{C{eps}}, otherwise a neg- or positive C{float}
548 if in- respectively outside this triaxial.
550 @see: Methods L{Triaxial.height4} and L{Triaxial.normal3d}.
551 '''
552 return _sideOf(_otherV3d_(x_xyz, y, z).xyz, self._abc3, eps=eps)
554 def _sqrt(self, x):
555 '''(INTERNAL) Helper.
556 '''
557 if x < 0:
558 raise TriaxialError(Fmt.PAREN(sqrt=x))
559 return _0_0 if x < EPS02 else sqrt(x)
561 def toEllipsoid(self, name=NN):
562 '''Convert this triaxial to an L{Ellipsoid}, provided C{a == b} or C{b == c}.
564 @return: An L{Ellipsoid} with north along this C{Z} axis if C{a == b},
565 this C{Y} axis if C{a == c} or this C{X} axis if C{b == c}.
567 @raise TriaxialError: This C{a != b}, C{b != c} and C{c != a}.
569 @see: Method L{Ellipsoid.toTriaxial}.
570 '''
571 a, b, c = self._abc3
572 if a == b: # N = Z
573 b = c
574 elif b == c: # N = X
575 a, b = b, a
576 elif a != c:
577 t = _SPACE_(_a_, _NOTEQUAL_, _b_, _NOTEQUAL_, _c_)
578 raise TriaxialError(a=a, b=b, c=c, txt=t)
579 return Ellipsoid(a, b=b, name=name or self.name)
581 def toStr(self, prec=9, name=NN, **unused): # PYCHOK signature
582 '''Return this C{Triaxial} as a string.
584 @kwarg prec: Precision, number of decimal digits (0..9).
585 @kwarg name: Override name (C{str}) or C{None} to exclude
586 this triaxial's name.
588 @return: This C{Triaxial}'s attributes (C{str}).
589 '''
590 T = Triaxial_
591 t = T.a, T.b, T.c, T.e2ab, T.e2bc, T.e2ac
592 if isinstance(self, JacobiConformal):
593 t += JacobiConformal.xyQ2,
594 t += T.volume, T.area
595 return self._instr(name, prec, props=t, area_p=self.area_p())
597 @Property_RO
598 def volume(self):
599 '''Get the volume (C{meter**3}), M{4 / 3 * PI * a * b * c}.
600 '''
601 a, b, c = self._abc3
602 return Meter3(volume=a * b * c * PI_3 * _4_0)
605class Triaxial(Triaxial_):
606 '''I{Ordered} triaxial ellipsoid.
608 @see: L{Triaxial_} for more information.
609 '''
610 _unordered = False
612 def __init__(self, a_triaxial, b=None, c=None, name=NN):
613 '''New I{ordered} L{Triaxial}.
615 @arg a_triaxial: Largest semi-axis (C{scalar}, conventionally in C{meter})
616 or an other L{Triaxial} or L{Triaxial_} instance.
617 @kwarg b: Middle semi-axis (C{meter}, same units as B{C{a}}), required
618 if C{B{a_triaxial} is scalar}, ignored otherwise.
619 @kwarg c: Smallest semi-axis (C{meter}, same units as B{C{a}}), required
620 if C{B{a_triaxial} is scalar}, ignored otherwise.
621 @kwarg name: Optional name (C{str}).
623 @note: The semi-axes must be ordered as C{B{a} >= B{b} >= B{c} > 0} and
624 must be ellipsoidal, C{B{a} > B{c}}.
626 @raise TriaxialError: Semi-axes not ordered, spherical or invalid.
627 '''
628 Triaxial_.__init__(self, a_triaxial, b=b, c=c, name=name)
630 @Property_RO
631 def _a2b2_a2c2(self):
632 '''@see: Methods C{.forwardBetaOmega} and C{._k2_kp2}.
633 '''
634 return self._a2b2 / self._a2c2
636 @Property_RO
637 def area(self):
638 '''Get the surface area (C{meter} I{squared}).
640 @see: U{Surface area<https://WikiPedia.org/wiki/Ellipsoid#Surface_area>}.
641 '''
642 a, b, c = self._abc3
643 if a != b:
644 kp2, k2 = self._k2_kp2 # swapped!
645 aE = self._Elliptic(k2, _0_0, kp2, _1_0)
646 c2 = self._1e2ac # cos(phi)**2 == (c/a)**2
647 s2 = self. e2ac # sin(phi)**2 == 1 - c2
648 s = sqrt(s2)
649 r = asin1(s) # phi == atan2(sqrt(c2), s)
650 b *= fsum_(aE.fE(r) * s, c / a * c / b,
651 aE.fF(r) * c2 / s, floats=True)
652 a = Meter2(area=a * b * PI2)
653 else: # a == b > c
654 a = Ellipsoid(a, b=c).areax
655 return a
657 def _exyz3(self, u):
658 '''(INTERNAL) Helper for C{.forwardBetOmg}.
659 '''
660 if u > 0:
661 u2 = u**2
662 x = u * self._sqrt(_1_0 + self._a2c2 / u2)
663 y = u * self._sqrt(_1_0 + self._b2c2 / u2)
664 else:
665 x = y = u = _0_0
666 return x, y, u
668 def forwardBetaOmega(self, beta, omega, height=0, name=NN):
669 '''Convert I{ellipsoidal} lat- and longitude C{beta}, C{omega}
670 and height to cartesian.
672 @arg beta: Ellipsoidal latitude (C{radians} or L{Degrees}).
673 @arg omega: Ellipsoidal longitude (C{radians} or L{Degrees}).
674 @arg height: Height above or below the ellipsoid's surface (C{meter}, same
675 units as this triaxial's C{a}, C{b} and C{c} semi-axes).
676 @kwarg name: Optional name (C{str}).
678 @return: A L{Vector3Tuple}C{(x, y, z)}.
680 @see: Method L{Triaxial.reverseBetaOmega} and U{Expressions (23-25)<https://
681 www.Topo.Auth.GR/wp-content/uploads/sites/111/2021/12/09_Panou.pdf>}.
682 '''
683 if height:
684 h = Height_(height=height, low=-self.c, Error=TriaxialError)
685 x, y, z = self._exyz3(h + self.c)
686 else:
687 x, y, z = self._abc3 # == self._exyz3(self.c)
688 if z: # and x and y:
689 sa, ca = SinCos2(beta)
690 sb, cb = SinCos2(omega)
692 r = self._a2b2_a2c2
693 x *= cb * self._sqrt(ca**2 + r * sa**2)
694 y *= ca * sb
695 z *= sa * self._sqrt(_1_0 - r * cb**2)
696 return Vector3Tuple(x, y, z, name=name)
698 def forwardBetaOmega_(self, sbeta, cbeta, somega, comega, name=NN):
699 '''Convert I{ellipsoidal} lat- and longitude C{beta} and C{omega}
700 to cartesian coordinates I{on the triaxial's surface}.
702 @arg sbeta: Ellipsoidal latitude C{sin(beta)} (C{scalar}).
703 @arg cbeta: Ellipsoidal latitude C{cos(beta)} (C{scalar}).
704 @arg somega: Ellipsoidal longitude C{sin(omega)} (C{scalar}).
705 @arg comega: Ellipsoidal longitude C{cos(omega)} (C{scalar}).
706 @kwarg name: Optional name (C{str}).
708 @return: A L{Vector3Tuple}C{(x, y, z)} on the surface.
710 @raise TriaxialError: This triaxial is near-spherical.
712 @see: Method L{Triaxial.reverseBetaOmega}, U{Triaxial ellipsoid coordinate
713 system<https://WikiPedia.org/wiki/Geodesics_on_an_ellipsoid#
714 Triaxial_ellipsoid_coordinate_system>} and U{expressions (23-25)<https://
715 www.Topo.Auth.GR/wp-content/uploads/sites/111/2021/12/09_Panou.pdf>}.
716 '''
717 t = self._radialTo3(sbeta, cbeta, somega, comega)
718 return Vector3Tuple(*t, name=name)
720 def forwardCartesian(self, x_xyz, y=None, z=None, name=NN, **normal_eps):
721 '''Project a cartesian on this triaxial.
723 @arg x_xyz: X component (C{scalar}) or a cartesian (C{Cartesian},
724 L{Ecef9Tuple}, L{Vector3d}, L{Vector3Tuple} or L{Vector4Tuple}).
725 @kwarg y: Y component (C{scalar}), required if B{C{x_xyz}} if C{scalar}.
726 @kwarg z: Z component (C{scalar}), required if B{C{x_xyz}} if C{scalar}.
727 @kwarg name: Optional name (C{str}).
728 @kwarg normal_eps: Optional keyword arguments C{B{normal}=True} and
729 C{B{eps}=EPS}, see method L{Triaxial.height4}.
731 @see: Method L{Triaxial.height4} for further information and method
732 L{Triaxial.reverseCartesian} to reverse the projection.
733 '''
734 t = self.height4(x_xyz, y, z, **normal_eps)
735 _ = t.rename(name)
736 return t
738 def forwardLatLon(self, lat, lon, height=0, name=NN):
739 '''Convert I{geodetic} lat-, longitude and heigth to cartesian.
741 @arg lat: Geodetic latitude (C{degrees}).
742 @arg lon: Geodetic longitude (C{degrees}).
743 @arg height: Height above the ellipsoid (C{meter}, same units
744 as this triaxial's C{a}, C{b} and C{c} axes).
745 @kwarg name: Optional name (C{str}).
747 @return: A L{Vector3Tuple}C{(x, y, z)}.
749 @see: Method L{Triaxial.reverseLatLon} and U{Expressions (9-11)<https://
750 www.Topo.Auth.GR/wp-content/uploads/sites/111/2021/12/09_Panou.pdf>}.
751 '''
752 return self._forwardLatLon3(height, name, *sincos2d_(lat, lon))
754 def forwardLatLon_(self, slat, clat, slon, clon, height=0, name=NN):
755 '''Convert I{geodetic} lat-, longitude and heigth to cartesian.
757 @arg slat: Geodetic latitude C{sin(lat)} (C{scalar}).
758 @arg clat: Geodetic latitude C{cos(lat)} (C{scalar}).
759 @arg slon: Geodetic longitude C{sin(lon)} (C{scalar}).
760 @arg clon: Geodetic longitude C{cos(lon)} (C{scalar}).
761 @arg height: Height above the ellipsoid (C{meter}, same units
762 as this triaxial's axes C{a}, C{b} and C{c}).
763 @kwarg name: Optional name (C{str}).
765 @return: A L{Vector3Tuple}C{(x, y, z)}.
767 @see: Method L{Triaxial.reverseLatLon} and U{Expressions (9-11)<https://
768 www.Topo.Auth.GR/wp-content/uploads/sites/111/2021/12/09_Panou.pdf>}.
769 '''
770 sa, ca = self._norm2(slat, clat)
771 sb, cb = self._norm2(slon, clon)
772 return self._forwardLatLon3(height, name, sa, ca, sb, cb)
774 def _forwardLatLon3(self, h, name, sa, ca, sb, cb):
775 '''(INTERNAL) Helper for C{.forwardLatLon} and C{.forwardLatLon_}.
776 '''
777 ca_x_sb = ca * sb
778 # 1 - (1 - (c/a)**2) * sa**2 - (1 - (b/a)**2) * ca**2 * sb**2
779 t = fsum_(_1_0, -self.e2ac * sa**2, -self.e2ab * ca_x_sb**2, floats=True)
780 n = self.a / self._sqrt(t) # prime vertical
781 x = (h + n) * ca * cb
782 y = (h + n * self._1e2ab) * ca_x_sb
783 z = (h + n * self._1e2ac) * sa
784 return Vector3Tuple(x, y, z, name=name)
786 @Property_RO
787 def _k2_kp2(self):
788 '''(INTERNAL) Get C{k2} and C{kp2} for C{._xE}, C{._yE} and C{.area}.
789 '''
790 # k2 = a2b2 / a2c2 * c2_b2
791 # kp2 = b2c2 / a2c2 * a2_b2
792 # b2 = b**2
793 # xE = Elliptic(k2, -a2b2 / b2, kp2, a2_b2)
794 # yE = Elliptic(kp2, +b2c2 / b2, k2, c2_b2)
795 # aE = Elliptic(kp2, 0, k2, 1)
796 return (self._a2b2_a2c2 * self._c2_b2,
797 self._b2c2 / self._a2c2 * self._a2_b2)
799 def _radialTo3(self, sbeta, cbeta, somega, comega):
800 '''(INTERNAL) Convert I{ellipsoidal} lat- and longitude C{beta} and
801 C{omega} to cartesian coordinates I{on the triaxial's surface},
802 also I{ordered} helper for C{.height4}.
803 '''
804 sa, ca = self._norm2(sbeta, cbeta)
805 sb, cb = self._norm2(somega, comega)
807 b2_a2 = self._1e2ab # == (b/a)**2
808 c2_a2 = -self._1e2ac # == -(c/a)**2
809 a2c2_a2 = self. e2ac # (a**2 - c**2) / a**2 == 1 - (c/a)**2
811 x = Fsum(_1_0, -b2_a2 * sa**2, c2_a2 * ca**2).fover(a2c2_a2)
812 z = Fsum(c2_a2, sb**2, b2_a2 * cb**2).fover(a2c2_a2)
814 x = self.a * cb * self._sqrt(x)
815 y = self.b * ca * sb
816 z = self.c * sa * self._sqrt(z)
817 return x, y, z
819 def reverseBetaOmega(self, x_xyz, y=None, z=None, name=NN):
820 '''Convert cartesian to I{ellipsoidal} lat- and longitude, C{beta}, C{omega}
821 and height.
823 @arg x_xyz: X component (C{scalar}) or a cartesian (C{Cartesian},
824 L{Ecef9Tuple}, L{Vector3d}, L{Vector3Tuple} or L{Vector4Tuple}).
825 @kwarg y: Y component (C{scalar}), required if B{C{x_xyz}} if C{scalar}.
826 @kwarg z: Z component (C{scalar}), required if B{C{x_xyz}} if C{scalar}.
827 @kwarg name: Optional name (C{str}).
829 @return: A L{BetaOmega3Tuple}C{(beta, omega, height)} with C{beta} and
830 C{omega} in C{Radians} and (radial) C{height} in C{meter}, same
831 units as this triaxial's axes.
833 @see: Methods L{Triaxial.forwardBetaOmega} and L{Triaxial.forwardBetaOmega_}
834 and U{Expressions (21-22)<https://www.Topo.Auth.GR/wp-content/uploads/
835 sites/111/2021/12/09_Panou.pdf>}.
836 '''
837 v = _otherV3d_(x_xyz, y, z)
838 a, b, h = self._reverseLatLon3(v, atan2, v, self.forwardBetaOmega_)
839 return BetaOmega3Tuple(Radians(beta=a), Radians(omega=b), h, name=name)
841 def reverseCartesian(self, x_xyz, y=None, z=None, h=0, normal=True, eps=_EPS2e4, name=NN):
842 '''"Unproject" a cartesian on to a cartesion I{off} this triaxial's surface.
844 @arg x_xyz: X component (C{scalar}) or a cartesian (C{Cartesian},
845 L{Ecef9Tuple}, L{Vector3d}, L{Vector3Tuple} or L{Vector4Tuple}).
846 @kwarg y: Y component (C{scalar}), required if B{C{x_xyz}} if C{scalar}.
847 @kwarg z: Z component (C{scalar}), required if B{C{x_xyz}} if C{scalar}.
848 @arg h: Height above or below this triaxial's surface (C{meter}, same units
849 as the axes).
850 @kwarg normal: If C{True} the height is C{normal} to the surface, otherwise
851 C{radially} to the center of this triaxial (C{bool}).
852 @kwarg eps: Tolerance for surface test (C{scalar}).
853 @kwarg name: Optional name (C{str}).
855 @return: A L{Vector3Tuple}C{(x, y, z)}.
857 @raise TrialError: Cartesian C{(x, y, z)} not on this triaxial's surface.
859 @see: Methods L{Triaxial.forwardCartesian} and L{Triaxial.height4}.
860 '''
861 v = _otherV3d_(x_xyz, y, z, name=name)
862 s = _sideOf(v.xyz, self._abc3, eps=eps)
863 if s: # PYCHOK no cover
864 t = _SPACE_((_inside_ if s < 0 else _outside_), self.toRepr())
865 raise TriaxialError(eps=eps, sideOf=s, x=v.x, y=v.y, z=v.z, txt=t)
867 if h:
868 if normal:
869 v = v.plus(self.normal3d(*v.xyz, length=h))
870 elif v.length > EPS0:
871 v = v.times(_1_0 + (h / v.length))
872 return v.xyz # Vector3Tuple
874 def reverseLatLon(self, x_xyz, y=None, z=None, name=NN):
875 '''Convert cartesian to I{geodetic} lat-, longitude and height.
877 @arg x_xyz: X component (C{scalar}) or a cartesian (C{Cartesian},
878 L{Ecef9Tuple}, L{Vector3d}, L{Vector3Tuple} or L{Vector4Tuple}).
879 @kwarg y: Y component (C{scalar}), required if B{C{x_xyz}} if C{scalar}.
880 @kwarg z: Z component (C{scalar}), required if B{C{x_xyz}} if C{scalar}.
881 @kwarg name: Optional name (C{str}).
883 @return: A L{LatLon3Tuple}C{(lat, lon, height)} with C{lat} and C{lon}
884 in C{degrees} and (radial) C{height} in C{meter}, same units
885 as this triaxial's axes.
887 @see: Methods L{Triaxial.forwardLatLon} and L{Triaxial.forwardLatLon_}
888 and U{Expressions (4-5)<https://www.Topo.Auth.GR/wp-content/uploads/
889 sites/111/2021/12/09_Panou.pdf>}.
890 '''
891 v = _otherV3d_(x_xyz, y, z)
892 s = v.times_(self._1e2ac, # == 1 - e_sub_x**2
893 self._1e2bc, # == 1 - e_sub_y**2
894 _1_0)
895 t = self._reverseLatLon3(s, atan2d, v, self.forwardLatLon_)
896 return LatLon3Tuple(*t, name=name)
898 def _reverseLatLon3(self, s, atan2_, v, forward_):
899 '''(INTERNAL) Helper for C{.reverseBetOmg} and C{.reverseLatLon}.
900 '''
901 x, y, z = s.xyz
902 d = hypot( x, y)
903 a = atan2_(z, d)
904 b = atan2_(y, x)
905 h = v.minus_(*forward_(z, d, y, x)).length
906 return a, b, h
909class JacobiConformal(Triaxial):
910 '''This is a conformal projection of a triaxial ellipsoid to a plane in which the
911 C{X} and C{Y} grid lines are straight.
913 Ellipsoidal coordinates I{beta} and I{omega} are converted to Jacobi Conformal
914 I{y} respectively I{x} separately. Jacobi's coordinates have been multiplied
915 by C{sqrt(B{a}**2 - B{c}**2) / (2 * B{b})} so that the customary results are
916 returned in the case of an ellipsoid of revolution (or a sphere, I{currently
917 not supported}).
919 Copyright (C) U{Charles Karney<mailto:Charles@Karney.com>} (2014-2020) and
920 licensed under the MIT/X11 License.
922 @note: This constructor can not be used to specify a sphere.
924 @see: L{Triaxial}, C++ class U{JacobiConformal<https://GeographicLib.SourceForge.io/
925 C++/doc/classGeographicLib_1_1JacobiConformal.html#details>}, U{Jacobi's conformal
926 projection<https://GeographicLib.SourceForge.io/C++/doc/jacobi.html>} and Jacobi,
927 C. G. J. I{U{Vorlesungen über Dynamik<https://Books.Google.com/books?
928 id=ryEOAAAAQAAJ&pg=PA212>}}, page 212ff,
929 '''
930# @Property_RO
931# def ab(self):
932# '''Get relative magnitude C{ab} (C{None} or C{meter}, same units as B{C{a}}).
933# '''
934# return self._ab
936# @Property_RO
937# def bc(self):
938# '''Get relative magnitude C{bc} (C{None} or C{meter}, same units as B{C{a}}).
939# '''
940# return self._bc
942 @Property_RO
943 def _xE(self):
944 '''(INTERNAL) Get the x-elliptic function.
945 '''
946 k2, kp2 = self._k2_kp2
947 # -a2b2 / b2 == (b2 - a2) / b2 == 1 - a2 / b2 == 1 - a2_b2
948 return self._Elliptic(k2, _1_0 - self._a2_b2, kp2, self._a2_b2)
950 def xR(self, omega):
951 '''Compute a Jacobi Conformal C{x} projection.
953 @arg omega: Ellipsoidal longitude (C{radians} or L{Degrees}).
955 @return: The C{x} projection (C{Radians}).
956 '''
957 return self.xR_(*SinCos2(omega))
959 def xR_(self, somega, comega):
960 '''Compute a Jacobi Conformal C{x} projection.
962 @arg somega: Ellipsoidal longitude C{sin(omega)} (C{scalar}).
963 @arg comega: Ellipsoidal longitude C{cos(omega)} (C{scalar}).
965 @return: The C{x} projection (C{Radians}).
966 '''
967 s, c = self._norm2(somega, comega, self.a)
968 return Radians(x=self._xE.fPi(s, c) * self._a2_b2)
970 @Property_RO
971 def xyQ2(self):
972 '''Get the Jacobi Conformal quadrant size (L{Jacobi2Tuple}C{(x, y)}).
973 '''
974 return Jacobi2Tuple(Radians(x=self._a2_b2 * self._xE.cPi),
975 Radians(y=self._c2_b2 * self._yE.cPi),
976 name=JacobiConformal.xyQ2.name)
978 def xyR2(self, beta, omega, name=NN):
979 '''Compute a Jacobi Conformal C{x} and C{y} projection.
981 @arg beta: Ellipsoidal latitude (C{radians} or L{Degrees}).
982 @arg omega: Ellipsoidal longitude (C{radians} or L{Degrees}).
983 @kwarg name: Optional name (C{str}).
985 @return: A L{Jacobi2Tuple}C{(x, y)}.
986 '''
987 return self.xyR2_(*(SinCos2(beta) + SinCos2(omega)),
988 name=name or self.xyR2.__name__)
990 def xyR2_(self, sbeta, cbeta, somega, comega, name=NN):
991 '''Compute a Jacobi Conformal C{x} and C{y} projection.
993 @arg sbeta: Ellipsoidal latitude C{sin(beta)} (C{scalar}).
994 @arg cbeta: Ellipsoidal latitude C{cos(beta)} (C{scalar}).
995 @arg somega: Ellipsoidal longitude C{sin(omega)} (C{scalar}).
996 @arg comega: Ellipsoidal longitude C{cos(omega)} (C{scalar}).
997 @kwarg name: Optional name (C{str}).
999 @return: A L{Jacobi2Tuple}C{(x, y)}.
1000 '''
1001 return Jacobi2Tuple(self.xR_(somega, comega),
1002 self.yR_(sbeta, cbeta),
1003 name=name or self.xyR2_.__name__)
1005 @Property_RO
1006 def _yE(self):
1007 '''(INTERNAL) Get the x-elliptic function.
1008 '''
1009 kp2, k2 = self._k2_kp2 # swapped!
1010 # b2c2 / b2 == (b2 - c2) / b2 == 1 - c2 / b2 == e2bc
1011 return self._Elliptic(k2, self.e2bc, kp2, self._c2_b2)
1013 def yR(self, beta):
1014 '''Compute a Jacobi Conformal C{y} projection.
1016 @arg beta: Ellipsoidal latitude (C{radians} or L{Degrees}).
1018 @return: The C{y} projection (C{Radians}).
1019 '''
1020 return self.yR_(*SinCos2(beta))
1022 def yR_(self, sbeta, cbeta):
1023 '''Compute a Jacobi Conformal C{y} projection.
1025 @arg sbeta: Ellipsoidal latitude C{sin(beta)} (C{scalar}).
1026 @arg cbeta: Ellipsoidal latitude C{cos(beta)} (C{scalar}).
1028 @return: The C{y} projection (C{Radians}).
1029 '''
1030 s, c = self._norm2(sbeta, cbeta, self.c)
1031 return Radians(y=self._yE.fPi(s, c) * self._c2_b2)
1034class TriaxialError(_ValueError):
1035 '''Raised for L{Triaxial} issues.
1036 '''
1037 pass # ...
1040class Triaxials(_NamedEnum):
1041 '''(INTERNAL) L{Triaxial} registry, I{must} be a sub-class
1042 to accommodate the L{_LazyNamedEnumItem} properties.
1043 '''
1044 def _Lazy(self, *abc, **name):
1045 '''(INTERNAL) Instantiate the C{Triaxial}.
1046 '''
1047 a, b, c = map(km2m, abc)
1048 return Triaxial(a, b, c, **name)
1050Triaxials = Triaxials(Triaxial, Triaxial_) # PYCHOK singleton
1051'''Some pre-defined L{Triaxial}s, all I{lazily} instantiated.'''
1052# <https://ArxIV.org/pdf/1909.06452.pdf> Table 1 Semi-axes in km
1053# <https://www.JPS.NASA.gov/education/images/pdf/ss-moons.pdf>
1054# <https://link.Springer.com/article/10.1007/s00190-022-01650-9>
1055_E = _WGS84.ellipsoid
1056Triaxials._assert( # a (km) b (km) c (km) planet
1057 Amalthea = _lazy('Amalthea', 125.0, 73.0, 64), # Jupiter
1058 Ariel = _lazy('Ariel', 581.1, 577.9, 577.7), # Uranus
1059 Earth = _lazy('Earth', 6378.173435, 6378.1039, 6356.7544),
1060 Enceladus = _lazy('Enceladus', 256.6, 251.4, 248.3), # Saturn
1061 Europa = _lazy('Europa', 1564.13, 1561.23, 1560.93), # Jupiter
1062 Io = _lazy('Io', 1829.4, 1819.3, 1815.7), # Jupiter
1063 Mars = _lazy('Mars', 3394.6, 3393.3, 3376.3),
1064 Mimas = _lazy('Mimas', 207.4, 196.8, 190.6), # Saturn
1065 Miranda = _lazy('Miranda', 240.4, 234.2, 232.9), # Uranus
1066 Moon = _lazy('Moon', 1735.55, 1735.324, 1734.898), # Earth
1067 Tethys = _lazy('Tethys', 535.6, 528.2, 525.8), # Saturn
1068 WGS84_35 = _lazy('WGS84_35', *map1(m2km, _E.a + 35, _E.a - 35, _E.b)))
1069del _E
1072def _getitems(items, *indices):
1073 '''(INTERNAL) Get the C{items} at the given I{indices}.
1075 @return: C{Type(items[i] for i in indices)} with
1076 C{Type = type(items)}, any C{type} having
1077 the special method C{__getitem__}.
1078 '''
1079 return type(items)(map(items.__getitem__, indices))
1082def _hartzell3d2(pov, los, Tun): # MCCABE 13 in .ellipsoidal.hartzell4, .formy.hartzell
1083 '''(INTERNAL) Hartzell's "Satellite Line-of-Sight Intersection ...",
1084 formula for I{un-/ordered} triaxials.
1085 '''
1086 a, b, c, T = Tun._ordered4
1088 a2 = a**2 # largest, factored out
1089 b2, p2 = (b**2, T._1e2ab) if b != a else (a2, _1_0)
1090 c2, q2 = (c**2, T._1e2ac) if c != a else (a2, _1_0)
1092 p3 = T._order3d(_otherV3d(pov=pov))
1093 u3 = T._order3d(_otherV3d(los=los)) if los else p3.negate()
1094 u3 = u3.unit() # unit vector, opposing signs
1096 x2, y2, z2 = p3.x2y2z2 # p3.times_(p3).xyz
1097 ux, vy, wz = u3.times_(p3).xyz
1098 u2, v2, w2 = u3.x2y2z2 # u3.times_(u3).xyz
1100 t = (p2 * c2), c2, b2
1101 m = fdot(t, u2, v2, w2) # a2 factored out
1102 if m < EPS0: # zero or near-null LOS vector
1103 raise _ValueError(_near_(_null_))
1105 r = fsum_(b2 * w2, c2 * v2, -v2 * z2, vy * wz * 2,
1106 -w2 * y2, b2 * u2 * q2, -u2 * z2 * p2, ux * wz * 2 * p2,
1107 -w2 * x2 * p2, -u2 * y2 * q2, -v2 * x2 * q2, ux * vy * 2 * q2, floats=True)
1108 if r > 0: # a2 factored out
1109 r = sqrt(r) * b * c # == a * a * b * c / a2
1110 elif r < 0: # LOS pointing away from or missing the triaxial
1111 raise _ValueError(_opposite_ if max(ux, vy, wz) > 0 else _outside_)
1113 d = Fdot(t, ux, vy, wz).fadd_(r).fover(m) # -r for antipode, a2 factored out
1114 if d > 0: # POV inside or LOS missing, outside the triaxial
1115 s = fsum_(_1_0, x2 / a2, y2 / b2, z2 / c2, _N_2_0, floats=True) # like _sideOf
1116 raise _ValueError(_outside_ if s > 0 else _inside_)
1117 elif fsum_(x2, y2, z2, floats=True) < d**2: # d past triaxial's center
1118 raise _ValueError(_too_(_distant_))
1120 v = p3.minus(u3.times(d)) # Vector3d
1121 h = p3.minus(v).length # distance to triaxial
1122 return T._order3d(v, reverse=True), h
1125def hartzell4(pov, los=None, tri_biax=_WGS84, name=NN):
1126 '''Compute the intersection of a tri-/biaxial ellipsoid and a Line-Of-Sight
1127 from a Point-Of-View outside.
1129 @arg pov: Point-Of-View outside the tri-/biaxial (C{Cartesian}, L{Ecef9Tuple}
1130 or L{Vector3d}).
1131 @kwarg los: Line-Of-Sight, I{direction} to the tri-/biaxial (L{Vector3d}) or
1132 C{None} to point to the tri-/biaxial's center.
1133 @kwarg tri_biax: A triaxial (L{Triaxial}, L{Triaxial_}, L{JacobiConformal})
1134 or biaxial ellipsoid (L{Datum}, L{Ellipsoid}, L{Ellipsoid2},
1135 L{a_f2Tuple} or C{scalar} radius, conventionally in C{meter}).
1136 @kwarg name: Optional name (C{str}).
1138 @return: L{Vector4Tuple}C{(x, y, z, h)} on the tri-/biaxial's surface, with C{h}
1139 the distance from B{C{pov}} to C{(x, y, z)} along B{C{los}}, all in
1140 C{meter}, conventionally.
1142 @raise TriaxialError: Null B{C{pov}} or B{C{los}}, or B{C{pov}} is inside the
1143 tri-/biaxial or B{C{los}} points outside the tri-/biaxial
1144 or points in an opposite direction.
1146 @raise TypeError: Invalid B{C{pov}} or B{C{los}}.
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 T = D.ellipsoid._triaxial
1159 try:
1160 v, h = _hartzell3d2(pov, los, T)
1161 except Exception as x:
1162 raise TriaxialError(pov=pov, los=los, tri_biax=tri_biax, cause=x)
1163 return Vector4Tuple(v.x, v.y, v.z, h, name=name or hartzell4.__name__)
1166def _normalTo4(x, y, a, b, eps=EPS): # MCCABE 14
1167 '''(INTERNAL) Nearest point on and distance to a 2-D ellipse, I{unordered}.
1169 @see: Function C{pygeodesy.ellipsoids._normalTo3} and I{Eberly}'s U{Distance
1170 from a Point to ... an Ellipsoid ...<https://www.GeometricTools.com/
1171 Documentation/DistancePointEllipseEllipsoid.pdf>}.
1172 '''
1173 def _root2d(r, u, v, g, eps):
1174 # robust root finder
1175 _1, __2 = _1_0, _0_5
1176 _a, _h2 = fabs, _hypot21_
1177 u *= r
1178 t0 = v - _1
1179 t1 = _0_0 if g < 0 else _h2(u, v)
1180 for i in range(1, _TRIPS):
1181 t = (t0 + t1) * __2
1182 if t in (t0, t1) or _a(t0 - t1) < eps:
1183 break
1184 g = _h2(u / (t + r), v / (t + _1))
1185 if g > 0:
1186 t0 = t
1187 elif g < 0:
1188 t1 = t
1189 else:
1190 break
1191 else: # PYCHOK no cover
1192 e = _a(t0 - t1)
1193 t = _root2d.__name__
1194 raise _ValueError(Fmt.no_convergence(e, eps), txt=t)
1195 return t, i
1197 if a < b:
1198 b, a, d, i = _normalTo4(y, x, b, a, eps=eps)
1199 return a, b, d, i
1201 if not (isfinite(a) and b > 0):
1202 raise _ValueError(a=a, b=b)
1204 i = None
1205 if y:
1206 if x:
1207 u = fabs(x / a)
1208 v = fabs(y / b)
1209 g = _hypot21_(u, v)
1210 if g:
1211 r = (a / b)**2
1212 t, i = _root2d(r, u, v, g, eps)
1213 a = x / (t / r + _1_0)
1214 b = y / (t + _1_0)
1215 d = hypot(x - a, y - b)
1216 else: # on the ellipse
1217 a, b, d = x, y, _0_0
1218 else: # x == 0
1219 if y < 0:
1220 b = -b
1221 a, d = x, fabs(y - b)
1223 else: # y == 0
1224 n = a * x
1225 d = (a + b) * (a - b)
1226 if d > fabs(n): # PYCHOK no cover
1227 r = n / d
1228 a *= r
1229 b *= sqrt(_1_0 - r**2)
1230 d = hypot(x - a, b)
1231 else:
1232 if x < 0:
1233 a = -a
1234 b, d = y, fabs(x - a)
1235 return a, b, d, i
1238def _normalTo5(x, y, z, Tun, eps=EPS): # MCCABE 24
1239 '''(INTERNAL) Nearest point on and distance to an I{un- or ordered} triaxial.
1241 @see: I{Eberly}'s U{Distance from a Point to ... an Ellipsoid ...<https://
1242 www.GeometricTools.com/Documentation/DistancePointEllipseEllipsoid.pdf>}.
1243 '''
1244 def _root3d(r, s, u, v, w, g, eps):
1245 # robust root finder
1246 _1, __2 = _1_0, _0_5
1247 _a, _h2 = fabs, _hypot21_
1248 u *= r
1249 v *= s
1250 t0 = w - _1
1251 t1 = _0_0 if g < 0 else _h2(u, v, w)
1252 for i in range(1, _TRIPS):
1253 t = (t0 + t1) * __2
1254 if t in (t0, t1) or _a(t0 - t1) < eps:
1255 break
1256 g = _h2(u / (t + r), v / (t + s), w / (t + _1))
1257 if g > 0:
1258 t0 = t
1259 elif g < 0:
1260 t1 = t
1261 else:
1262 break
1263 else: # PYCHOK no cover
1264 e = _a(t0 - t1)
1265 t = _root3d.__name__
1266 raise _ValueError(Fmt.no_convergence(e, eps), txt=t)
1267 return t, i
1269 a, b, c, T = Tun._ordered4
1270 if Tun is not T: # T is ordered, Tun isn't
1271 t = T._order3(x, y, z) + (T,)
1272 a, b, c, d, i = _normalTo5(*t, eps=eps)
1273 return T._order3(a, b, c, reverse=True) + (d, i)
1275 if not (isfinite(a) and c > 0):
1276 raise _ValueError(a=a, b=b, c=c)
1278 if eps > 0:
1279 val = max(eps * 1e8, EPS)
1280 else: # no validation
1281 val, eps = 0, -eps
1283 i = None
1284 if z:
1285 if y:
1286 if x:
1287 u = fabs(x / a)
1288 v = fabs(y / b)
1289 w = fabs(z / c)
1290 g = _hypot21_(u, v, w)
1291 if g:
1292 r = T._1e2ac # (c / a)**2
1293 s = T._1e2bc # (c / b)**2
1294 t, i = _root3d(_1_0 / r, _1_0 / s, u, v, w, g, eps)
1295 a = x / (t * r + _1_0)
1296 b = y / (t * s + _1_0)
1297 c = z / (t + _1_0)
1298 d = hypot_(x - a, y - b, z - c)
1299 else: # on the ellipsoid
1300 a, b, c, d = x, y, z, _0_0
1301 else: # x == 0
1302 a = x # 0
1303 b, c, d, i = _normalTo4(y, z, b, c, eps=eps)
1304 elif x: # y == 0
1305 b = y # 0
1306 a, c, d, i = _normalTo4(x, z, a, c, eps=eps)
1307 else: # x == y == 0
1308 if z < 0:
1309 c = -c
1310 a, b, d = x, y, fabs(z - c)
1312 else: # z == 0
1313 t = False
1314 n = a * x
1315 d = T._a2c2 # (a + c) * (a - c)
1316 if d > fabs(n):
1317 u = n / d
1318 n = b * y
1319 d = T._b2c2 # (b + c) * (b - c)
1320 if d > fabs(n):
1321 v = n / d
1322 n = _hypot21_(u, v)
1323 if n < 0:
1324 a *= u
1325 b *= v
1326 c *= sqrt(-n)
1327 d = hypot_(x - a, y - b, c)
1328 t = True
1329 if not t:
1330 c = z # 0
1331 a, b, d, i = _normalTo4(x, y, a, b, eps=eps)
1333 if val > 0: # validate
1334 e = T.sideOf(a, b, c, eps=val)
1335 if e: # not near the ellipsoid's surface
1336 raise _ValueError(a=a, b=b, c=c, d=d,
1337 sideOf=e, eps=val)
1338 if d: # angle of delta and normal vector
1339 m = Vector3d(x, y, z).minus_(a, b, c)
1340 if m.euclid > val:
1341 m = m.unit()
1342 n = T.normal3d(a, b, c)
1343 e = n.dot(m) # n.negate().dot(m)
1344 if not isnear1(fabs(e), eps1=val):
1345 raise _ValueError(n=n, m=m,
1346 dot=e, eps=val)
1347 return a, b, c, d, i
1350def _otherV3d_(x_xyz, y, z, name=NN):
1351 '''(INTERNAL) Get a Vector3d from C{x_xyz}, C{y} and C{z}.
1352 '''
1353 return Vector3d(x_xyz, y, z, name=name) if isscalar(x_xyz) else \
1354 _otherV3d(x_xyz=x_xyz)
1357def _sideOf(xyz, abc, eps=EPS): # in .formy
1358 '''(INTERNAL) Helper for C{_hartzell3d2}, M{.sideOf} and M{.reverseCartesian}.
1360 @return: M{sum((x / a)**2 for x, a in zip(xyz, abc)) - 1} or C{INT0},
1361 '''
1362 s = _hypot21_(*((x / a) for x, a in _zip(xyz, abc) if a)) # strict=True
1363 return s if fabs(s) > eps else INT0
1366if __name__ == '__main__':
1368 from pygeodesy import printf
1370 # __doc__ of this file, force all into registery
1371 t = [NN] + Triaxials.toRepr(all=True, asorted=True).split(_NL_)
1372 printf(_NLATvar_.join(i.strip(_COMMA_) for i in t))
1374# **) MIT License
1375#
1376# Copyright (C) 2022-2023 -- mrJean1 at Gmail -- All Rights Reserved.
1377#
1378# Permission is hereby granted, free of charge, to any person obtaining a
1379# copy of this software and associated documentation files (the "Software"),
1380# to deal in the Software without restriction, including without limitation
1381# the rights to use, copy, modify, merge, publish, distribute, sublicense,
1382# and/or sell copies of the Software, and to permit persons to whom the
1383# Software is furnished to do so, subject to the following conditions:
1384#
1385# The above copyright notice and this permission notice shall be included
1386# in all copies or substantial portions of the Software.
1387#
1388# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
1389# OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
1390# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
1391# THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR
1392# OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE,
1393# ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR
1394# OTHER DEALINGS IN THE SOFTWARE.