refactor attitude from framelib to sgp4lib

Author: jamesgrimmett

skyfield/framelib.py

@@ -1,12 +1,12 @@
 # -*- coding: utf-8 -*-
 """Raw transforms between coordinate frames, as NumPy matrices."""
 
-from numpy import array, cos, sin
+from numpy import array
 
 from .constants import ANGVEL, ASEC2RAD, DAY_S, pi, tau
 from .data.spice import inertial_frames as _inertial_frames
 from .elementslib import osculating_elements_of
-from .functions import mxm, mxmxm, mxv, rot_x, rot_y, rot_z
+from .functions import mxm, mxmxm, rot_x, rot_y, rot_z
 
 
 def build_matrix():
@@ -155,15 +155,23 @@ def __repr__(self):
     def rotation_at(self, t):
         return self._matrix
 
-class LVLH(InertialFrame):
+class LVLH(object):
     """The Local Vertical Local Horizontal (LVLH) reference frame.
 
-    An inertial frame for a body in orbit, according to CCSDS conventions.
+    A reference frame for a body in orbit, according to CCSDS conventions.
     The frame origin is centered on the body `position`, with z-axis in the
     direction of the `position.center` (i.e., the negative position vector),
     and y-axis opposite the orbital momentum vector.
     """
-    def __init__(self, position):
+    def __init__(self, satellite):
+        self.satellite = satellite
+        self.__doc__ = (
+            "Center ({0}) Pointing Local Vertical Local Horizontal"
+            " reference frame."
+            ).format(satellite.center)
+
+    def rotation_at(self, t):
+        position = self.satellite.at(t)
         elements = osculating_elements_of(position)
         i = elements.inclination.radians
         u = elements.argument_of_latitude.radians
@@ -172,34 +180,7 @@ def __init__(self, position):
         matrix = mxmxm(rot_x(pi / 2.0), rot_y(pi / 2.0), rot_y(om - pi))
         matrix = mxmxm(matrix, rot_z(pi - i), rot_y(-u))
 
-        doc = (
-            "Center ({0}) Pointing Local Vertical Local Horizontal"
-            " reference frame."
-            ).format(position.center)
-        super(LVLH, self).__init__(doc, matrix)
-
-    def local_looking_vector(self, pitch, roll):
-        """Return a unit vector in the attitude direction in local coordinates.
-
-        Rotations should be provided as radians; `roll` is rotation about the
-        x-axis and `pitch` is rotation about the y-axis.
-        """
-        local_looking_vector = array(
-            [sin(pitch), -sin(roll), cos(pitch) * cos(roll)]
-            )
-
-        return local_looking_vector
-
-    def icrf_looking_vector(self, pitch, roll):
-        """Return a unit vector in the attitude direction, in ICRF coordinates.
-
-        Rotations should be provided in radians; `roll` is rotation about the
-        x-axis and `pitch` is rotation about the y-axis.
-        """
-        local_looking_vector = self.local_looking_vector(pitch, roll)
-        icrf_looking_vector = mxv(self._matrix, local_looking_vector)
-
-        return icrf_looking_vector
+        return matrix
 
 equatorial_B1950_frame = InertialFrame(
     'Reference frame of the Earth’s mean equator and equinox at B1950.',

skyfield/sgp4lib.py

@@ -4,9 +4,11 @@
 from numpy import (
     array, concatenate, identity, multiply, ones_like, repeat,
 )
+from collections import namedtuple
 from sgp4.api import SGP4_ERRORS, Satrec
 
 from .constants import AU_KM, DAY_S, T0, tau
+from .framelib import LVLH
 from .functions import _T, mxm, mxv, rot_x, rot_y, rot_z
 from .searchlib import _find_discrete, find_maxima
 from .timelib import compute_calendar_date
@@ -281,6 +283,30 @@ def below_horizon_at(t):
         i = jd.argsort()
         return ts.tt_jd(jd[i]), v[i]
 
+    def _attitude(self, t, roll=0.0, pitch=0.0, yaw=0.0, rotation_order="xyz"):
+        """Return a unit vector in the attitude direction in local coordinates.
+
+        Rotations are applied to an initial vector pointing along the z-axis.
+        Rotations should be provided as radians; `roll` is rotation about the
+        x-axis, `pitch` is about the y-axis, and `yaw` is about the z-axis.
+        """
+        rotations = {
+            "x": rot_x(roll) if roll != 0 else _identity,
+            "y": rot_y(pitch) if pitch != 0 else _identity,
+            "z": rot_z(yaw) if yaw != 0 else _identity,
+        }
+
+        lvlh_attitude = array([0, 0, 1])
+        for axis in rotation_order:
+            lvlh_attitude = mxv(rotations[axis], lvlh_attitude)
+
+        lvlh_frame = LVLH(self)
+        direction = mxv(lvlh_frame.rotation_at(t), lvlh_attitude)
+        Attitude = namedtuple("Attitude", "center, target, t")
+        attitude = Attitude(self.at(t).position, direction, t)
+
+        return attitude
+
 class TEME(object):
     """The SGP4-specific True Equator Mean Equinox frame of reference.
 

skyfield/tests/test_earth_satellites.py

@@ -1,9 +1,10 @@
 # -*- coding: utf-8 -*-
 
-from numpy import array
+from numpy import allclose, array, arange, isclose
 from skyfield import api
 from skyfield.api import EarthSatellite, load
-from skyfield.constants import AU_KM, AU_M
+from skyfield.constants import AU_KM, AU_M, pi
+from skyfield.functions import angle_between, length_of
 from skyfield.sgp4lib import TEME_to_ITRF, VectorFunction
 from skyfield.timelib import julian_date
 
@@ -207,3 +208,41 @@ def _at(self, t):
         True, True, True, False,
         True, True, True, False,  # Last two fake sats can see each other.
     ]
+
+def test_neutral_attitude():
+    ts = api.load.timescale()
+    times = [ts.utc(2023, 4, 8, 1 + i, 25) for i in range(10)]
+    sat = EarthSatellite(line1, line2)
+
+    for t in times:
+        a = EarthSatellite(line1, line2)._attitude(t)
+        p = sat.at(t).position.au
+        a = sat._attitude(t)
+        # Attitude should be the negative of position vector.
+        neg_p = -array(p) / length_of(array(p))
+        assert allclose(a.center.au, p)
+        assert allclose(a.target, neg_p)
+
+def test_attitude():
+    ts = api.load.timescale()
+    t = ts.utc(2023, 4, 8, 1, 25)
+    sat = EarthSatellite(line1, line2)
+    pos = sat.at(t).position.au
+    angles = array(arange(-pi, pi, pi / 8))
+
+    for angle in angles:
+        abs_angle = abs(angle)
+        a = sat._attitude(t, roll=angle)
+        assert isclose(angle_between(a.target, -pos), abs_angle)
+
+        a = sat._attitude(t, pitch=angle)
+        assert isclose(angle_between(a.target, -pos), abs_angle)
+
+        a = sat._attitude(t, yaw=angle)
+        assert isclose(angle_between(a.target, -pos), 0.0)
+
+        a = sat._attitude(t, pitch=angle, yaw=angle)
+        assert allclose(angle_between(a.target, -pos), abs_angle)
+
+        a = sat._attitude(t, roll=angle, yaw=angle)
+        assert allclose(angle_between(a.target, -pos), abs_angle)

skyfield/tests/test_topos.py

@@ -1,4 +1,5 @@
 from assay import assert_raises
+from collections import namedtuple
 from numpy import abs, arange, sqrt
 
 from skyfield import constants
@@ -266,8 +267,10 @@ def test_deprecated_position_subpoint_method(ts, angle):
 def test_intersection_from_pole(ts):
     t = ts.utc(2018, 1, 19, 14, 37, 55)
     p = wgs84.latlon(90.0, 0.0, 1234.0).at(t)
-    attitude = -p.xyz.au / length_of(p.xyz.au)
-    earth_point = wgs84.intersection_of(p, attitude)
+    direction = -p.xyz.au / length_of(p.xyz.au)
+    Vector = namedtuple("Vector", "center, target, t")
+    vector = Vector(p, direction, t)
+    earth_point = wgs84.intersection_of(vector)
 
     error_degrees = abs(earth_point.latitude.degrees - 90.0)
     error_mas = 60.0 * 60.0 * 1000.0 * error_degrees
@@ -277,8 +280,10 @@ def test_intersection_from_pole(ts):
 def test_intersection_from_equator(ts):
     t = ts.utc(2018, 1, 19, 14, 37, 55)
     p = wgs84.latlon(0.0, 0.0, 1234.0).at(t)
-    attitude = -p.xyz.au / length_of(p.xyz.au)
-    earth_point = wgs84.intersection_of(p, attitude)
+    direction = -p.xyz.au / length_of(p.xyz.au)
+    Vector = namedtuple("Vector", "center, target, t")
+    vector = Vector(p, direction, t)
+    earth_point = wgs84.intersection_of(vector)
 
     error_degrees = abs(earth_point.latitude.degrees - 0.0)
     error_mas = 60.0 * 60.0 * 1000.0 * error_degrees
@@ -294,22 +299,26 @@ def test_limb_intersection_points(ts):
     d = 100.0
     a = wgs84.radius.au
     c = a * (1.0 - 1.0 / wgs84.inverse_flattening)
+    pos = ICRF(position_au=[d, 0.0, 0.0], t=t, center=399)
 
-    # Attitude vectors pointing to the polar and equatorial limbs of the Earth
-    attitude_bottom_tangent = [-d, 0.0, -c] / sqrt(d**2 + c**2)
-    attitude_top_tangent = [-d, 0.0, c] / sqrt(d**2 + c**2)
-    attitude_left_tangent = [-d, -a, 0.0] / sqrt(d**2 + c**2)
-    attitude_right_tangent = [-d, a, 0.0] / sqrt(d**2 + c**2)
+    # Vectors pointing to the polar and equatorial limbs of the Earth
+    direction_bottom_tangent = [-d, 0.0, -c] / sqrt(d**2 + c**2)
+    direction_top_tangent = [-d, 0.0, c] / sqrt(d**2 + c**2)
+    direction_left_tangent = [-d, -a, 0.0] / sqrt(d**2 + c**2)
+    direction_right_tangent = [-d, a, 0.0] / sqrt(d**2 + c**2)
+    Vector = namedtuple("Vector", "center, target, t")
+    bottom_tangent = Vector(pos, direction_bottom_tangent, t)
+    top_tangent = Vector(pos, direction_top_tangent, t)
+    left_tangent = Vector(pos, direction_left_tangent, t)
+    right_tangent = Vector(pos, direction_right_tangent, t)
     # Attitude vector pointing straight down
-    attitude_zenith = [-1.0, 0.0, 0.0]
-
-    pos = ICRF(position_au=[d, 0.0, 0.0], t=t, center=399)
+    zenith = Vector(pos, [-1.0, 0.0, 0.0], t)
 
-    intersection_bottom = wgs84.intersection_of(pos, attitude_bottom_tangent)
-    intersection_top = wgs84.intersection_of(pos, attitude_top_tangent)
-    intersection_left = wgs84.intersection_of(pos, attitude_left_tangent)
-    intersection_right = wgs84.intersection_of(pos, attitude_right_tangent)
-    intersection_zenith = wgs84.intersection_of(pos, attitude_zenith)
+    intersection_bottom = wgs84.intersection_of(bottom_tangent)
+    intersection_top = wgs84.intersection_of(top_tangent)
+    intersection_left = wgs84.intersection_of(left_tangent)
+    intersection_right = wgs84.intersection_of(right_tangent)
+    intersection_zenith = wgs84.intersection_of(zenith)
 
     # Viewed from sufficient distance, points of intersection should be nearly
     # tangent to the north and south poles, and the zenith longitude +/- 90.0

skyfield/toposlib.py

@@ -277,32 +277,18 @@ def _compute_latitude(self, position):
 
     subpoint = geographic_position_of  # deprecated method name
 
-    def intersection_of(self, position, direction):
+    def intersection_of(self, vector):
         """Return the surface point intersected by a vector.
 
-        The ``position`` and ``direction`` describe the tail and orientation
-        of the vector to intersect the ellipsoid, respectively. The vector must
-        be provided in ICRF coordinates with a ``.center``` of 399, the center
-        of the Earth. The direction should be an |xyz| unit vector. Returns a
-        `GeographicPosition` giving the geodetic ``latitude`` and ``longitude``
-        at the point that the ray intersects the surface of the Earth.
+        TODO ...
+        Returns a `GeographicPosition` giving the geodetic ``latitude`` and
+        ``longitude`` at the point that the ray intersects the surface of the
+        Earth.
 
         The main calculation implemented here is based on JPL's NAIF toolkit;
         https://naif.jpl.nasa.gov/pub/naif/toolkit_docs/C/req/ellipses.html
         """
-        if not isinstance(position, ICRF):
-            raise ValueError(
-                'you can only calculate an intersection for a vector that is'
-                ' defined by coordinates along the ICRF axes'
-            )
-        if position.center != 399:
-            raise ValueError(
-                'you can only calculate an intersection for a vector that is'
-                ' defined relative to a position that is geocentric'
-                ' (center=399), but this position has a center of {0}'
-                .format(position.center)
-            )
-        if position.t is None:
+        if vector.t is None:
             raise ValueError(
                 'you can only calculate an intersection for a vector that is'
                 ' defined at a fixed point in time (must have a non-null `.t`)'
@@ -317,9 +303,9 @@ def intersection_of(self, position, direction):
         inv_d_matrix = array([[a, 0.0, 0.0], [0.0, a, 0.0], [0.0, 0.0, c]])
 
         # Apply distortion matrix to the vector
-        position_xyz = position.xyz.au
+        position_xyz = vector.center.xyz.au
         d_position_xyz = mxv(d_matrix, position_xyz)
-        d_direction = mxv(d_matrix, direction)
+        d_direction = mxv(d_matrix, vector.target)
         # Rescale the vector to unit length
         d_direction = d_direction / length_of(d_direction)
 
@@ -341,7 +327,7 @@ def intersection_of(self, position, direction):
 
         # Apply inverse distortion and convert back to geocentric coords
         intersection = Geocentric(
-            mxv(inv_d_matrix, d_intersection), t=position.t
+            mxv(inv_d_matrix, d_intersection), t=vector.t
             )
 
         # Retrieve latitude and longitude