Coverage for pygeodesy/triaxials.py : 96%

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, I{ordered} L{Triaxial} and I{unordered} L{Triaxial_} and miscellaneous classes L{BetaOmega2Tuple}, L{BetaOmega3Tuple}, 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, map1, _zip # from .fsums, .namedTuples, .streprs _0_0, _0_5, _1_0, _N_2_0, isfinite, isnear1, \ _4_0 # PYCHOK used! # from pygeodesy.ellipsoids import Ellipsoid # from .datums # from pygeodesy.elliptic import Elliptic # ._MODS # from pygeodesy.errors import _ValueError # from .streprs _near_, _not_, _NOTEQUAL_, _null_, _opposite_, _outside_, \ _SPACE_, _spherical_, _too_, _x_, _y_ # from pygeodesy.lazily import _ALL_LAZY, _ALL_MODS as _MODS # from .vector3d # from pygeodesy.props import Property_RO # from .fsums
'''(INTERNAL) C{-.toDegrees}, C{-.toRadians} base. ''' isinstance(b, Degrees): return self else:
return self if isinstance(a, Radians) and \ isinstance(b, Radians) else \ self.classof(a.toRadians(), b.toRadians(), *c, name=self.name)
'''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 _ToNamedBase._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 _ToNamedBase._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 _ToNamedBase._toRadians(self, *self)
'''Reduce this L{BetaOmega3Tuple} to a L{BetaOmega2Tuple}. ''' return BetaOmega2Tuple(*self[:2])
'''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 _ToNamedBase._toRadians(self, *self)
'''I{Unordered} triaxial ellipsoid and base class.
Triaxial ellipsoids with right-handed 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, plane 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: I{Geodetic} C{lat}- and C{lon}gitudes are in C{degrees}, I{geodetic} C{phi} and C{lam}bda are in C{radians}, but I{ellipsoidal} lat- and longitude C{beta} and C{omega} are in C{Radians} by default (or in C{Degrees} if converted). '''
'''New I{unordered} L{Triaxial_}.
@arg a_triax: C{X} semi-axis (C{scalar}, conventionally in C{meter}) or an other L{Triaxial} or L{Triaxial_} instance. @kwarg b: C{Y} semi-axis (C{meter}, same units as B{C{a}}), required if C{B{a_triax} is scalar}, ignored otherwise. @kwarg c: C{Z} semi-axis (C{meter}, same units as B{C{a}}), required if C{B{a_triax} is scalar}, ignored otherwise. @kwarg name: Optional name (C{str}).
@raise TriaxialError: Invalid semi-axis or -axes. ''' except (TypeError, ValueError) as x: raise TriaxialError(a=a, b=b, c=c, cause=x)
raise TriaxialError(a=a, b=b, c=c) # txt=_invalid_ raise TriaxialError(a=a, b=b, c=c, txt=_not_ordered_) raise TriaxialError(a=a, c=c, e2ac=self.e2ac, txt=_spherical_)
'''Get the (largest) C{x} semi-axis (C{meter}, conventionally). '''
'''(INTERNAL) Get C{a**2 - b**2} == E_sub_e**2. '''
'''(INTERNAL) Get C{(a/b)**2}. '''
'''(INTERNAL) Get C{a**2 - c**2} == E_sub_x**2. '''
'''Get the surface area (C{meter} I{squared}). ''' Ellipsoid(a, b=c).areax # a == b else: # a == c == b a = Meter2(area=a**2 * PI4)
'''I{Approximate} the surface area (C{meter} I{squared}).
@kwarg p: Exponent (C{scalar} > 0), 1.6 for near-spherical or 1.5849625007 for "near-flat" triaxials.
@see: U{Surface area<https://WikiPedia.org/wiki/Ellipsoid#Approximate_formula>}. ''' a *= a else:
'''Get the (middle) C{y} semi-axis (C{meter}, same units as B{C{a}}). '''
'''(INTERNAL) Get C{b**2 - c**2} == E_sub_y**2. '''
'''Get the (smallest) C{z} semi-axis (C{meter}, same units as B{C{a}}). '''
'''(INTERNAL) Get C{(c/b)**2}. '''
'''Get the C{ab} ellipse' I{(1st) eccentricity squared} (C{scalar}), M{1 - (b/a)**2}. '''
'''(INTERNAL) Get C{1 - e2ab} == C{(b/a)**2}. '''
'''Get the C{bc} ellipse' I{(1st) eccentricity squared} (C{scalar}), M{1 - (c/b)**2}. '''
'''Get the C{ac} ellipse' I{(1st) eccentricity squared} (C{scalar}), M{1 - (c/a)**2}. '''
'''(INTERNAL) Get C{1 - e2ac} == C{(c/a)**2}. '''
'''(INTERNAL) Get class L{Elliptic} once. '''
'''Compute the intersection of this triaxial's surface with a Line-Of-Sight from a Point-Of-View in space.
@see: Function L{pygeodesy.hartzell4} for further details. '''
'''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}), required if B{C{x_xyz}} if C{scalar}. @kwarg z: Z component (C{scalar}), required if B{C{x_xyz}} if C{scalar}. @kwarg normal: If C{True} the projection is perpendicular to (the nearest point on) this triaxial's surface, otherwise the C{radial} line to this triaxial's center (C{bool}). @kwarg eps: Tolerance for root finding and validation (C{scalar}), use a negative value to skip validation.
@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}}, no convergence in root finding or validation failed.
@see: Method L{Ellipsoid.height4} and I{Eberly}'s U{Distance from a Point to ... an Ellipsoid ...<https://www.GeometricTools.com/Documentation/ DistancePointEllipseEllipsoid.pdf>}. '''
else: else: # radially to triaxial's center except Exception as e: raise TriaxialError(x=x, y=y, z=z, cause=e)
'''Is this triaxial I{ordered} and not I{spherical} (C{bool})? '''
'''Is this triaxial I{spherical} (C{Radius} or INT0)? '''
'''Get a 3-D vector at a cartesian on, perpendicular to 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}), required if B{C{x_xyz}} if C{scalar}. @kwarg z: Z component (C{scalar}), required if B{C{x_xyz}} if C{scalar}. @kwarg length: Optional length and in-/outward direction (C{scalar}).
@return: A C{Vector3d(x_, y_, z_)} normalized to B{C{length}}, pointing in- or outward for neg- respectively positive B{C{length}}.
@note: Cartesian location C{(B{x}, B{y}, B{z})} must be on this triaxial's surface, use method L{Triaxial.sideOf} to validate. ''' # n = 2 * (x / a2, y / b2, z / c2) # == 2 * (x, y * a2 / b2, z * a2 / c2) / a2 # iff ordered # == 2 * (x, y / _1e2ab, z / _1e2ac) / a2 # == unit(x, y / _1e2ab, z / _1e2ac).times(length) raise TriaxialError(x=x_xyz, y=y, z=z, txt=_null_)
'''(INTERNAL) Get M{Vector3d((d/a)**2, (d/b)**2, (d/c)**2)}, M{d = max(a, b, c)}. '''
'''(INTERNAL) Normalize C{s} and C{c} iff not already. ''' fabs(_hypot21_(s, c)) > EPS0:
'''(INTERNAL) Un-/Order C{a}, C{b} and C{c}.
@return: 3-Tuple C{(a, b, c)} ordered by or un-ordered (reverse-ordered) C{ijk} if C{B{reverse}=True}. '''
'''(INTERNAL) Un-/Order a C{Vector3d}.
@return: Vector3d(x, y, z) un-/ordered. '''
'''(INTERNAL) Helper for C{_hartzell3d2} and C{_normalTo5}. ''' '''(INTERNAL) Un-Order C{a}, C{b} and C{c}.
@return: 2-Tuple C{((a, b, c), ijk)} with C{a} >= C{b} >= C{c} and C{ijk} a 3-tuple with the initial indices. ''' # reverse (k, j, i) since (a, b, c) is reversed-sorted
'''(INTERNAL) Get the un-/order indices. '''
'''(INTERNAL) I{Unordered} helper for C{.height4}. ''' # <https://WikiPedia.org/wiki/Ellipse#Polar_form_relative_to_focus> # polar form: radius(phi) = a * b / hypot(a * sphi, b * cphi) (a / hypot(cphi, a / b * sphi)) if a < b else a)
'''Is a cartesian above, below or on the surface of this triaxial?
@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}), required if B{C{x_xyz}} if C{scalar}. @kwarg z: Z component (C{scalar}), required if B{C{x_xyz}} if C{scalar}. @kwarg eps: Near surface tolerance(C{scalar}).
@return: C{INT0} if C{(B{x}, B{y}, B{z})} is near this triaxial's surface within tolerance B{C{eps}}, otherwise a neg- or positive C{float} if in- respectively outside this triaxial.
@see: Methods L{Triaxial.height4} and L{Triaxial.normal3d}. '''
'''(INTERNAL) Helper. ''' raise TriaxialError(Fmt.PAREN(sqrt=x))
'''Convert this triaxial to an L{Ellipsoid}, provided C{a == b} or C{b == c}.
@return: An L{Ellipsoid} with north along this C{Z} axis if C{a == b}, this C{Y} axis if C{a == c} or this C{X} axis if C{b == c}.
@raise TriaxialError: This C{a} != C{b}, C{b} != C{c} and C{c} != C{a}.
@see: Method L{Ellipsoid.toTriaxial}. ''' elif a != c: t = _SPACE_(_a_, _NOTEQUAL_, _b_, _NOTEQUAL_, _c_) raise TriaxialError(a=a, b=b, c=c, txt=t)
'''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 triaxial's name.
@return: This C{Triaxial}'s attributes (C{str}). '''
'''Get the volume (C{meter**3}), M{4 / 3 * PI * a * b * c}. '''
'''I{Ordered} triaxial ellipsoid.
@see: L{Triaxial_} for more information. '''
'''New I{ordered} L{Triaxial}.
@arg a_triax: Largest semi-axis (C{scalar}, conventionally in C{meter}) or an other L{Triaxial} or L{Triaxial_} instance. @kwarg b: Middle semi-axis (C{meter}, same units as B{C{a}}), required if C{B{a_triax} is scalar}, ignored otherwise. @kwarg c: Smallest semi-axis (C{meter}, same units as B{C{a}}), required if C{B{a_triax} is scalar}, ignored otherwise. @kwarg name: Optional name (C{str}).
@note: The semi-axes must be ordered as C{B{a} >= B{b} >= B{c} > 0} and must be ellipsoidal, C{B{a} > B{c}}.
@raise TriaxialError: Semi-axes not ordered, spherical or invalid. '''
''' @see: Method C{forwardBetaOmega}. '''
'''Get the surface area (C{meter} I{squared}).
@see: U{Surface area<https://WikiPedia.org/wiki/Ellipsoid#Surface_area>}. ''' else: # a == b > c
'''(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: Method L{Triaxial.reverseBetaOmega} and 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 triaxial'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)} on the surface.
@raise TriaxialError: This triaxial is near-spherical.
@see: Method L{Triaxial.reverseBetaOmega}, 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>}. '''
'''Project a cartesian on this triaxial.
@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}), required if B{C{x_xyz}} if C{scalar}. @kwarg z: Z component (C{scalar}), required if B{C{x_xyz}} if C{scalar}. @kwarg name: Optional name (C{str}). @kwarg normal_eps: Optional keyword arguments C{B{normal}=True} and C{B{eps}=EPS}, see method L{Triaxial.height4}.
@see: Method L{Triaxial.height4} for further information and method L{Triaxial.reverseCartesian} to reverse the projection. '''
'''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: Method L{Triaxial.reverseLatLon} and 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: Method L{Triaxial.reverseLatLon} and 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
'''(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) Convert I{ellipsoidal} lat- and longitude C{beta} and C{omega} to cartesian coordinates I{on the triaxial's surface}, also I{ordered} helper for C{.height4}. '''
'''Convert cartesian to I{ellipsoidal} lat- and longitude, C{beta}, C{omega} and height.
@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}), required if B{C{x_xyz}} if C{scalar}. @kwarg z: Z component (C{scalar}), required if B{C{x_xyz}} if C{scalar}. @kwarg name: Optional name (C{str}).
@return: A L{BetaOmega3Tuple}C{(beta, omega, height)} with C{beta} and C{omega} in C{Radians} and (radial) C{height} in C{meter}, same units as this triaxial's axes.
@see: Methods L{Triaxial.forwardBetaOmega} and L{Triaxial.forwardBetaOmega_} and U{Expressions (21-22)<https://www.Topo.Auth.GR/wp-content/uploads/ sites/111/2021/12/09_Panou.pdf>}. '''
'''"Unproject" a cartesian on to a cartesion I{off} 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}), required if B{C{x_xyz}} if C{scalar}. @kwarg z: Z component (C{scalar}), required if B{C{x_xyz}} if C{scalar}. @arg h: Height above or below this triaxial's surface (C{meter}, same units as the axes). @kwarg normal: If C{True} the height is C{normal} to the surface, otherwise C{radially} to the center of this triaxial (C{bool}). @kwarg eps: Tolerance for surface test (C{scalar}). @kwarg name: Optional name (C{str}).
@return: A L{Vector3Tuple}C{(x, y, z)}.
@raise TrialError: Cartesian C{(x, y, z)} not on this triaxial's surface.
@see: Methods L{Triaxial.forwardCartesian} and L{Triaxial.height4}. ''' if s: # PYCHOK no cover t = _SPACE_((_inside_ if s < 0 else _outside_), self.toRepr()) raise TriaxialError(eps=eps, sideOf=s, x=v.x, y=v.y, z=v.z, txt=t)
'''Convert cartesian to I{geodetic} lat-, longitude and height.
@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}), required if B{C{x_xyz}} if C{scalar}. @kwarg z: Z component (C{scalar}), required if B{C{x_xyz}} if C{scalar}. @kwarg name: Optional name (C{str}).
@return: A L{LatLon3Tuple}C{(lat, lon, height)} with C{lat} and C{lon} in C{degrees} and (radial) C{height} in C{meter}, same units as this triaxial's axes.
@see: Methods L{Triaxial.forwardLatLon} and L{Triaxial.forwardLatLon_} and U{Expressions (4-5)<https://www.Topo.Auth.GR/wp-content/uploads/ sites/111/2021/12/09_Panou.pdf>}. ''' self._1e2bc, # == 1 - e_sub_y**2 _1_0)
'''(INTERNAL) Helper for C{.reverseBetOmg} and C{.reverseLatLon}. '''
'''This is a conformal projection of a triaxial ellipsoid to a plane in which the C{X} and C{Y} grid lines are straight.
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>}, U{Jacobi's conformal projection<https://GeographicLib.SourceForge.io/C++/doc/jacobi.html>} and Jacobi, C. G. J. I{U{Vorlesungen über Dynamik<https://Books.Google.com/books? id=ryEOAAAAQAAJ&pg=PA212>}}, page 212ff, ''' # @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}). '''
'''Get the Jacobi Conformal quadrant size (L{Jacobi2Tuple}C{(x, y)}). ''' Radians(y=self._c2_b2 * self._yE.cPi), name=JacobiConformal.xyQ2.name)
'''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__)
'''(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}). '''
'''Raised for L{Triaxial} issues. '''
'''(INTERNAL) Get the C{items} at the given I{indices}.
@return: C{Type(items[i] for i in indices)} with C{Type = type(items)}, any C{type} having the special method C{__getitem__}. '''
'''(INTERNAL) Hartzell's "Satellite Lin-of-Sight Intersection ...", I{un- or ordered}. '''
raise _ValueError(_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 _ValueError(_opposite_ if max(ux, vy, wz) > 0 else _outside_)
s = fsum_(_1_0, x2 / a2, y2 / b2, z2 / c2, _N_2_0, floats=True) # like _sideOf raise _ValueError(_outside_ if s > 0 else _inside_) raise _ValueError(_too_(_distant_))
'''Compute the intersection of a tri-/biaxial ellipsoid and a Line-Of-Sight from a Point-Of-View outside.
@arg pov: Point-Of-View outside the tri-/biaxial (C{Cartesian}, L{Ecef9Tuple} or L{Vector3d}). @kwarg los: Line-Of-Sight, I{direction} to the tri-/biaxial (L{Vector3d}) or C{None} to point to the tri-/biaxial's center. @kwarg tri_biax: A triaxial (L{Triaxial}, L{Triaxial_}, L{JacobiConformal}) or biaxial ellipsoid (L{Datum}, L{Ellipsoid}, L{Ellipsoid2}, L{a_f2Tuple} or C{scalar} radius in C{meter}). @kwarg name: Optional name (C{str}).
@return: L{Vector4Tuple}C{(x, y, z, h)} on the tri-/biaxial's surface, with C{h} the distance from B{C{pov}} to C{(x, y, z)} along B{C{los}}.
@raise TriaxialError: Null B{C{pov}} or B{C{los}} vector, B{C{pov}} is inside the tri-/biaxial or B{C{los}} points outside the tri-/biaxial or points in an opposite direction.
@raise TypeError: Invalid B{C{pov}} or B{C{los}}.
@see: Function L{pygeodesy.hartzell}, L{pygeodesy.tyr3d} for B{C{los}} and U{I{Satellite Line-of-Sight Intersection with Earth}<https://StephenHartzell. Medium.com/satellite-line-of-sight-intersection-with-earth-d786b4a6a9b6>}. ''' else: D = tri_biax if isinstance(tri_biax, Datum) else \ _spherical_datum(tri_biax, name=hartzell4.__name__) E = D.ellipsoid T = Triaxial_(E.a, E.a, E.b, name=E.name)
except Exception as x: raise TriaxialError(pov=pov, los=los, tri_biax=tri_biax, cause=x)
'''(INTERNAL) Nearest point on and distance to a 2-D ellipse, I{unordered}.
@see: Function C{pygeodesy.ellipsoids._normalTo3} and I{Eberly}'s U{Distance from a Point to ... an Ellipsoid ...<https://www.GeometricTools.com/ Documentation/DistancePointEllipseEllipsoid.pdf>}. ''' # robust root finder else: else: # PYCHOK no cover e = _a(t0 - t1) t = _root2d.__name__ raise _ValueError(Fmt.no_convergence(e, eps), txt=t)
b, a, d, i = _normalTo4(y, x, b, a, eps=eps) return a, b, d, i
raise _ValueError(a=a, b=b)
else: # on the ellipse a, b, d = x, y, _0_0 else: # x == 0
else: # y == 0 if d > fabs(n): # PYCHOK no cover r = n / d a *= r b *= sqrt(_1_0 - r**2) d = hypot(x - a, b) else: a = -a
'''(INTERNAL) Nearest point on and distance to an I{un- or ordered} triaxial.
@see: I{Eberly}'s U{Distance from a Point to ... an Ellipsoid ...<https:// www.GeometricTools.com/Documentation/DistancePointEllipseEllipsoid.pdf>}. ''' # robust root finder else: else: # PYCHOK no cover e = _a(t0 - t1) t = _root3d.__name__ raise _ValueError(Fmt.no_convergence(e, eps), txt=t)
raise _ValueError(a=a, b=b, c=c)
else: # no validation val, eps = 0, -eps
else: # on the ellipsoid a, b, c, d = x, y, z, _0_0 else: # x == 0 else: # x == y == 0
else: # z == 0
raise _ValueError(a=a, b=b, c=c, d=d, sideOf=e, eps=val) raise _ValueError(n=n, m=m, dot=e, eps=val)
'''(INTERNAL) Get a Vector3d from C{x_xyz}, C{y} and C{z}. ''' _otherV3d(x_xyz=x_xyz)
'''(INTERNAL) Helper for C{_hartzell3d2}, M{.sideOf} and M{.reverseCartesian}.
@return: M{sum((x / a)**2 for x, a in zip(xyz, abc)) - 1} or C{INT0}, '''
'''Get C{sin} and C{cos} of C{x} in C{Degrees}, C{Radians} or {radians}. ''' sincos2(x) if isinstance(x, Radians) else sincos2(float(x))) # assume C{radians}
if __name__ == '__main__':
from pygeodesy import printf
# <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), km2m(b), km2m(c), name=n) printf('# %r', t) if n == 'Earth': printf('# %r', JacobiConformal(t.a, t.b, t.c, name=n)) printf('# %r', _WGS84)
# **) 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.
# % python3 -m pygeodesy.triaxials
# 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) # 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) # 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) # JacobiConformal(name='Earth', a=6378173.435, b=6378103.9, c=6356754.399999999, e2ab=0.000021804, e2bc=0.006683418, e2ac=0.006705077, xyQ2=xyQ2(x=1.572084, y=4.249876), volume=1083208241574987694080, area=510065911057441.0625, area_p=510065915922713.6875) # 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) # 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) # 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) # 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) # 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) # 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) # 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) # 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) # 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) |