Coverage for pygeodesy/triaxials.py : 95%

Hot-keys on this page
r m x p toggle line displays
j k next/prev highlighted chunk
0 (zero) top of page
1 (one) first highlighted chunk
# -*- coding: utf-8 -*-
from I{Charles Karney}'s C++ class U{JacobiConformal<https://GeographicLib.SourceForge.io/C++/ doc/classGeographicLib_1_1JacobiConformal.html#details>} to pure Python, class L{Triaxial} and classes L{BetaOmega3Tuple},L{JacobiError}, L{Jacobi2Tuple} and L{TriaxialError}.
@see: U{Geodesics on a triaxial ellipsoid<https://WikiPedia.org/wiki/Geodesics_on_an_ellipsoid# Geodesics_on_a_triaxial_ellipsoid>} and U{Triaxial coordinate systems and their geometrical interpretation<https://www.Topo.Auth.GR/wp-content/uploads/sites/111/2021/12/09_Panou.pdf>}. ''' # make sure int/int division yields float quotient, see .basics
# from pygeodesy.basics import isscalar # from .fsums _1_0, _3_0, _4_0, isfinite # from pygeodesy.errors import _ValueError # from .fsums _outside_, _spherical_, _too_, _x_, _y_ # from pygeodesy.lazily import _ALL_LAZY, _ALL_MODS as _MODS # from .vector3d # from pygeodesy.streprs import Fmt # from .fmath
'''2-Tuple C{(beta, omega)} with I{ellipsoidal} lat- and longitude C{beta} and C{omega} both in C{Radians} (or C{Degrees}). '''
'''Convert this L{BetaOmega2Tuple} to C{Degrees} or C{toDMS}.
@return: L{BetaOmega2Tuple}C{(beta, omega)} with C{beta} and C{omega} both in C{Degrees} or as an L{toDMS} string provided some B{C{toDMS_kwds}} are supplied. ''' return _toDegrees(self, *self, **toDMS_kwds)
'''Convert this L{BetaOmega2Tuple} to C{Radians}.
@return: L{BetaOmega2Tuple}C{(beta, omega)} with C{beta} and C{omega} both in C{Radians}. ''' return _toRadians(self, *self)
'''3-Tuple C{(beta, omega, height)} with I{ellipsoidal} lat- and longitude C{beta} and C{omega} both in C{Radians} (or C{Degrees}) and the C{height}, rather the (signed) I{distance} to the triaxial's surface (measured along the radial line to the triaxial's center) in C{meter}, conventionally. '''
'''Convert this L{BetaOmega3Tuple} to C{Degrees} or C{toDMS}.
@return: L{BetaOmega3Tuple}C{(beta, omega, height)} with C{beta} and C{omega} both in C{Degrees} or as an L{toDMS} string provided some B{C{toDMS_kwds}} are supplied. '''
'''Convert this L{BetaOmega3Tuple} to C{Radians}.
@return: L{BetaOmega3Tuple}C{(beta, omega, height)} with C{beta} and C{omega} both in C{Radians}. ''' return _toRadians(self, *self)
'''Reduce this L{BetaOmega3Tuple} to a L{BetaOmega2Tuple}. ''' return BetaOmega2Tuple(*self[:2])
'''Raised for L{JacobiConformal} issues. '''
'''2-Tuple C{(x, y)} with a Jacobi Conformal C{x} and C{y} projection, both in C{Radians} (or C{Degrees}). '''
'''Convert this L{Jacobi2Tuple} to C{Degrees} or C{toDMS}.
@return: L{Jacobi2Tuple}C{(x, y)} with C{x} and C{y} both in C{Degrees} or as an L{toDMS} string provided some B{C{toDMS_kwds}} are supplied. '''
'''Convert this L{Jacobi2Tuple} to C{Radians}.
@return: L{Jacobi2Tuple}C{(x, y)} with C{x} and C{y} both in C{Radians}. ''' return _toRadians(self, *self)
'''Raised for L{Triaxial} issues. '''
'''Triaxial ellipsoid with semi-axes C{a}, C{b} and C{c}, oriented such that the large principal ellipse C{ab} is the equator I{Z}=0, I{beta}=0, while the small principal ellipse C{ac} is the prime meridian, I{Y}=0, I{omega}=0. The four umbilic points, C{abs}(I{omega}) = C{abs}(I{beta}) = C{PI/2}, lie on the middle principal ellipse C{bc} in plane I{X}=0, I{omega}=C{PI/2}.
@note: Geodetic lat- and longitudes are in C{degrees}, I{ellipsoidal} lat- and longitude C{beta} and C{omega} are in C{Radians} by default (or C{Degrees}. '''
'''New L{Triaxial}.
@arg a: The largest semi-axis (C{meter}, conventionally). @arg b: The middle semi-axis (C{meter}, same units as B{C{a}}). @arg c: The smallest semi-axis (C{meter}, same units as B{C{a}}). @kwarg name: Optional name (C{str}).
@note: The semi-axes must be ordered as C{B{a} >= B{b} >= B{c} > 0} and be ellipsoidal, C{B{a} > B{c}}.
@raise TriaxialError: Semi-axes not ordered, spherical or invalid. '''
raise E(a=a, b=b, c=c, txt=_not_ordered_)
raise E(a=a, c=c, e2ac=self.e2ac, txt=_spherical_)
return self.toStr()
'''Get the largest semi-axis (C{meter}, conventionally). '''
'''Get the surface area (C{meter} I{squared}).
@see: U{Surface area<https://WikiPedia.org/wiki/Ellipsoid#Surface_area>}. '''
'''I{Approximate} the surface area (C{meter} I{squared}).
@kwarg p: Exponent (C{scalar}), 1.6 for near-spherical or 1.5849625007 for "near-flat" triaxials
@see: U{Surface area<https://WikiPedia.org/wiki/Ellipsoid#Approximate_formula>}. '''
'''Get the middle semi-axis (C{meter}, same units as B{C{a}}). '''
'''Get the smallest semi-axis (C{meter}, same units as B{C{a}}). '''
'''Get the C{ab} ellipse' I{(1st) eccentricity squared} (C{scalar}), M{1 - (b / a)**2}. '''
'''(INTERNAL) Get C{1 - e2ab}. '''
'''Get the C{bc} ellipse' I{(1st) eccentricity squared} (C{scalar}), M{1 - (c / b)**2}. '''
'''(INTERNAL) Get C{1 - e2bc}. '''
'''Get the C{ac} ellipse' I{(1st) eccentricity squared} (C{scalar}), M{1 - (c / a)**2}. '''
'''(INTERNAL) Get C{1 - e2ac}. '''
'''(INTERNAL) Helper for C{.forwardBetOmg}. '''
'''(INTERNAL) Helper for C{.forwardBetOmg}. ''' else: x = y = u = _0_0
'''Convert I{ellipsoidal} lat- and longitude C{beta}, C{omega} and height to cartesian.
@arg beta: Ellipsoidal latitude (C{radians} or L{Degrees}). @arg omega: Ellipsoidal longitude (C{radians} or L{Degrees}). @arg height: Height above or below the ellipsoid's surface (C{meter}, same units as this triaxial's C{a}, C{b} and C{c} semi-axes). @kwarg name: Optional name (C{str}).
@return: A L{Vector3Tuple}C{(x, y, z)}.
@see: U{Expressions (23-25)<https://www.Topo.Auth.GR/wp-content/ uploads/sites/111/2021/12/09_Panou.pdf>}. '''
else:
'''Convert I{ellipsoidal} lat- and longitude C{beta} and C{omega} to cartesian coordinates I{on the ellipsoid's surface}.
@arg sbeta: Ellipsoidal latitude C{sin(beta)} (C{scalar}). @arg cbeta: Ellipsoidal latitude C{cos(beta)} (C{scalar}). @arg somega: Ellipsoidal longitude C{sin(omega)} (C{scalar}). @arg comega: Ellipsoidal longitude C{cos(omega)} (C{scalar}). @kwarg name: Optional name (C{str}).
@return: A L{Vector3Tuple}C{(x, y, z)}.
@see: U{Triaxial ellipsoid coordinate system<https://WikiPedia.org/wiki/ Geodesics_on_an_ellipsoid#Triaxial_ellipsoid_coordinate_system>} and U{expressions (23-25)<https://www.Topo.Auth.GR/wp-content/ uploads/sites/111/2021/12/09_Panou.pdf>}. '''
'''Convert I{geodetic} lat-, longitude and heigth to cartesian.
@arg lat: Geodetic latitude (C{degrees}). @arg lon: Geodetic longitude (C{degrees}). @arg height: Height above the ellipsoid (C{meter}, same units as this triaxial's C{a}, C{b} and C{c} axes). @kwarg name: Optional name (C{str}).
@return: A L{Vector3Tuple}C{(x, y, z)}.
@see: U{Expressions (9-11)<https://www.Topo.Auth.GR/wp-content/ uploads/sites/111/2021/12/09_Panou.pdf>}. '''
'''Convert I{geodetic} lat-, longitude and heigth to cartesian.
@arg slat: Geodetic latitude C{sin(lat)} (C{scalar}). @arg clat: Geodetic latitude C{cos(lat)} (C{scalar}). @arg slon: Geodetic longitude C{sin(lon)} (C{scalar}). @arg clon: Geodetic longitude C{cos(lon)} (C{scalar}). @arg height: Height above the ellipsoid (C{meter}, same units as this triaxial's axes C{a}, C{b} and C{c}). @kwarg name: Optional name (C{str}).
@return: A L{Vector3Tuple}C{(x, y, z)}.
@see: U{Expressions (9-11)<https://www.Topo.Auth.GR/wp-content/ uploads/sites/111/2021/12/09_Panou.pdf>}. '''
'''(INTERNAL) Helper for C{.forwardLatLon} and C{.forwardLatLon_}. ''' # 1 - (1 - (c/a)**2) * sa**2 - (1 - (b / a)**2) * ca**2 * sb**2
'''Compute the intersection of a Line-Of-Sight from a Point-Of-View in space with the surface of this triaxial.
@arg pov: Point-Of-View outside this triaxial (C{Cartesian}, L{Ecef9Tuple} or L{Vector3d}). @kwarg los: Line-Of-Sight, I{direction} to this triaxial (L{Vector3d}) or C{None} to point to this triaxial's center. @kwarg name: Optional name (C{str}).
@return: L{Vector4Tuple}C{(x, y, z, h)} on this triaxial's surface.
@raise TriaxialError: Null B{C{pov}} or B{C{los}} vector, B{C{pov}} is inside this triaxial or B{C{los}} points outside this triaxial or points in an opposite direction.
@raise TypeError: Invalid B{C{pov}} or B{C{los}}.
@see: Function L{pygeodesy.tyr3d} for B{C{los}} and Hartzell, S. U{I{Satellite Line-of-Sight Intersection with Earth}<https://StephenHartzell.Medium.com/ satellite-line-of-sight-intersection-with-earth-d786b4a6a9b6>}. ''' return TriaxialError(pov=pov, los=los, triaxial=self, txt=txt)
'''Compute the projection on and the height of a cartesian above or below this triaxial's surface.
@arg x_xyz: X component (C{scalar}) or a cartesian (C{Cartesian}, L{Ecef9Tuple}, L{Vector3d}, L{Vector3Tuple} or L{Vector4Tuple}). @kwarg y: Y component (C{scalar}), needed if B{C{x_xyz}} if C{scalar}. @kwarg z: Z component (C{scalar}), needed if B{C{x_xyz}} if C{scalar}. @kwarg normal: If C{True} the projection is the nearest point on this triaxial's surface, otherwise the intersection of the radial line to the center and this triaxial's surface. @kwarg eps: Tolerance for root finding (C{scalar}).
@return: L{Vector4Tuple}C{(x, y, z, h)} with the cartesian coordinates C{x}, C{y} and C{z} of the projection on or the intersection with and with the height C{h} above or below the triaxial's surface in C{meter}, conventionally.
@raise TriaxialError: Non-cartesian B{C{xyz}}, invalid B{C{eps}} or no convergence in root finding.
@see: U{Eberly<https://www.GeometricTools.com/Documentation/DistancePointEllipseEllipsoid.pdf>}. ''' _otherV3d(x_xyz=x_xyz)
else: # radial to triaxial's center x, y, z = self.forwardBetaOmega_(v.z, hypot(v.x, v.y), v.y, v.x) h, i = v.minus_(x, y, z).length, None
h = -h # below the surface
'''(INTERNAL) Get C{k2} and C{kp2} for C{._xE}, C{._yE} and C{.area}. ''' # k2 = a2b2 / a2c2 * c2_b2 # kp2 = b2c2 / a2c2 * a2_b2 # b2 = b**2 # xE = Elliptic(k2, -a2b2 / b2, kp2, a2_b2) # yE = Elliptic(kp2, +b2c2 / b2, k2, c2_b2) # aE = Elliptic(kp2, 0, k2, 1) self._b2c2 / self._a2c2 * self._a2_b2)
'''(INTERNAL) Normalize C{s} and C{c} iff needed. '''
'''Convert cartesian to I{ellipsoidal} lat- and longitude, C{beta}, C{omega} and height.
@arg x: X coordinate along C{a}-axis (C{meter}, same units as the axes). @arg y: Y coordinate along C{b}-axis (C{meter}, same units as the axes). @arg z: Z coordinate along C{c}-axis (C{meter}, same units as the axes). @kwarg name: Optional name (C{str}).
@return: A L{BetaOmega3Tuple}C{(beta, omega, height)} with C{beta} and C{omega} in C{Radians} and C{height} in C{meter}, same units as this triaxial's axes.
@see: U{Expressions (21-22)<https://www.Topo.Auth.GR/wp-content/uploads/ sites/111/2021/12/09_Panou.pdf>}. '''
'''Convert cartesian to I{geodetic} lat-, longitude and height.
@arg x: X coordinate along C{a}-axis (C{meter}, same units as the axes). @arg y: Y coordinate along C{b}-axis (C{meter}, same units as the axes). @arg z: Z coordinate along C{c}-axis (C{meter}, same units as the axes). @kwarg name: Optional name (C{str}).
@return: A L{LatLon3Tuple}C{(lat, lon, height)} with C{lat} and C{lon} in C{degrees} and C{height} in C{meter}, same units as this triaxial's axes.
@see: U{Expressions (4-5)<https://www.Topo.Auth.GR/wp-content/uploads/ sites/111/2021/12/09_Panou.pdf>}. '''
'''(INTERNAL) Helper for C{.reverseBetOmg} and C{.reverseLatLon}. '''
'''(INTERNAL) Helper. ''' raise self._Error(Fmt.PAREN(sqrt=x))
'''Return this C{Triaxial} as a string.
@kwarg prec: Precision, number of decimal digits (0..9). @kwarg name: Override name (C{str}) or C{None} to exclude this ellipsoid's name.
@return: This C{Triaxial}'s attributes (C{str}). ''' except AttributeError: t = () T.e2ab.name, T.e2bc.name, T.e2ac.name, T.area.name, T.volume.name, *t)
'''Get the volume (C{meter**3}), M{4 / 3 * PI * a * b * c}. '''
'''This is a conformal projection of the ellipsoid to a plane in which the grid lines are straight, see Jacobi, C. G. J. U{I{Vorlesungen über Dynamik} <https://Books.Google.com/books?id=ryEOAAAAQAAJ&pg=PA212>}, page 212ff.
Ellipsoidal coordinates I{beta} and I{omega} are converted to Jacobi Conformal I{y} respectively I{x} separately. Jacobi's coordinates have been multiplied by C{sqrt(B{a}**2 - B{c}**2) / (2 * B{b})} so that the customary results are returned in the case of an ellipsoid of revolution (or a sphere, I{currently not supported}).
Copyright (C) U{Charles Karney<mailto:Charles@Karney.com>} (2014-2020) and licensed under the MIT/X11 License.
@note: This constructor can not be used to specify a sphere.
@see: L{Triaxial}, C++ class U{JacobiConformal<https://GeographicLib.SourceForge.io/ C++/doc/classGeographicLib_1_1JacobiConformal.html#details>} and U{Jacobi's conformal projection<https://GeographicLib.SourceForge.io/C++/doc/jacobi.html>}. '''
# @Property_RO # def ab(self): # '''Get relative magnitude C{ab} (C{None} or C{meter}, same units as B{C{a}}). # ''' # return self._ab
# @Property_RO # def bc(self): # '''Get relative magnitude C{bc} (C{None} or C{meter}, same units as B{C{a}}). # ''' # return self._bc
'''(INTERNAL) Get the x-elliptic function. ''' # -a2b2 / b2 == (b2 - a2) / b2 == 1 - a2 / b2 == 1 - a2_b2
'''Compute a Jacobi Conformal C{x} projection.
@arg omega: Ellipsoidal longitude (C{radians} or L{Degrees}).
@return: The C{x} projection (C{Radians}). '''
'''Compute a Jacobi Conformal C{x} projection.
@arg somega: Ellipsoidal longitude C{sin(omega)} (C{scalar}). @arg comega: Ellipsoidal longitude C{cos(omega)} (C{scalar}).
@return: The C{x} projection (C{Radians}). '''
'''Compute a Jacobi Conformal C{x} and C{y} projection.
@arg beta: Ellipsoidal latitude (C{radians} or L{Degrees}). @arg omega: Ellipsoidal longitude (C{radians} or L{Degrees}). @kwarg name: Optional name (C{str}).
@return: A L{Jacobi2Tuple}C{(x, y)}. ''' name=name or self.xyR2.__name__)
'''Compute a Jacobi Conformal C{x} and C{y} projection.
@arg sbeta: Ellipsoidal latitude C{sin(beta)} (C{scalar}). @arg cbeta: Ellipsoidal latitude C{cos(beta)} (C{scalar}). @arg somega: Ellipsoidal longitude C{sin(omega)} (C{scalar}). @arg comega: Ellipsoidal longitude C{cos(omega)} (C{scalar}). @kwarg name: Optional name (C{str}).
@return: A L{Jacobi2Tuple}C{(x, y)}. ''' self.yR_(sbeta, cbeta), name=name or self.xyR2_.__name__)
'''Get the Jacobi Conformal quadrant size (L{Jacobi2Tuple}C{(x, y)}). ''' Radians(y=self._c2_b2 * self._yE.cPi), name=JacobiConformal.xyQ2.name)
'''(INTERNAL) Get the x-elliptic function. ''' # b2c2 / b2 == (b2 - c2) / b2 == 1 - c2 / b2 == e2bc
'''Compute a Jacobi Conformal C{y} projection.
@arg beta: Ellipsoidal latitude (C{radians} or L{Degrees}).
@return: The C{y} projection (C{Radians}). '''
'''Compute a Jacobi Conformal C{y} projection.
@arg sbeta: Ellipsoidal latitude C{sin(beta)} (C{scalar}). @arg cbeta: Ellipsoidal latitude C{cos(beta)} (C{scalar}).
@return: The C{y} projection (C{Radians}). '''
'''(INTERNAL) Hartzell's "Satellite Lin-of-Sight Intersection ..." ''' else:
raise Error(_near_(_null_))
-w2 * y2, b2 * u2 * q2, -u2 * z2 * p2, ux * wz * 2 * p2, -w2 * x2 * p2, -u2 * y2 * q2, -v2 * x2 * q2, ux * vy * 2 * q2, floats=True) elif r < 0: # LOS pointing away from or missing the triaxial raise Error(_opposite_ if max(ux, vy, wz) > 0 else _outside_)
raise Error(_inside_ if min(x2 - a2, y2 - b2, z2 - c2) < EPS else _outside_) raise Error(_too_(_distant_))
'''(INTERNAL) Nearest point on and distance to a 2-D ellipse.
@see: Function C{.ellipsoids._normalTo3} and U{Eberly<https:// www.GeometricTools.com/Documentation/DistancePointEllipseEllipsoid.pdf>}. ''' # robust root finder else: else: # PYCHOK no cover t = _root2d.__name__ raise TriaxialError(Fmt.no_convergence(e, eps), txt=t)
raise TriaxialError(a=a, b=b, eps=eps)
else: a, b, d = x, y, _0_0 else: # x == 0 a, b, d = _0_0, y, fabs(b - y)
else: # PYCHOK no cover n = a * x d = (a + b) * (a - b) if d > fabs(n): r = n / d a *= r b *= sqrt(_1_0 - r**2) d = hypot(a - x, b) else: a, b, d = x, _0_0, fabs(a - x)
'''(INTERNAL) Nearest point on and distance to a 3-D triaxial.
@see: U{Eberly<https://www.GeometricTools.com/Documentation/ DistancePointEllipseEllipsoid.pdf>}. ''' # robust root finder break else: else: # PYCHOK no cover t = _root3d.__name__ raise TriaxialError(Fmt.no_convergence(e, eps), txt=t)
raise TriaxialError(a=a, b=b, c=c, eps=eps)
else: # on the ellipsoid a, b, c, d = x, y, z, _0_0 else: # x == 0 else: # x == y == 0 a, b, c, d = x, y, z, fabs(c - z)
else: # z == 0 if d > fabs(n): # PYCHOK no cover v = n / d d = _1_0 - hypot2(u, v) if d > 0: a *= u b *= v c *= sqrt(d) d = hypot_(a - x, b - y, c) t = True
raise TriaxialError(x=a, y=b, z=c, d=d, triaxial=T, e=e, eps=eps, txt=_not_('on'))
'''Get C{sin} and C{cos} of C{x} in C{radians} or {degrees}. ''' sincos2(x) if isinstance(x, Radians) else sincos2(float(x))) # assume C{radians}
'''Helper for L{BetaOmega3Tuple} and L{Jacobi2Tuple} ''' isinstance(b, Degrees): return inst else:
'''Helper for L{BetaOmega3Tuple} and L{Jacobi2Tuple} ''' return inst if isinstance(a, Radians) and \ isinstance(b, Radians) else \ inst.classof(a.toRadians(), b.toRadians(), *c, name=inst.name)
if __name__ == '__main__':
from pygeodesy.lazily import printf
def _km2m(*abc): for m in abc: yield m * 1e3
# <https://ArxIV.org/pdf/1909.06452.pdf> Table 1 Semi-axes in km # Planet # <https://www.JPS.NASA.gov/education/images/pdf/ss-moons.pdf> # <https://link.Springer.com/article/10.1007/s00190-022-01650-9> for n, a, b, c in (('Amalthea', 125.0, 73.0, 64), # Jupiter ('Ariel', 581.1, 577.9, 577.7), # Uranus ('Earth', 6378.173435, 6378.1039, 6356.7544), ('Enceladus', 256.6, 251.4, 248.3), # Saturn ('Europa', 1564.13, 1561.23, 1560.93), # Jupiter ('Io', 1829.4, 1819.3, 1815.7), # Jupiter ('Mars', 3394.6, 3393.3, 3376.3), ('Mimas', 207.4, 196.8, 190.6), # Saturn ('Miranda', 240.4, 234.2, 232.9), # Uranus ('Moon', 1735.55, 1735.324, 1734.898), # Earth ('Tethys', 535.6, 528.2, 525.8)): # Saturn t = Triaxial(*_km2m(a, b, c), name=n) printf('%r, area_p=%g', t, t.area_p())
# **) MIT License # # Copyright (C) 2022-2023 -- mrJean1 at Gmail -- All Rights Reserved. # # Permission is hereby granted, free of charge, to any person obtaining a # copy of this software and associated documentation files (the "Software"), # to deal in the Software without restriction, including without limitation # the rights to use, copy, modify, merge, publish, distribute, sublicense, # and/or sell copies of the Software, and to permit persons to whom the # Software is furnished to do so, subject to the following conditions: # # The above copyright notice and this permission notice shall be included # in all copies or substantial portions of the Software. # # THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS # OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, # FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL # THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR # OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, # ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR # OTHER DEALINGS IN THE SOFTWARE. |