Coverage for pygeodesy/triaxials.py: 96%
584 statements
« prev ^ index » next coverage.py v7.2.2, created at 2024-01-10 13:43 -0500
« prev ^ index » next coverage.py v7.2.2, created at 2024-01-10 13:43 -0500
2# -*- coding: utf-8 -*-
4u'''Triaxal ellipsoid classes I{ordered} L{Triaxial} and I{unordered} L{Triaxial_} and Jacobi
5conformal projections L{JacobiConformal} and L{JacobiConformalSpherical}, transcoded from
6I{Charles Karney}'s C++ class U{JacobiConformal<https://GeographicLib.SourceForge.io/C++/doc/
7classGeographicLib_1_1JacobiConformal.html#details>} to pure Python and miscellaneous classes
8L{BetaOmega2Tuple}, L{BetaOmega3Tuple}, L{Jacobi2Tuple} and L{TriaxialError}.
10Copyright (C) U{Charles Karney<mailto:Karney@Alum.MIT.edu>} (2008-2023). For more information,
11see the U{GeographicLib<https://GeographicLib.SourceForge.io>} documentation.
13@see: U{Geodesics on a triaxial ellipsoid<https://WikiPedia.org/wiki/Geodesics_on_an_ellipsoid#
14 Geodesics_on_a_triaxial_ellipsoid>} and U{Triaxial coordinate systems and their geometrical
15 interpretation<https://www.Topo.Auth.GR/wp-content/uploads/sites/111/2021/12/09_Panou.pdf>}.
17@var Triaxials.Amalthea: Triaxial(name='Amalthea', a=125000, b=73000, c=64000, e2ab=0.658944, e2bc=0.231375493, e2ac=0.737856, volume=2446253479595252, area=93239507787.490371704, area_p=93212299402.670425415)
18@var Triaxials.Ariel: Triaxial(name='Ariel', a=581100, b=577900, c=577700, e2ab=0.01098327, e2bc=0.000692042, e2ac=0.011667711, volume=812633172614203904, area=4211301462766.580078125, area_p=4211301574065.829589844)
19@var Triaxials.Earth: Triaxial(name='Earth', a=6378173.435, b=6378103.9, c=6356754.399999999, e2ab=0.000021804, e2bc=0.006683418, e2ac=0.006705077, volume=1083208241574987694080, area=510065911057441.0625, area_p=510065915922713.6875)
20@var Triaxials.Enceladus: Triaxial(name='Enceladus', a=256600, b=251400, c=248300, e2ab=0.040119337, e2bc=0.024509841, e2ac=0.06364586, volume=67094551514082248, area=798618496278.596679688, area_p=798619018175.109863281)
21@var Triaxials.Europa: Triaxial(name='Europa', a=1564130, b=1561230, c=1560930, e2ab=0.003704694, e2bc=0.000384275, e2ac=0.004087546, volume=15966575194402123776, area=30663773697323.51953125, area_p=30663773794562.45703125)
22@var Triaxials.Io: Triaxial(name='Io', a=1829400, b=1819300, c=1815700, e2ab=0.011011391, e2bc=0.003953651, e2ac=0.014921506, volume=25313121117889765376, area=41691875849096.7421875, area_p=41691877397441.2109375)
23@var Triaxials.Mars: Triaxial(name='Mars', a=3394600, b=3393300, c=3376300, e2ab=0.000765776, e2bc=0.009994646, e2ac=0.010752768, volume=162907283585817247744, area=144249140795107.4375, area_p=144249144150662.15625)
24@var Triaxials.Mimas: Triaxial(name='Mimas', a=207400, b=196800, c=190600, e2ab=0.09960581, e2bc=0.062015624, e2ac=0.155444317, volume=32587072869017956, area=493855762247.691894531, area_p=493857714107.9375)
25@var Triaxials.Miranda: Triaxial(name='Miranda', a=240400, b=234200, c=232900, e2ab=0.050915557, e2bc=0.011070811, e2ac=0.061422691, volume=54926187094835456, area=698880863325.756958008, area_p=698881306767.950317383)
26@var Triaxials.Moon: Triaxial(name='Moon', a=1735550, b=1735324, c=1734898, e2ab=0.000260419, e2bc=0.000490914, e2ac=0.000751206, volume=21886698675223740416, area=37838824729886.09375, area_p=37838824733332.2265625)
27@var Triaxials.Tethys: Triaxial(name='Tethys', a=535600, b=528200, c=525800, e2ab=0.027441672, e2bc=0.009066821, e2ac=0.036259685, volume=623086233855821440, area=3528073490771.394042969, area_p=3528074261832.738769531)
28@var Triaxials.WGS84_35: Triaxial(name='WGS84_35', a=6378172, b=6378102, c=6356752.314245179, e2ab=0.00002195, e2bc=0.006683478, e2ac=0.006705281, volume=1083207319768789942272, area=510065621722018.125, area_p=510065626587483.3125)
29'''
30# make sure int/int division yields float quotient, see .basics
31from __future__ import division as _; del _ # PYCHOK semicolon
33from pygeodesy.basics import isLatLon, isscalar, map2, _zip, _ValueError
34from pygeodesy.constants import EPS, EPS0, EPS02, EPS4, _EPS2e4, INT0, PI2, PI_3, PI4, \
35 _0_0, _0_5, _1_0, _N_2_0, float0_, isfinite, isnear1, \
36 _4_0 # PYCHOK used!
37from pygeodesy.datums import Datum, _spherical_datum, _WGS84, Ellipsoid, _EWGS84, Fmt
38# from pygeodesy.dms import toDMS # _MODS
39# from pygeodesy.ellipsoids import Ellipsoid, _EWGS84 # from .datums
40# from pygeodesy.elliptic import Elliptic # _MODS
41# from pygeodesy.errors import _ValueError # from .basics
42from pygeodesy.fmath import Fdot, fdot, fmean_, hypot, hypot_, norm2, sqrt0
43from pygeodesy.fsums import Fsum, fsumf_, fsum1f_
44from pygeodesy.interns import NN, _a_, _b_, _beta_, _c_, _distant_, _finite_, \
45 _height_, _inside_, _near_, _not_, _NOTEQUAL_, _null_, \
46 _opposite_, _outside_, _SPACE_, _spherical_, _too_, \
47 _x_, _y_
48# from pygeodesy.lazily import _ALL_LAZY, _ALL_MODS as _MODS # from .vector3d
49from pygeodesy.named import _NamedEnum, _NamedEnumItem, _NamedTuple, _Pass, \
50 _lazyNamedEnumItem as _lazy
51from pygeodesy.namedTuples import LatLon3Tuple, Vector3Tuple, Vector4Tuple
52from pygeodesy.props import Property_RO, property_RO
53# from pygeodesy.streprs import Fmt # from .datums
54from pygeodesy.units import Float, Height_, Meter, Meter2, Meter3, Radians, \
55 Radius, Scalar_, _toDegrees, _toRadians
56from pygeodesy.utily import asin1, atan2d, km2m, m2km, SinCos2, sincos2d_
57from pygeodesy.vector3d import _otherV3d, Vector3d, _ALL_LAZY, _MODS
59from math import atan2, fabs, sqrt
61__all__ = _ALL_LAZY.triaxials
62__version__ = '24.01.06'
64_not_ordered_ = _not_('ordered')
65_omega_ = 'omega'
66_TRIPS = 537 # 52..58, Eberly 1074?
69class _NamedTupleTo(_NamedTuple): # in .testNamedTuples
70 '''(INTERNAL) Base for C{-.toDegrees}, C{-.toRadians}.
71 '''
72 def _toDegrees(self, a, b, *c, **toDMS_kwds):
73 a, b, _ = _toDegrees(self, a, b, **toDMS_kwds)
74 return _ or self.classof(a, b, *c, name=self.name)
76 def _toRadians(self, a, b, *c):
77 a, b, _ = _toRadians(self, a, b)
78 return _ or self.classof(a, b, *c, name=self.name)
81class BetaOmega2Tuple(_NamedTupleTo):
82 '''2-Tuple C{(beta, omega)} with I{ellipsoidal} lat- and
83 longitude C{beta} and C{omega} both in L{Radians} (or
84 L{Degrees}).
85 '''
86 _Names_ = (_beta_, _omega_)
87 _Units_ = (_Pass, _Pass)
89 def toDegrees(self, **toDMS_kwds):
90 '''Convert this L{BetaOmega2Tuple} to L{Degrees} or C{toDMS}.
92 @return: L{BetaOmega2Tuple}C{(beta, omega)} with
93 C{beta} and C{omega} both in L{Degrees}
94 or as a L{toDMS} string provided some
95 B{C{toDMS_kwds}} keyword arguments are
96 specified.
97 '''
98 return _NamedTupleTo._toDegrees(self, *self, **toDMS_kwds)
100 def toRadians(self):
101 '''Convert this L{BetaOmega2Tuple} to L{Radians}.
103 @return: L{BetaOmega2Tuple}C{(beta, omega)} with
104 C{beta} and C{omega} both in L{Radians}.
105 '''
106 return _NamedTupleTo._toRadians(self, *self)
109class BetaOmega3Tuple(_NamedTupleTo):
110 '''3-Tuple C{(beta, omega, height)} with I{ellipsoidal} lat- and
111 longitude C{beta} and C{omega} both in L{Radians} (or L{Degrees})
112 and the C{height}, rather the (signed) I{distance} to the triaxial's
113 surface (measured along the radial line to the triaxial's center)
114 in C{meter}, conventionally.
115 '''
116 _Names_ = BetaOmega2Tuple._Names_ + (_height_,)
117 _Units_ = BetaOmega2Tuple._Units_ + ( Meter,)
119 def toDegrees(self, **toDMS_kwds):
120 '''Convert this L{BetaOmega3Tuple} to L{Degrees} or C{toDMS}.
122 @return: L{BetaOmega3Tuple}C{(beta, omega, height)} with
123 C{beta} and C{omega} both in L{Degrees} or as a
124 L{toDMS} string provided some B{C{toDMS_kwds}}
125 keyword arguments are specified.
126 '''
127 return _NamedTupleTo._toDegrees(self, *self, **toDMS_kwds)
129 def toRadians(self):
130 '''Convert this L{BetaOmega3Tuple} to L{Radians}.
132 @return: L{BetaOmega3Tuple}C{(beta, omega, height)} with
133 C{beta} and C{omega} both in L{Radians}.
134 '''
135 return _NamedTupleTo._toRadians(self, *self)
137 def to2Tuple(self):
138 '''Reduce this L{BetaOmega3Tuple} to a L{BetaOmega2Tuple}.
139 '''
140 return BetaOmega2Tuple(*self[:2])
143class Jacobi2Tuple(_NamedTupleTo):
144 '''2-Tuple C{(x, y)} with a Jacobi Conformal C{x} and C{y}
145 projection, both in L{Radians} (or L{Degrees}).
146 '''
147 _Names_ = (_x_, _y_)
148 _Units_ = (_Pass, _Pass)
150 def toDegrees(self, **toDMS_kwds):
151 '''Convert this L{Jacobi2Tuple} to L{Degrees} or C{toDMS}.
153 @return: L{Jacobi2Tuple}C{(x, y)} with C{x} and C{y}
154 both in L{Degrees} or as a L{toDMS} string
155 provided some B{C{toDMS_kwds}} keyword
156 arguments are specified.
157 '''
158 return _NamedTupleTo._toDegrees(self, *self, **toDMS_kwds)
160 def toRadians(self):
161 '''Convert this L{Jacobi2Tuple} to L{Radians}.
163 @return: L{Jacobi2Tuple}C{(x, y)} with C{x}
164 and C{y} both in L{Radians}.
165 '''
166 return _NamedTupleTo._toRadians(self, *self)
169class Triaxial_(_NamedEnumItem):
170 '''I{Unordered} triaxial ellipsoid and base class.
172 Triaxial ellipsoids with right-handed semi-axes C{a}, C{b} and C{c}, oriented
173 such that the large principal ellipse C{ab} is the equator I{Z}=0, I{beta}=0,
174 while the small principal ellipse C{ac} is the prime meridian, plane I{Y}=0,
175 I{omega}=0.
177 The four umbilic points, C{abs}(I{omega}) = C{abs}(I{beta}) = C{PI/2}, lie on
178 the middle principal ellipse C{bc} in plane I{X}=0, I{omega}=C{PI/2}.
180 @note: I{Geodetic} C{lat}- and C{lon}gitudes are in C{degrees}, I{geodetic}
181 C{phi} and C{lam}bda are in C{radians}, but I{ellipsoidal} lat- and
182 longitude C{beta} and C{omega} are in L{Radians} by default (or in
183 L{Degrees} if converted).
184 '''
185 _ijk = _kji = None
186 _unordered = True
188 def __init__(self, a_triaxial, b=None, c=None, name=NN):
189 '''New I{unordered} L{Triaxial_}.
191 @arg a_triaxial: Large, C{X} semi-axis (C{scalar}, conventionally in
192 C{meter}) or an other L{Triaxial} or L{Triaxial_} instance.
193 @kwarg b: Middle, C{Y} semi-axis (C{meter}, same units as B{C{a}}), required
194 if C{B{a_triaxial} is scalar}, ignored otherwise.
195 @kwarg c: Small, C{Z} semi-axis (C{meter}, same units as B{C{a}}), required
196 if C{B{a_triaxial} is scalar}, ignored otherwise.
197 @kwarg name: Optional name (C{str}).
199 @raise TriaxialError: Invalid semi-axis or -axes.
200 '''
201 try:
202 a = a_triaxial
203 t = a._abc3 if isinstance(a, Triaxial_) else (
204 Radius(a=a), Radius(b=b), Radius(c=c))
205 except (TypeError, ValueError) as x:
206 raise TriaxialError(a=a, b=b, c=c, cause=x)
207 if name:
208 self.name = name
210 a, b, c = self._abc3 = t
211 if self._unordered: # == not isinstance(self, Triaxial)
212 s, _, t = sorted(t)
213 if not (isfinite(t) and s > 0):
214 raise TriaxialError(a=a, b=b, c=c) # txt=_invalid_
215 elif not (isfinite(a) and a >= b >= c > 0):
216 raise TriaxialError(a=a, b=b, c=c, txt=_not_ordered_)
217 elif not (a > c and self._a2c2 > 0 and self.e2ac > 0):
218 raise TriaxialError(a=a, c=c, e2ac=self.e2ac, txt=_spherical_)
220 def __str__(self):
221 return self.toStr()
223 @Property_RO
224 def a(self):
225 '''Get the (largest) C{x} semi-axis (C{meter}, conventionally).
226 '''
227 a, _, _ = self._abc3
228 return a
230 @Property_RO
231 def _a2b2(self):
232 '''(INTERNAL) Get C{a**2 - b**2} == E_sub_e**2.
233 '''
234 a, b, _ = self._abc3
235 return ((a - b) * (a + b)) if a != b else _0_0
237 @Property_RO
238 def _a2_b2(self):
239 '''(INTERNAL) Get C{(a/b)**2}.
240 '''
241 a, b, _ = self._abc3
242 return (a / b)**2 if a != b else _1_0
244 @Property_RO
245 def _a2c2(self):
246 '''(INTERNAL) Get C{a**2 - c**2} == E_sub_x**2.
247 '''
248 a, _, c = self._abc3
249 return ((a - c) * (a + c)) if a != c else _0_0
251 @Property_RO
252 def area(self):
253 '''Get the surface area (C{meter} I{squared}).
254 '''
255 c, b, a = sorted(self._abc3)
256 if a > c:
257 a = Triaxial(a, b, c).area if a > b else \
258 Ellipsoid(a, b=c).areax # a == b
259 else: # a == c == b
260 a = Meter2(area=a**2 * PI4)
261 return a
263 def area_p(self, p=1.6075):
264 '''I{Approximate} the surface area (C{meter} I{squared}).
266 @kwarg p: Exponent (C{scalar} > 0), 1.6 for near-spherical or 1.5849625007
267 for "near-flat" triaxials.
269 @see: U{Surface area<https://WikiPedia.org/wiki/Ellipsoid#Approximate_formula>}.
270 '''
271 a, b, c = self._abc3
272 if a == b == c:
273 a *= a
274 else:
275 _p = pow
276 a = _p(fmean_(_p(a * b, p), _p(a * c, p), _p(b * c, p)), _1_0 / p)
277 return Meter2(area_p=a * PI4)
279 @Property_RO
280 def b(self):
281 '''Get the (middle) C{y} semi-axis (C{meter}, same units as B{C{a}}).
282 '''
283 _, b, _ = self._abc3
284 return b
286 @Property_RO
287 def _b2c2(self):
288 '''(INTERNAL) Get C{b**2 - c**2} == E_sub_y**2.
289 '''
290 _, b, c = self._abc3
291 return ((b - c) * (b + c)) if b != c else _0_0
293 @Property_RO
294 def c(self):
295 '''Get the (smallest) C{z} semi-axis (C{meter}, same units as B{C{a}}).
296 '''
297 _, _, c = self._abc3
298 return c
300 @Property_RO
301 def _c2_b2(self):
302 '''(INTERNAL) Get C{(c/b)**2}.
303 '''
304 _, b, c = self._abc3
305 return (c / b)**2 if b != c else _1_0
307 @Property_RO
308 def e2ab(self):
309 '''Get the C{ab} ellipse' I{(1st) eccentricity squared} (C{scalar}), M{1 - (b/a)**2}.
310 '''
311 return Float(e2ab=(_1_0 - self._1e2ab) or _0_0)
313 @Property_RO
314 def _1e2ab(self):
315 '''(INTERNAL) Get C{1 - e2ab} == C{(b/a)**2}.
316 '''
317 a, b, _ = self._abc3
318 return (b / a)**2 if a != b else _1_0
320 @Property_RO
321 def e2ac(self):
322 '''Get the C{ac} ellipse' I{(1st) eccentricity squared} (C{scalar}), M{1 - (c/a)**2}.
323 '''
324 return Float(e2ac=(_1_0 - self._1e2ac) or _0_0)
326 @Property_RO
327 def _1e2ac(self):
328 '''(INTERNAL) Get C{1 - e2ac} == C{(c/a)**2}.
329 '''
330 a, _, c = self._abc3
331 return (c / a)**2 if a != c else _1_0
333 @Property_RO
334 def e2bc(self):
335 '''Get the C{bc} ellipse' I{(1st) eccentricity squared} (C{scalar}), M{1 - (c/b)**2}.
336 '''
337 return Float(e2bc=(_1_0 - self._1e2bc) or _0_0)
339 _1e2bc = _c2_b2 # C{1 - e2bc} == C{(c/b)**2}
341 @property_RO
342 def _Elliptic(self):
343 '''(INTERNAL) Get class L{Elliptic}, I{once}.
344 '''
345 Triaxial_._Elliptic = E = _MODS.elliptic.Elliptic # overwrite property_RO
346 return E
348 def hartzell4(self, pov, los=None, name=NN):
349 '''Compute the intersection of this triaxial's surface with a Line-Of-Sight
350 from a Point-Of-View in space.
352 @see: Function L{pygeodesy.hartzell4} for further details.
353 '''
354 return hartzell4(pov, los=los, tri_biax=self, name=name)
356 def height4(self, x_xyz, y=None, z=None, normal=True, eps=EPS):
357 '''Compute the projection on and the height of a cartesian above or below
358 this triaxial's surface.
360 @arg x_xyz: X component (C{scalar}) or a cartesian (C{Cartesian},
361 L{Ecef9Tuple}, L{Vector3d}, L{Vector3Tuple} or L{Vector4Tuple}).
362 @kwarg y: Y component (C{scalar}), required if B{C{x_xyz}} if C{scalar}.
363 @kwarg z: Z component (C{scalar}), required if B{C{x_xyz}} if C{scalar}.
364 @kwarg normal: If C{True} the projection is perpendicular to (the nearest
365 point on) this triaxial's surface, otherwise the C{radial}
366 line to this triaxial's center (C{bool}).
367 @kwarg eps: Tolerance for root finding and validation (C{scalar}), use a
368 negative value to skip validation.
370 @return: L{Vector4Tuple}C{(x, y, z, h)} with the cartesian coordinates
371 C{x}, C{y} and C{z} of the projection on or the intersection
372 with and with the height C{h} above or below the triaxial's
373 surface in C{meter}, conventionally.
375 @raise TriaxialError: Non-cartesian B{C{xyz}}, invalid B{C{eps}}, no
376 convergence in root finding or validation failed.
378 @see: Method L{Ellipsoid.height4} and I{Eberly}'s U{Distance from a Point
379 to ... an Ellipsoid ...<https://www.GeometricTools.com/Documentation/
380 DistancePointEllipseEllipsoid.pdf>}.
381 '''
382 v, r = _otherV3d_(x_xyz, y, z), self.isSpherical
384 i, h = None, v.length
385 if h < EPS0: # EPS
386 x = y = z = _0_0
387 h -= min(self._abc3) # nearest
388 elif r: # .isSpherical
389 x, y, z = v.times(r / h).xyz
390 h -= r
391 else:
392 x, y, z = v.xyz
393 try:
394 if normal: # perpendicular to triaxial
395 x, y, z, h, i = _normalTo5(x, y, z, self, eps=eps)
396 else: # radially to triaxial's center
397 x, y, z = self._radialTo3(z, hypot(x, y), y, x)
398 h = v.minus_(x, y, z).length
399 except Exception as e:
400 raise TriaxialError(x=x, y=y, z=z, cause=e)
401 if h > 0 and self.sideOf(v, eps=EPS0) < 0:
402 h = -h # below the surface
403 return Vector4Tuple(x, y, z, h, iteration=i, name=self.height4.__name__)
405 @Property_RO
406 def isOrdered(self):
407 '''Is this triaxial I{ordered} and I{not spherical} (C{bool})?
408 '''
409 a, b, c = self._abc3
410 return bool(a >= b > c) # b > c!
412 @Property_RO
413 def isSpherical(self):
414 '''Is this triaxial I{spherical} (C{Radius} or INT0)?
415 '''
416 a, b, c = self._abc3
417 return a if a == b == c else INT0
419 def _norm2(self, s, c, *a):
420 '''(INTERNAL) Normalize C{s} and C{c} iff not already.
421 '''
422 if fabs(_hypot21(s, c)) > EPS02:
423 s, c = norm2(s, c)
424 if a:
425 s, c = norm2(s * self.b, c * a[0])
426 return float0_(s, c)
428 def normal3d(self, x_xyz, y=None, z=None, length=_1_0):
429 '''Get a 3-D vector perpendicular to at a cartesian on this triaxial's surface.
431 @arg x_xyz: X component (C{scalar}) or a cartesian (C{Cartesian},
432 L{Ecef9Tuple}, L{Vector3d}, L{Vector3Tuple} or L{Vector4Tuple}).
433 @kwarg y: Y component (C{scalar}), required if B{C{x_xyz}} if C{scalar}.
434 @kwarg z: Z component (C{scalar}), required if B{C{x_xyz}} if C{scalar}.
435 @kwarg length: Optional length and in-/outward direction (C{scalar}).
437 @return: A C{Vector3d(x_, y_, z_)} normalized to B{C{length}}, pointing
438 in- or outward for neg- respectively positive B{C{length}}.
440 @note: Cartesian location C{(B{x}, B{y}, B{z})} must be on this triaxial's
441 surface, use method L{Triaxial.sideOf} to validate.
442 '''
443 # n = 2 * (x / a2, y / b2, z / c2)
444 # == 2 * (x, y * a2 / b2, z * a2 / c2) / a2 # iff ordered
445 # == 2 * (x, y / _1e2ab, z / _1e2ac) / a2
446 # == unit(x, y / _1e2ab, z / _1e2ac).times(length)
447 n = self._normal3d.times_(*_otherV3d_(x_xyz, y, z).xyz)
448 if n.length < EPS0:
449 raise TriaxialError(x=x_xyz, y=y, z=z, txt=_null_)
450 return n.times(length / n.length)
452 @Property_RO
453 def _normal3d(self):
454 '''(INTERNAL) Get M{Vector3d((d/a)**2, (d/b)**2, (d/c)**2)}, M{d = max(a, b, c)}.
455 '''
456 d = max(self._abc3)
457 t = tuple(((d / x)**2 if x != d else _1_0) for x in self._abc3)
458 return Vector3d(*t, name=self.normal3d.__name__)
460 def _order3(self, *abc, **reverse): # reverse=False
461 '''(INTERNAL) Un-/Order C{a}, C{b} and C{c}.
463 @return: 3-Tuple C{(a, b, c)} ordered by or un-ordered
464 (reverse-ordered) C{ijk} if C{B{reverse}=True}.
465 '''
466 ijk = self._order_ijk(**reverse)
467 return _getitems(abc, *ijk) if ijk else abc
469 def _order3d(self, v, **reverse): # reverse=False
470 '''(INTERNAL) Un-/Order a C{Vector3d}.
472 @return: Vector3d(x, y, z) un-/ordered.
473 '''
474 ijk = self._order_ijk(**reverse)
475 return v.classof(*_getitems(v.xyz, *ijk)) if ijk else v
477 @Property_RO
478 def _ordered4(self):
479 '''(INTERNAL) Helper for C{_hartzell2} and C{_normalTo5}.
480 '''
481 def _order2(reverse, a, b, c):
482 '''(INTERNAL) Un-Order C{a}, C{b} and C{c}.
484 @return: 2-Tuple C{((a, b, c), ijk)} with C{a} >= C{b} >= C{c}
485 and C{ijk} a 3-tuple with the initial indices.
486 '''
487 i, j, k = 0, 1, 2 # range(3)
488 if a < b:
489 a, b, i, j = b, a, j, i
490 if a < c:
491 a, c, i, k = c, a, k, i
492 if b < c:
493 b, c, j, k = c, b, k, j
494 # reverse (k, j, i) since (a, b, c) is reversed-sorted
495 ijk = (k, j, i) if reverse else (None if i < j < k else (i, j, k))
496 return (a, b, c), ijk
498 abc, T = self._abc3, self
499 if not self.isOrdered:
500 abc, ijk = _order2(False, *abc)
501 if ijk:
502 _, kji = _order2(True, *ijk)
503 T = Triaxial_(*abc)
504 T._ijk, T._kji = ijk, kji
505 return abc + (T,)
507 def _order_ijk(self, reverse=False):
508 '''(INTERNAL) Get the un-/order indices.
509 '''
510 return self._kji if reverse else self._ijk
512 def _radialTo3(self, sbeta, cbeta, somega, comega):
513 '''(INTERNAL) I{Unordered} helper for C{.height4}.
514 '''
515 def _rphi(a, b, sphi, cphi):
516 # <https://WikiPedia.org/wiki/Ellipse#Polar_form_relative_to_focus>
517 # polar form: radius(phi) = a * b / hypot(a * sphi, b * cphi)
518 return (b / hypot(sphi, b / a * cphi)) if a > b else (
519 (a / hypot(cphi, a / b * sphi)) if a < b else a)
521 sa, ca = self._norm2(sbeta, cbeta)
522 sb, cb = self._norm2(somega, comega)
524 a, b, c = self._abc3
525 if a != b:
526 a = _rphi(a, b, sb, cb)
527 if a != c:
528 c = _rphi(a, c, sa, ca)
529 z, r = c * sa, c * ca
530 x, y = r * cb, r * sb
531 return x, y, z
533 def sideOf(self, x_xyz, y=None, z=None, eps=EPS4):
534 '''Is a cartesian above, below or on the surface of this triaxial?
536 @arg x_xyz: X component (C{scalar}) or a cartesian (C{Cartesian},
537 L{Ecef9Tuple}, L{Vector3d}, L{Vector3Tuple} or L{Vector4Tuple}).
538 @kwarg y: Y component (C{scalar}), required if B{C{x_xyz}} if C{scalar}.
539 @kwarg z: Z component (C{scalar}), required if B{C{x_xyz}} if C{scalar}.
540 @kwarg eps: Near surface tolerance(C{scalar}).
542 @return: C{INT0} if C{(B{x}, B{y}, B{z})} is near this triaxial's surface
543 within tolerance B{C{eps}}, otherwise a neg- or positive C{float}
544 if in- respectively outside this triaxial.
546 @see: Methods L{Triaxial.height4} and L{Triaxial.normal3d}.
547 '''
548 return _sideOf(_otherV3d_(x_xyz, y, z).xyz, self._abc3, eps=eps)
550 def toEllipsoid(self, name=NN):
551 '''Convert this triaxial to an L{Ellipsoid}, provided 2 axes match.
553 @return: An L{Ellipsoid} with north along this C{Z} axis if C{a == b},
554 this C{Y} axis if C{a == c} or this C{X} axis if C{b == c}.
556 @raise TriaxialError: This C{a != b}, C{b != c} and C{c != a}.
558 @see: Method L{Ellipsoid.toTriaxial}.
559 '''
560 a, b, c = self._abc3
561 if a == b:
562 b = c # N = c-Z
563 elif b == c: # N = a-X
564 a, b = b, a
565 elif a != c: # N = b-Y
566 t = _SPACE_(_a_, _NOTEQUAL_, _b_, _NOTEQUAL_, _c_)
567 raise TriaxialError(a=a, b=b, c=c, txt=t)
568 return Ellipsoid(a, b=b, name=name or self.name)
570 def toStr(self, prec=9, name=NN, **unused): # PYCHOK signature
571 '''Return this C{Triaxial} as a string.
573 @kwarg prec: Precision, number of decimal digits (0..9).
574 @kwarg name: Override name (C{str}) or C{None} to exclude
575 this triaxial's name.
577 @return: This C{Triaxial}'s attributes (C{str}).
578 '''
579 T = Triaxial_
580 t = T.a,
581 J = JacobiConformalSpherical
582 t += (J.ab, J.bc) if isinstance(self, J) else (T.b, T.c)
583 t += T.e2ab, T.e2bc, T.e2ac
584 J = JacobiConformal
585 if isinstance(self, J):
586 t += J.xyQ2,
587 t += T.volume, T.area
588 return self._instr(name, prec, props=t, area_p=self.area_p())
590 @Property_RO
591 def volume(self):
592 '''Get the volume (C{meter**3}), M{4 / 3 * PI * a * b * c}.
593 '''
594 a, b, c = self._abc3
595 return Meter3(volume=a * b * c * PI_3 * _4_0)
598class Triaxial(Triaxial_):
599 '''I{Ordered} triaxial ellipsoid.
601 @see: L{Triaxial_} for more information.
602 '''
603 _unordered = False
605 def __init__(self, a_triaxial, b=None, c=None, name=NN):
606 '''New I{ordered} L{Triaxial}.
608 @arg a_triaxial: Largest semi-axis (C{scalar}, conventionally in C{meter})
609 or an other L{Triaxial} or L{Triaxial_} instance.
610 @kwarg b: Middle semi-axis (C{meter}, same units as B{C{a}}), required
611 if C{B{a_triaxial} is scalar}, ignored otherwise.
612 @kwarg c: Smallest semi-axis (C{meter}, same units as B{C{a}}), required
613 if C{B{a_triaxial} is scalar}, ignored otherwise.
614 @kwarg name: Optional name (C{str}).
616 @note: The semi-axes must be ordered as C{B{a} >= B{b} >= B{c} > 0} and
617 must be ellipsoidal, C{B{a} > B{c}}.
619 @raise TriaxialError: Semi-axes not ordered, spherical or invalid.
620 '''
621 Triaxial_.__init__(self, a_triaxial, b=b, c=c, name=name)
623 @Property_RO
624 def _a2b2_a2c2(self):
625 '''@see: Methods C{.forwardBetaOmega} and C{._k2_kp2}.
626 '''
627 return self._a2b2 / self._a2c2
629 @Property_RO
630 def area(self):
631 '''Get the surface area (C{meter} I{squared}).
633 @see: U{Surface area<https://WikiPedia.org/wiki/Ellipsoid#Surface_area>}.
634 '''
635 a, b, c = self._abc3
636 if a != b:
637 kp2, k2 = self._k2_kp2 # swapped!
638 aE = self._Elliptic(k2, _0_0, kp2, _1_0)
639 c2 = self._1e2ac # cos(phi)**2 = (c/a)**2
640 s = sqrt(self.e2ac) # sin(phi)**2 = 1 - c2
641 r = asin1(s) # phi = atan2(sqrt(c2), s)
642 b *= fsum1f_(aE.fE(r) * s, c / a * c / b,
643 aE.fF(r) * c2 / s)
644 a = Meter2(area=a * b * PI2)
645 else: # a == b > c
646 a = Ellipsoid(a, b=c).areax
647 return a
649 def _exyz3(self, u):
650 '''(INTERNAL) Helper for C{.forwardBetOmg}.
651 '''
652 if u > 0:
653 u2 = u**2
654 x = u * sqrt0(_1_0 + self._a2c2 / u2, Error=TriaxialError)
655 y = u * sqrt0(_1_0 + self._b2c2 / u2, Error=TriaxialError)
656 else:
657 x = y = u = _0_0
658 return x, y, u
660 def forwardBetaOmega(self, beta, omega, height=0, name=NN):
661 '''Convert I{ellipsoidal} lat- and longitude C{beta}, C{omega}
662 and height to cartesian.
664 @arg beta: Ellipsoidal latitude (C{radians} or L{Degrees}).
665 @arg omega: Ellipsoidal longitude (C{radians} or L{Degrees}).
666 @arg height: Height above or below the ellipsoid's surface (C{meter}, same
667 units as this triaxial's C{a}, C{b} and C{c} semi-axes).
668 @kwarg name: Optional name (C{str}).
670 @return: A L{Vector3Tuple}C{(x, y, z)}.
672 @see: Method L{Triaxial.reverseBetaOmega} and U{Expressions (23-25)<https://
673 www.Topo.Auth.GR/wp-content/uploads/sites/111/2021/12/09_Panou.pdf>}.
674 '''
675 if height:
676 h = self._Height(height)
677 x, y, z = self._exyz3(h + self.c)
678 else:
679 x, y, z = self._abc3 # == self._exyz3(self.c)
680 if z: # and x and y:
681 sa, ca = SinCos2(beta)
682 sb, cb = SinCos2(omega)
684 r = self._a2b2_a2c2
685 x *= cb * sqrt0(ca**2 + r * sa**2, Error=TriaxialError)
686 y *= ca * sb
687 z *= sa * sqrt0(_1_0 - r * cb**2, Error=TriaxialError)
688 return Vector3Tuple(x, y, z, name=name)
690 def forwardBetaOmega_(self, sbeta, cbeta, somega, comega, name=NN):
691 '''Convert I{ellipsoidal} lat- and longitude C{beta} and C{omega}
692 to cartesian coordinates I{on the triaxial's surface}.
694 @arg sbeta: Ellipsoidal latitude C{sin(beta)} (C{scalar}).
695 @arg cbeta: Ellipsoidal latitude C{cos(beta)} (C{scalar}).
696 @arg somega: Ellipsoidal longitude C{sin(omega)} (C{scalar}).
697 @arg comega: Ellipsoidal longitude C{cos(omega)} (C{scalar}).
698 @kwarg name: Optional name (C{str}).
700 @return: A L{Vector3Tuple}C{(x, y, z)} on the surface.
702 @raise TriaxialError: This triaxial is near-spherical.
704 @see: Method L{Triaxial.reverseBetaOmega}, U{Triaxial ellipsoid coordinate
705 system<https://WikiPedia.org/wiki/Geodesics_on_an_ellipsoid#
706 Triaxial_ellipsoid_coordinate_system>} and U{expressions (23-25)<https://
707 www.Topo.Auth.GR/wp-content/uploads/sites/111/2021/12/09_Panou.pdf>}.
708 '''
709 t = self._radialTo3(sbeta, cbeta, somega, comega)
710 return Vector3Tuple(*t, name=name)
712 def forwardCartesian(self, x_xyz, y=None, z=None, name=NN, **normal_eps):
713 '''Project a cartesian on this triaxial.
715 @arg x_xyz: X component (C{scalar}) or a cartesian (C{Cartesian},
716 L{Ecef9Tuple}, L{Vector3d}, L{Vector3Tuple} or L{Vector4Tuple}).
717 @kwarg y: Y component (C{scalar}), required if B{C{x_xyz}} if C{scalar}.
718 @kwarg z: Z component (C{scalar}), required if B{C{x_xyz}} if C{scalar}.
719 @kwarg name: Optional name (C{str}).
720 @kwarg normal_eps: Optional keyword arguments C{B{normal}=True} and
721 C{B{eps}=EPS}, see method L{Triaxial.height4}.
723 @see: Method L{Triaxial.height4} for further information and method
724 L{Triaxial.reverseCartesian} to reverse the projection.
725 '''
726 t = self.height4(x_xyz, y, z, **normal_eps)
727 _ = t.rename(name)
728 return t
730 def forwardLatLon(self, lat, lon, height=0, name=NN):
731 '''Convert I{geodetic} lat-, longitude and heigth to cartesian.
733 @arg lat: Geodetic latitude (C{degrees}).
734 @arg lon: Geodetic longitude (C{degrees}).
735 @arg height: Height above the ellipsoid (C{meter}, same units
736 as this triaxial's C{a}, C{b} and C{c} axes).
737 @kwarg name: Optional name (C{str}).
739 @return: A L{Vector3Tuple}C{(x, y, z)}.
741 @see: Method L{Triaxial.reverseLatLon} and U{Expressions (9-11)<https://
742 www.Topo.Auth.GR/wp-content/uploads/sites/111/2021/12/09_Panou.pdf>}.
743 '''
744 return self._forwardLatLon3(height, name, *sincos2d_(lat, lon))
746 def forwardLatLon_(self, slat, clat, slon, clon, height=0, name=NN):
747 '''Convert I{geodetic} lat-, longitude and heigth to cartesian.
749 @arg slat: Geodetic latitude C{sin(lat)} (C{scalar}).
750 @arg clat: Geodetic latitude C{cos(lat)} (C{scalar}).
751 @arg slon: Geodetic longitude C{sin(lon)} (C{scalar}).
752 @arg clon: Geodetic longitude C{cos(lon)} (C{scalar}).
753 @arg height: Height above the ellipsoid (C{meter}, same units
754 as this triaxial's axes C{a}, C{b} and C{c}).
755 @kwarg name: Optional name (C{str}).
757 @return: A L{Vector3Tuple}C{(x, y, z)}.
759 @see: Method L{Triaxial.reverseLatLon} and U{Expressions (9-11)<https://
760 www.Topo.Auth.GR/wp-content/uploads/sites/111/2021/12/09_Panou.pdf>}.
761 '''
762 sa, ca = self._norm2(slat, clat)
763 sb, cb = self._norm2(slon, clon)
764 return self._forwardLatLon3(height, name, sa, ca, sb, cb)
766 def _forwardLatLon3(self, height, name, sa, ca, sb, cb):
767 '''(INTERNAL) Helper for C{.forwardLatLon} and C{.forwardLatLon_}.
768 '''
769 ca_x_sb = ca * sb
770 h = self._Height(height)
771 # 1 - (1 - (c/a)**2) * sa**2 - (1 - (b/a)**2) * ca**2 * sb**2
772 t = fsumf_(_1_0, -self.e2ac * sa**2, -self.e2ab * ca_x_sb**2)
773 n = self.a / sqrt0(t, Error=TriaxialError) # prime vertical
774 x = (h + n) * ca * cb
775 y = (h + n * self._1e2ab) * ca_x_sb
776 z = (h + n * self._1e2ac) * sa
777 return Vector3Tuple(x, y, z, name=name)
779 def _Height(self, height):
780 '''(INTERNAL) Validate a C{height}.
781 '''
782 return Height_(height=height, low=-self.c, Error=TriaxialError)
784 @Property_RO
785 def _k2_kp2(self):
786 '''(INTERNAL) Get C{k2} and C{kp2} for C{._xE}, C{._yE} and C{.area}.
787 '''
788 # k2 = a2b2 / a2c2 * c2_b2
789 # kp2 = b2c2 / a2c2 * a2_b2
790 # b2 = b**2
791 # xE = Elliptic(k2, -a2b2 / b2, kp2, a2_b2)
792 # yE = Elliptic(kp2, +b2c2 / b2, k2, c2_b2)
793 # aE = Elliptic(kp2, 0, k2, 1)
794 return (self._a2b2_a2c2 * self._c2_b2,
795 self._b2c2 / self._a2c2 * self._a2_b2)
797 def _radialTo3(self, sbeta, cbeta, somega, comega):
798 '''(INTERNAL) Convert I{ellipsoidal} lat- and longitude C{beta} and
799 C{omega} to cartesian coordinates I{on the triaxial's surface},
800 also I{ordered} helper for C{.height4}.
801 '''
802 sa, ca = self._norm2(sbeta, cbeta)
803 sb, cb = self._norm2(somega, comega)
805 b2_a2 = self._1e2ab # == (b/a)**2
806 c2_a2 = -self._1e2ac # == -(c/a)**2
807 a2c2_a2 = self. e2ac # (a**2 - c**2) / a**2 == 1 - (c/a)**2
809 x2 = Fsum(_1_0, -b2_a2 * sa**2, c2_a2 * ca**2).fover(a2c2_a2)
810 z2 = Fsum(c2_a2, sb**2, b2_a2 * cb**2).fover(a2c2_a2)
812 x, y, z = self._abc3
813 x *= cb * sqrt0(x2, Error=TriaxialError)
814 y *= ca * sb
815 z *= sa * sqrt0(z2, Error=TriaxialError)
816 return x, y, z
818 def reverseBetaOmega(self, x_xyz, y=None, z=None, name=NN):
819 '''Convert cartesian to I{ellipsoidal} lat- and longitude, C{beta}, C{omega}
820 and height.
822 @arg x_xyz: X component (C{scalar}) or a cartesian (C{Cartesian},
823 L{Ecef9Tuple}, L{Vector3d}, L{Vector3Tuple} or L{Vector4Tuple}).
824 @kwarg y: Y component (C{scalar}), required if B{C{x_xyz}} if C{scalar}.
825 @kwarg z: Z component (C{scalar}), required if B{C{x_xyz}} if C{scalar}.
826 @kwarg name: Optional name (C{str}).
828 @return: A L{BetaOmega3Tuple}C{(beta, omega, height)} with C{beta} and
829 C{omega} in L{Radians} and (radial) C{height} in C{meter}, same
830 units as this triaxial's axes.
832 @see: Methods L{Triaxial.forwardBetaOmega} and L{Triaxial.forwardBetaOmega_}
833 and U{Expressions (21-22)<https://www.Topo.Auth.GR/wp-content/uploads/
834 sites/111/2021/12/09_Panou.pdf>}.
835 '''
836 v = _otherV3d_(x_xyz, y, z)
837 a, b, h = self._reverseLatLon3(v, atan2, v, self.forwardBetaOmega_)
838 return BetaOmega3Tuple(Radians(beta=a), Radians(omega=b), h, name=name)
840 def reverseCartesian(self, x_xyz, y=None, z=None, h=0, normal=True, eps=_EPS2e4, name=NN):
841 '''"Unproject" a cartesian on to a cartesion I{off} this triaxial's surface.
843 @arg x_xyz: X component (C{scalar}) or a cartesian (C{Cartesian},
844 L{Ecef9Tuple}, L{Vector3d}, L{Vector3Tuple} or L{Vector4Tuple}).
845 @kwarg y: Y component (C{scalar}), required if B{C{x_xyz}} if C{scalar}.
846 @kwarg z: Z component (C{scalar}), required if B{C{x_xyz}} if C{scalar}.
847 @arg h: Height above or below this triaxial's surface (C{meter}, same units
848 as the axes).
849 @kwarg normal: If C{True} the height is C{normal} to the surface, otherwise
850 C{radially} to the center of this triaxial (C{bool}).
851 @kwarg eps: Tolerance for surface test (C{scalar}).
852 @kwarg name: Optional name (C{str}).
854 @return: A L{Vector3Tuple}C{(x, y, z)}.
856 @raise TrialError: Cartesian C{(x, y, z)} not on this triaxial's surface.
858 @see: Methods L{Triaxial.forwardCartesian} and L{Triaxial.height4}.
859 '''
860 v = _otherV3d_(x_xyz, y, z, name=name)
861 s = _sideOf(v.xyz, self._abc3, eps=eps)
862 if s: # PYCHOK no cover
863 t = _SPACE_((_inside_ if s < 0 else _outside_), self.toRepr())
864 raise TriaxialError(eps=eps, sideOf=s, x=v.x, y=v.y, z=v.z, txt=t)
866 if h:
867 if normal:
868 v = v.plus(self.normal3d(*v.xyz, length=h))
869 elif v.length > EPS0:
870 v = v.times(_1_0 + (h / v.length))
871 return v.xyz # Vector3Tuple
873 def reverseLatLon(self, x_xyz, y=None, z=None, name=NN):
874 '''Convert cartesian to I{geodetic} lat-, longitude and height.
876 @arg x_xyz: X component (C{scalar}) or a cartesian (C{Cartesian},
877 L{Ecef9Tuple}, L{Vector3d}, L{Vector3Tuple} or L{Vector4Tuple}).
878 @kwarg y: Y component (C{scalar}), required if B{C{x_xyz}} if C{scalar}.
879 @kwarg z: Z component (C{scalar}), required if B{C{x_xyz}} if C{scalar}.
880 @kwarg name: Optional name (C{str}).
882 @return: A L{LatLon3Tuple}C{(lat, lon, height)} with C{lat} and C{lon}
883 in C{degrees} and (radial) C{height} in C{meter}, same units
884 as this triaxial's axes.
886 @see: Methods L{Triaxial.forwardLatLon} and L{Triaxial.forwardLatLon_}
887 and U{Expressions (4-5)<https://www.Topo.Auth.GR/wp-content/uploads/
888 sites/111/2021/12/09_Panou.pdf>}.
889 '''
890 v = _otherV3d_(x_xyz, y, z)
891 s = v.times_(self._1e2ac, # == 1 - e_sub_x**2
892 self._1e2bc, # == 1 - e_sub_y**2
893 _1_0)
894 t = self._reverseLatLon3(s, atan2d, v, self.forwardLatLon_)
895 return LatLon3Tuple(*t, name=name)
897 def _reverseLatLon3(self, s, atan2_, v, forward_):
898 '''(INTERNAL) Helper for C{.reverseBetOmg} and C{.reverseLatLon}.
899 '''
900 x, y, z = s.xyz
901 d = hypot( x, y)
902 a = atan2_(z, d)
903 b = atan2_(y, x)
904 h = v.minus_(*forward_(z, d, y, x)).length
905 return a, b, h
908class JacobiConformal(Triaxial):
909 '''This is a conformal projection of a triaxial ellipsoid to a plane in which the
910 C{X} and C{Y} grid lines are straight.
912 Ellipsoidal coordinates I{beta} and I{omega} are converted to Jacobi Conformal
913 I{y} respectively I{x} separately. Jacobi's coordinates have been multiplied
914 by C{sqrt(B{a}**2 - B{c}**2) / (2 * B{b})} so that the customary results are
915 returned in the case of an ellipsoid of revolution.
917 Copyright (C) U{Charles Karney<mailto:Karney@Alum.MIT.edu>} (2014-2023) and
918 licensed under the MIT/X11 License.
920 @note: This constructor can I{not be used to specify a sphere}, see alternate
921 L{JacobiConformalSpherical}.
923 @see: L{Triaxial}, C++ class U{JacobiConformal<https://GeographicLib.SourceForge.io/
924 C++/doc/classGeographicLib_1_1JacobiConformal.html#details>}, U{Jacobi's conformal
925 projection<https://GeographicLib.SourceForge.io/C++/doc/jacobi.html>} and Jacobi,
926 C. G. J. I{U{Vorlesungen über Dynamik<https://Books.Google.com/books?
927 id=ryEOAAAAQAAJ&pg=PA212>}}, page 212ff.
928 '''
930 @Property_RO
931 def _xE(self):
932 '''(INTERNAL) Get the x-elliptic function.
933 '''
934 k2, kp2 = self._k2_kp2
935 # -a2b2 / b2 == (b2 - a2) / b2 == 1 - a2 / b2 == 1 - a2_b2
936 return self._Elliptic(k2, _1_0 - self._a2_b2, kp2, self._a2_b2)
938 def xR(self, omega):
939 '''Compute a Jacobi Conformal C{x} projection.
941 @arg omega: Ellipsoidal longitude (C{radians} or L{Degrees}).
943 @return: The C{x} projection (L{Radians}).
944 '''
945 return self.xR_(*SinCos2(omega))
947 def xR_(self, somega, comega):
948 '''Compute a Jacobi Conformal C{x} projection.
950 @arg somega: Ellipsoidal longitude C{sin(omega)} (C{scalar}).
951 @arg comega: Ellipsoidal longitude C{cos(omega)} (C{scalar}).
953 @return: The C{x} projection (L{Radians}).
954 '''
955 s, c = self._norm2(somega, comega, self.a)
956 return Radians(x=self._xE.fPi(s, c) * self._a2_b2)
958 @Property_RO
959 def xyQ2(self):
960 '''Get the Jacobi Conformal quadrant size (L{Jacobi2Tuple}C{(x, y)}).
961 '''
962 return Jacobi2Tuple(Radians(x=self._a2_b2 * self._xE.cPi),
963 Radians(y=self._c2_b2 * self._yE.cPi),
964 name=JacobiConformal.xyQ2.name)
966 def xyR2(self, beta, omega, name=NN):
967 '''Compute a Jacobi Conformal C{x} and C{y} projection.
969 @arg beta: Ellipsoidal latitude (C{radians} or L{Degrees}).
970 @arg omega: Ellipsoidal longitude (C{radians} or L{Degrees}).
971 @kwarg name: Optional name (C{str}).
973 @return: A L{Jacobi2Tuple}C{(x, y)}.
974 '''
975 return self.xyR2_(*(SinCos2(beta) + SinCos2(omega)),
976 name=name or self.xyR2.__name__)
978 def xyR2_(self, sbeta, cbeta, somega, comega, name=NN):
979 '''Compute a Jacobi Conformal C{x} and C{y} projection.
981 @arg sbeta: Ellipsoidal latitude C{sin(beta)} (C{scalar}).
982 @arg cbeta: Ellipsoidal latitude C{cos(beta)} (C{scalar}).
983 @arg somega: Ellipsoidal longitude C{sin(omega)} (C{scalar}).
984 @arg comega: Ellipsoidal longitude C{cos(omega)} (C{scalar}).
985 @kwarg name: Optional name (C{str}).
987 @return: A L{Jacobi2Tuple}C{(x, y)}.
988 '''
989 return Jacobi2Tuple(self.xR_(somega, comega),
990 self.yR_(sbeta, cbeta),
991 name=name or self.xyR2_.__name__)
993 @Property_RO
994 def _yE(self):
995 '''(INTERNAL) Get the x-elliptic function.
996 '''
997 kp2, k2 = self._k2_kp2 # swapped!
998 # b2c2 / b2 == (b2 - c2) / b2 == 1 - c2 / b2 == e2bc
999 return self._Elliptic(k2, self.e2bc, kp2, self._c2_b2)
1001 def yR(self, beta):
1002 '''Compute a Jacobi Conformal C{y} projection.
1004 @arg beta: Ellipsoidal latitude (C{radians} or L{Degrees}).
1006 @return: The C{y} projection (L{Radians}).
1007 '''
1008 return self.yR_(*SinCos2(beta))
1010 def yR_(self, sbeta, cbeta):
1011 '''Compute a Jacobi Conformal C{y} projection.
1013 @arg sbeta: Ellipsoidal latitude C{sin(beta)} (C{scalar}).
1014 @arg cbeta: Ellipsoidal latitude C{cos(beta)} (C{scalar}).
1016 @return: The C{y} projection (L{Radians}).
1017 '''
1018 s, c = self._norm2(sbeta, cbeta, self.c)
1019 return Radians(y=self._yE.fPi(s, c) * self._c2_b2)
1022class JacobiConformalSpherical(JacobiConformal):
1023 '''An alternate, I{spherical} L{JacobiConformal} projection.
1025 @see: L{JacobiConformal} for other and more details.
1026 '''
1027 _ab = _bc = 0
1029 def __init__(self, radius_triaxial, ab=0, bc=0, name=NN):
1030 '''New L{JacobiConformalSpherical}.
1032 @arg radius_triaxial: Radius (C{scalar}, conventionally in
1033 C{meter}) or an other L{JacobiConformalSpherical},
1034 L{JacobiConformal} or ordered L{Triaxial}.
1035 @kwarg ab: Relative magnitude of C{B{a} - B{b}} (C{meter},
1036 same units as C{scalar B{radius}}.
1037 @kwarg bc: Relative magnitude of C{B{b} - B{c}} (C{meter},
1038 same units as C{scalar B{radius}}.
1039 @kwarg name: Optional name (C{str}).
1041 @raise TriaxialError: Invalid B{C{radius_triaxial}}, negative
1042 B{C{ab}}, negative B{C{bc}} or C{(B{ab}
1043 + B{bc})} not positive.
1045 @note: If B{C{radius_triaxial}} is a L{JacobiConformalSpherical}
1046 and if B{C{ab}} and B{C{bc}} are both zero or C{None},
1047 the B{C{radius_triaxial}}'s C{ab}, C{bc}, C{a}, C{b}
1048 and C{c} are copied.
1049 '''
1050 try:
1051 r, j = radius_triaxial, False
1052 if isinstance(r, Triaxial): # ordered only
1053 if (not (ab or bc)) and isinstance(r, JacobiConformalSpherical):
1054 j = True
1055 t = r._abc3
1056 else:
1057 t = (Radius(radius=r),) * 3
1058 self._ab = r.ab if j else Scalar_(ab=ab) # low=0
1059 self._bc = r.bc if j else Scalar_(bc=bc) # low=0
1060 if (self.ab + self.bc) <= 0:
1061 raise ValueError('(ab + bc)')
1062 a, _, c = self._abc3 = t
1063 if not (a >= c and isfinite(self._a2b2)
1064 and isfinite(self._a2c2)):
1065 raise ValueError(_not_(_finite_))
1066 except (TypeError, ValueError) as x:
1067 raise TriaxialError(radius_triaxial=r, ab=ab, bc=bc, cause=x)
1068 if name:
1069 self.name = name
1071 @Property_RO
1072 def ab(self):
1073 '''Get relative magnitude C{ab} (C{meter}, same units as B{C{a}}).
1074 '''
1075 return self._ab
1077 @Property_RO
1078 def _a2b2(self):
1079 '''(INTERNAL) Get C{a**2 - b**2} == ab * (a + b).
1080 '''
1081 a, b, _ = self._abc3
1082 return self.ab * (a + b)
1084 @Property_RO
1085 def _a2c2(self):
1086 '''(INTERNAL) Get C{a**2 - c**2} == a2b2 + b2c2.
1087 '''
1088 return self._a2b2 + self._b2c2
1090 @Property_RO
1091 def bc(self):
1092 '''Get relative magnitude C{bc} (C{meter}, same units as B{C{a}}).
1093 '''
1094 return self._bc
1096 @Property_RO
1097 def _b2c2(self):
1098 '''(INTERNAL) Get C{b**2 - c**2} == bc * (b + c).
1099 '''
1100 _, b, c = self._abc3
1101 return self.bc * (b + c)
1103 @Property_RO
1104 def radius(self):
1105 '''Get radius (C{meter}, conventionally).
1106 '''
1107 return self.a
1110class TriaxialError(_ValueError):
1111 '''Raised for L{Triaxial} issues.
1112 '''
1113 pass # ...
1116class Triaxials(_NamedEnum):
1117 '''(INTERNAL) L{Triaxial} registry, I{must} be a sub-class
1118 to accommodate the L{_LazyNamedEnumItem} properties.
1119 '''
1120 def _Lazy(self, *abc, **name):
1121 '''(INTERNAL) Instantiate the C{Triaxial}.
1122 '''
1123 a, b, c = map(km2m, abc)
1124 return Triaxial(a, b, c, **name)
1126Triaxials = Triaxials(Triaxial, Triaxial_) # PYCHOK singleton
1127'''Some pre-defined L{Triaxial}s, all I{lazily} instantiated.'''
1128# <https://ArxIV.org/pdf/1909.06452.pdf> Table 1 Semi-axes in Km
1129# <https://www.JPS.NASA.gov/education/images/pdf/ss-moons.pdf>
1130# <https://link.Springer.com/article/10.1007/s00190-022-01650-9>
1131_EWGS84_35 = _EWGS84.a + 35, _EWGS84.a - 35, _EWGS84.b
1132Triaxials._assert( # a (Km) b (Km) c (Km) planet
1133 Amalthea = _lazy('Amalthea', 125.0, 73.0, 64), # Jupiter
1134 Ariel = _lazy('Ariel', 581.1, 577.9, 577.7), # Uranus
1135 Earth = _lazy('Earth', 6378.173435, 6378.1039, 6356.7544),
1136 Enceladus = _lazy('Enceladus', 256.6, 251.4, 248.3), # Saturn
1137 Europa = _lazy('Europa', 1564.13, 1561.23, 1560.93), # Jupiter
1138 Io = _lazy('Io', 1829.4, 1819.3, 1815.7), # Jupiter
1139 Mars = _lazy('Mars', 3394.6, 3393.3, 3376.3),
1140 Mimas = _lazy('Mimas', 207.4, 196.8, 190.6), # Saturn
1141 Miranda = _lazy('Miranda', 240.4, 234.2, 232.9), # Uranus
1142 Moon = _lazy('Moon', 1735.55, 1735.324, 1734.898), # Earth
1143 Tethys = _lazy('Tethys', 535.6, 528.2, 525.8), # Saturn
1144 WGS84_35 = _lazy('WGS84_35', *map2(m2km, _EWGS84_35)))
1145del _EWGS84_35
1148def _getitems(items, *indices):
1149 '''(INTERNAL) Get the C{items} at the given I{indices}.
1151 @return: C{Type(items[i] for i in indices)} with
1152 C{Type = type(items)}, any C{type} having
1153 the special method C{__getitem__}.
1154 '''
1155 return type(items)(map(items.__getitem__, indices))
1158def _hartzell2(pov, los, Tun): # in .ellipsoids.hartzell4, .formy.hartzell
1159 '''(INTERNAL) Hartzell's "Satellite Line-of-Sight Intersection ...",
1160 formula from a Point-Of-View to an I{un-/ordered} Triaxial.
1161 '''
1162 def _toUvwV3d(los, pov):
1163 try: # pov must be LatLon or Cartesian if los is a Los
1164 v = los.toUvw(pov)
1165 except (AttributeError, TypeError):
1166 v = _otherV3d(los=los)
1167 return v
1169 p3 = _otherV3d(pov=pov.toCartesian() if isLatLon(pov) else pov)
1170 u3 = _toUvwV3d(los, pov) if los else p3.negate()
1172 a, b, c, T = Tun._ordered4
1174 a2 = a**2 # largest, factored out
1175 b2, p2 = (b**2, T._1e2ab) if b != a else (a2, _1_0)
1176 c2, q2 = (c**2, T._1e2ac) if c != a else (a2, _1_0)
1178 p3 = T._order3d(p3)
1179 u3 = T._order3d(u3).unit() # unit vector, opposing signs
1181 x2, y2, z2 = p3.x2y2z2 # p3.times_(p3).xyz
1182 ux, vy, wz = u3.times_(p3).xyz
1183 u2, v2, w2 = u3.x2y2z2 # u3.times_(u3).xyz
1185 t = (p2 * c2), c2, b2
1186 m = fdot(t, u2, v2, w2) # a2 factored out
1187 if m < EPS0: # zero or near-null LOS vector
1188 raise _ValueError(_near_(_null_))
1190 r = fsumf_(b2 * w2, c2 * v2, -v2 * z2, vy * wz * 2,
1191 -w2 * y2, -u2 * y2 * q2, -u2 * z2 * p2, ux * wz * 2 * p2,
1192 -w2 * x2 * p2, b2 * u2 * q2, -v2 * x2 * q2, ux * vy * 2 * q2)
1193 if r > 0: # a2 factored out
1194 r = sqrt(r) * b * c # == a * a * b * c / a2
1195 elif r < 0: # LOS pointing away from or missing the triaxial
1196 raise _ValueError(_opposite_ if max(ux, vy, wz) > 0 else _outside_)
1198 d = Fdot(t, ux, vy, wz).fadd_(r).fover(m) # -r for antipode, a2 factored out
1199 if d > 0: # POV inside or LOS outside or missing the triaxial
1200 s = fsumf_(_1_0, x2 / a2, y2 / b2, z2 / c2, _N_2_0) # like _sideOf
1201 raise _ValueError(_outside_ if s > 0 else _inside_)
1202 elif fsum1f_(x2, y2, z2) < d**2: # d past triaxial's center
1203 raise _ValueError(_too_(_distant_))
1205 v = p3.minus(u3.times(d)) # cartesian type(pov) or Vector3d
1206 h = p3.minus(v).length # distance to pov == -d
1207 return T._order3d(v, reverse=True), h
1210def hartzell4(pov, los=None, tri_biax=_WGS84, name=NN):
1211 '''Compute the intersection of a tri-/biaxial ellipsoid and a Line-Of-Sight
1212 from a Point-Of-View outside.
1214 @arg pov: Point-Of-View outside the tri-/biaxial (C{Cartesian}, L{Ecef9Tuple}
1215 C{LatLon} or L{Vector3d}).
1216 @kwarg los: Line-Of-Sight, I{direction} to the tri-/biaxial (L{Los}, L{Vector3d})
1217 or C{None} to point to the tri-/biaxial's center.
1218 @kwarg tri_biax: A triaxial (L{Triaxial}, L{Triaxial_}, L{JacobiConformal} or
1219 L{JacobiConformalSpherical}) or biaxial ellipsoid (L{Datum},
1220 L{Ellipsoid}, L{Ellipsoid2}, L{a_f2Tuple} or C{scalar} radius,
1221 conventionally in C{meter}).
1222 @kwarg name: Optional name (C{str}).
1224 @return: L{Vector4Tuple}C{(x, y, z, h)} on the tri-/biaxial's surface, with C{h}
1225 the distance from B{C{pov}} to C{(x, y, z)} I{along the} B{C{los}}, all
1226 in C{meter}, conventionally.
1228 @raise TriaxialError: Null B{C{pov}} or B{C{los}}, or B{C{pov}} is inside the
1229 tri-/biaxial or B{C{los}} points outside the tri-/biaxial
1230 or points in an opposite direction.
1232 @raise TypeError: Invalid B{C{pov}} or B{C{los}}.
1234 @see: Function L{pygeodesy.hartzell}, L{pygeodesy.tyr3d}, class L{pygeodesy.Los}
1235 and U{lookAtSpheroid<https://PyPI.org/project/pymap3d>} and U{I{Satellite
1236 Line-of-Sight Intersection with Earth}<https://StephenHartzell.Medium.com/
1237 satellite-line-of-sight-intersection-with-earth-d786b4a6a9b6>}.
1238 '''
1239 if isinstance(tri_biax, Triaxial_):
1240 T = tri_biax
1241 else:
1242 D = tri_biax if isinstance(tri_biax, Datum) else \
1243 _spherical_datum(tri_biax, name=hartzell4.__name__)
1244 T = D.ellipsoid._triaxial
1246 try:
1247 v, h = _hartzell2(pov, los, T)
1248 except Exception as x:
1249 raise TriaxialError(pov=pov, los=los, tri_biax=tri_biax, cause=x)
1250 return Vector4Tuple(v.x, v.y, v.z, h, name=name or hartzell4.__name__)
1253def _hypot21(x, y, z=0):
1254 '''(INTERNAL) Compute M{x**2 + y**2 + z**2 - 1} with C{max(fabs(x),
1255 fabs(y), fabs(z))} rarely greater than 1.0.
1256 '''
1257 return fsumf_(_1_0, x**2, y**2, (z**2 if z else _0_0), _N_2_0)
1260def _normalTo4(x, y, a, b, eps=EPS):
1261 '''(INTERNAL) Nearest point on and distance to a 2-D ellipse, I{unordered}.
1263 @see: Function C{pygeodesy.ellipsoids._normalTo3} and I{Eberly}'s U{Distance
1264 from a Point to ... an Ellipse ...<https://www.GeometricTools.com/
1265 Documentation/DistancePointEllipseEllipsoid.pdf>}.
1266 '''
1267 if b > a:
1268 b, a, d, i = _normalTo4(y, x, b, a, eps=eps)
1269 return a, b, d, i
1271 if not (b > 0 and isfinite(a)):
1272 raise _ValueError(a=a, b=b)
1274 i, _a = None, fabs
1275 if y:
1276 if x:
1277 u = _a(x / a)
1278 v = _a(y / b)
1279 g = _hypot21(u, v)
1280 if g:
1281 r = (a / b)**2
1282 t, i = _rootXd(r, 0, u, 0, v, g, eps)
1283 a = x / (t / r + _1_0)
1284 b = y / (t + _1_0)
1285 d = hypot(x - a, y - b)
1286 else: # on the ellipse
1287 a, b, d = x, y, _0_0
1288 else: # x == 0
1289 if y < 0:
1290 b = -b
1291 a, d = x, _a(y - b)
1293 else: # y == 0
1294 n = a * x
1295 d = (a + b) * (a - b)
1296 if d > _a(n): # PYCHOK no cover
1297 r = n / d
1298 a *= r
1299 b *= sqrt(_1_0 - r**2)
1300 d = hypot(x - a, b)
1301 else:
1302 if x < 0:
1303 a = -a
1304 b, d = y, _a(x - a)
1305 return a, b, d, i
1308def _normalTo5(x, y, z, Tun, eps=EPS): # MCCABE 19
1309 '''(INTERNAL) Nearest point on and distance to an I{un-/ordered} triaxial.
1311 @see: I{Eberly}'s U{Distance from a Point to ... an Ellipsoid ...<https://
1312 www.GeometricTools.com/Documentation/DistancePointEllipseEllipsoid.pdf>}.
1313 '''
1314 a, b, c, T = Tun._ordered4
1315 if Tun is not T: # T is ordered, Tun isn't
1316 t = T._order3(x, y, z) + (T,)
1317 a, b, c, d, i = _normalTo5(*t, eps=eps)
1318 return T._order3(a, b, c, reverse=True) + (d, i)
1320 if not (c > 0 and isfinite(a)):
1321 raise _ValueError(a=a, b=b, c=c)
1323 if eps > 0:
1324 val = max(eps * 1e8, EPS)
1325 else: # no validation
1326 val, eps = 0, -eps
1328 i, _a = None, fabs
1329 if z:
1330 if y:
1331 if x:
1332 u = _a(x / a)
1333 v = _a(y / b)
1334 w = _a(z / c)
1335 g = _hypot21(u, v, w)
1336 if g:
1337 r = T._1e2ac # (c / a)**2
1338 s = T._1e2bc # (c / b)**2
1339 t, i = _rootXd(_1_0 / r, _1_0 / s, u, v, w, g, eps)
1340 a = x / (t * r + _1_0)
1341 b = y / (t * s + _1_0)
1342 c = z / (t + _1_0)
1343 d = hypot_(x - a, y - b, z - c)
1344 else: # on the ellipsoid
1345 a, b, c, d = x, y, z, _0_0
1346 else: # x == 0
1347 a = x # 0
1348 b, c, d, i = _normalTo4(y, z, b, c, eps=eps)
1349 elif x: # y == 0
1350 b = y # 0
1351 a, c, d, i = _normalTo4(x, z, a, c, eps=eps)
1352 else: # x == y == 0
1353 if z < 0:
1354 c = -c
1355 a, b, d = x, y, _a(z - c)
1357 else: # z == 0
1358 t = True
1359 d = T._a2c2 # (a + c) * (a - c)
1360 n = a * x
1361 if d > _a(n):
1362 u = n / d
1363 d = T._b2c2 # (b + c) * (b - c)
1364 n = b * y
1365 if d > _a(n):
1366 v = n / d
1367 n = _hypot21(u, v)
1368 if n < 0:
1369 a *= u
1370 b *= v
1371 c *= sqrt(-n)
1372 d = hypot_(x - a, y - b, c)
1373 t = False
1374 if t:
1375 c = z # signed-0
1376 a, b, d, i = _normalTo4(x, y, a, b, eps=eps)
1378 if val > 0: # validate
1379 e = T.sideOf(a, b, c, eps=val)
1380 if e: # not near the ellipsoid's surface
1381 raise _ValueError(a=a, b=b, c=c, d=d,
1382 sideOf=e, eps=val)
1383 if d: # angle of delta and normal vector
1384 m = Vector3d(x, y, z).minus_(a, b, c)
1385 if m.euclid > val:
1386 m = m.unit()
1387 n = T.normal3d(a, b, c)
1388 e = n.dot(m) # n.negate().dot(m)
1389 if not isnear1(_a(e), eps1=val):
1390 raise _ValueError(n=n, m=m,
1391 dot=e, eps=val)
1392 return a, b, c, d, i
1395def _otherV3d_(x_xyz, y, z, **name):
1396 '''(INTERNAL) Get a Vector3d from C{x_xyz}, C{y} and C{z}.
1397 '''
1398 return Vector3d(x_xyz, y, z, **name) if isscalar(x_xyz) else \
1399 _otherV3d(x_xyz=x_xyz)
1402def _rootXd(r, s, u, v, w, g, eps):
1403 '''(INTERNAL) Robust 2d- or 3d-root finder: 2d- if C{s == v == 0} else 3d-root.
1405 @see: I{Eberly}'s U{Robust Root Finders ...<https://www.GeometricTools.com/
1406 Documentation/DistancePointEllipseEllipsoid.pdf>}.
1407 '''
1408 _1, __2 = _1_0, _0_5
1409 _a, _h2 = fabs, _hypot21
1411 u *= r
1412 v *= s # 0 for 2d-root
1413 t0 = w - _1
1414 t1 = _0_0 if g < 0 else _h2(u, w, v)
1415 # assert t0 <= t1
1416 for i in range(1, _TRIPS): # 52-58
1417 e = _a(t0 - t1)
1418 if e < eps:
1419 break
1420 t = (t0 + t1) * __2
1421 if t in (t0, t1):
1422 break
1423 g = _h2(u / (t + r), w / (t + _1),
1424 (v / (t + s)) if v else 0)
1425 if g > 0:
1426 t0 = t
1427 elif g < 0:
1428 t1 = t
1429 else:
1430 break
1431 else: # PYCHOK no cover
1432 t = Fmt.no_convergence(e, eps)
1433 raise _ValueError(t, txt=_rootXd.__name__)
1434 return t, i
1437def _sideOf(xyz, abc, eps=EPS): # in .formy
1438 '''(INTERNAL) Helper for C{_hartzell2}, M{.sideOf} and M{.reverseCartesian}.
1440 @return: M{sum((x / a)**2 for x, a in zip(xyz, abc)) - 1} or C{INT0},
1441 '''
1442 s = _hypot21(*((x / a) for x, a in _zip(xyz, abc) if a)) # strict=True
1443 return s if fabs(s) > eps else INT0
1446if __name__ == '__main__':
1448 from pygeodesy import printf
1449 from pygeodesy.interns import _COMMA_, _NL_, _NLATvar_
1451 # __doc__ of this file, force all into registery
1452 t = [NN] + Triaxials.toRepr(all=True, asorted=True).split(_NL_)
1453 printf(_NLATvar_.join(i.strip(_COMMA_) for i in t))
1455# **) MIT License
1456#
1457# Copyright (C) 2022-2024 -- mrJean1 at Gmail -- All Rights Reserved.
1458#
1459# Permission is hereby granted, free of charge, to any person obtaining a
1460# copy of this software and associated documentation files (the "Software"),
1461# to deal in the Software without restriction, including without limitation
1462# the rights to use, copy, modify, merge, publish, distribute, sublicense,
1463# and/or sell copies of the Software, and to permit persons to whom the
1464# Software is furnished to do so, subject to the following conditions:
1465#
1466# The above copyright notice and this permission notice shall be included
1467# in all copies or substantial portions of the Software.
1468#
1469# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
1470# OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
1471# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
1472# THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR
1473# OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE,
1474# ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR
1475# OTHER DEALINGS IN THE SOFTWARE.