Add LVLH frame and geoid intersection method

skyfield/documentation/earth-satellites.rst

@@ -540,6 +540,65 @@ or else in the dynamical coordinate system of the date you specify.
 
 See :doc:`positions` to learn more about these possibilities.
 
+Where on Earth is a satellite looking?
+-----------------------------------------
+
+Rather than knowing where to look to observe a satellite in the sky, you might
+be interested in doing the reverse; finding the geographic location on Earth
+that a satellite would be observing for a given viewing geometry. To do this,
+you can define an attitude vector for the satellite position, and then use
+a :data:`~skyfield.toposlib.wgs84` intersection method to find the surface
+position in the line-of-sight:
+
+* :meth:`~skyfield.toposlib.Geoid._intersection_of()`
+
+For example, `roll`, `pitch`, and `yaw` angles can be provided to the
+:meth:`~skyfield.sgp4lib.EarthSatellite._attitude()` method, which will return
+a :class:`~skyfield.positionlib.ICRF` position vector representing
+the satellite's line-of-sight. The `target` and `center` attributes store the
+attitude vector and the satellite position, respectively;
+
+.. testcode::
+
+    attitude = satellite._attitude(t=t, roll=0.0, pitch=0.0, yaw=0.0)
+    print('Attitude.target:', attitude.target)
+    print('Attitude.center:', attitude.center)
+
+.. testoutput::
+
+    Attitude.target: [ 0.57745914  0.2781511  -0.76757599]
+    Attitude.center: <Geocentric ICRS position and velocity at date t center=399 target=-125544>
+
+In this case we have set the pitch and roll to zero, so the attitude vector
+is pointing toward the center of mass of the Earth, i.e., a unit vector
+anti-aligned with the position vector:
+
+.. testcode::
+
+    from skyfield.functions import length_of
+
+    neg_position = -geocentric.position.au / length_of(geocentric.position.au)
+    print(neg_position)
+
+.. testoutput::
+
+    [ 0.57745914  0.2781511  -0.76757599]
+
+Now, to find the intersection of this vector with the Earth’s surface we can
+use the :meth:`~skyfield.toposlib.Geoid._intersection_of()` method:
+
+.. testcode::
+
+    intersection = wgs84._intersection_of(attitude)
+    print('Latitude:', intersection.latitude)
+    print('Longitude:', intersection.longitude)
+
+.. testoutput::
+
+    Latitude: 50deg 15' 19.6"
+    Longitude: -86deg 23' 23.3"
+
+
 Find a satellite’s range rate
 -----------------------------
 

skyfield/framelib.py

@@ -1,10 +1,10 @@
 # -*- coding: utf-8 -*-
 """Raw transforms between coordinate frames, as NumPy matrices."""
 
-from numpy import array
+from numpy import array, cross
 from .constants import ANGVEL, ASEC2RAD, DAY_S, tau
 from .data.spice import inertial_frames as _inertial_frames
-from .functions import mxm, rot_x, rot_z
+from .functions import length_of, mxm, rot_x, rot_z
 
 def build_matrix():
     # 'xi0', 'eta0', and 'da0' are ICRS frame biases in arcseconds taken
@@ -149,6 +149,34 @@ def __init__(self, doc, matrix):
     def rotation_at(self, t):
         return self._matrix
 
+class LVLH(object):
+    """The Local Vertical Local Horizontal (LVLH) reference frame.
+
+    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, 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)
+        position_vec = position.position.km
+        velocity_vec = position.velocity.km_per_s
+        z_lvlh = -position_vec / length_of(position_vec)
+        L_vec = cross(velocity_vec, z_lvlh)
+        y_lvlh = -L_vec / length_of(L_vec)
+        x_lvlh = cross(y_lvlh, z_lvlh)
+
+        matrix = array([x_lvlh, y_lvlh, z_lvlh]).T
+
+        return matrix
+
 equatorial_B1950_frame = InertialFrame(
     'Reference frame of the Earth’s mean equator and equinox at B1950.',
     _inertial_frames['B1950'],

skyfield/sgp4lib.py

@@ -7,7 +7,9 @@
 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 .positionlib import ICRF
 from .searchlib import _find_discrete, find_maxima
 from .timelib import compute_calendar_date
 from .vectorlib import VectorFunction
@@ -281,6 +283,36 @@ 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 ICRF coordinates.
+
+        Rotations are applied extrinsically relative to the LVLH frame axes,
+        where the neutral attitude vector is pointing along the z-axis.
+        Rotation angles 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": (roll, rot_x),
+            "y": (pitch, rot_y),
+            "z": (yaw, rot_z),
+        }
+
+        lvlh_attitude = array([0, 0, 1])
+        for axis in rotation_order:
+            angle, rotation = rotations[axis]
+            if angle != 0:
+                lvlh_attitude = mxv(rotation(angle), lvlh_attitude)
+
+        R = LVLH(self).rotation_at(t)
+        direction = mxv(R, lvlh_attitude)
+        attitude = ICRF(
+            position_au=self.at(t).xyz.au,
+            velocity_au_per_d=self.at(t).velocity.au_per_d,
+            t=t, center=self.at(t), target=direction
+        )
+
+        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,61 @@ 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.xyz.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)
+        # Angle between attitude and position vector after pitch or roll should
+        # be equal to the provided 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)
+
+        # Yaw without other rotations should not change attitude.
+        a = sat._attitude(t, yaw=angle)
+        assert isclose(angle_between(a.target, -pos), 0.0)
+
+        # Angle between attitude and position vector after yaw with pitch or
+        # or roll should be equal to the provided angle, regardless of order.
+        a = sat._attitude(t, pitch=angle, yaw=angle, rotation_order="yz")
+        assert isclose(angle_between(a.target, -pos), abs_angle)
+
+        a = sat._attitude(t, pitch=angle, yaw=angle, rotation_order="zy")
+        assert isclose(angle_between(a.target, -pos), abs_angle)
+
+        a = sat._attitude(t, roll=angle, yaw=angle, rotation_order="xz")
+        assert isclose(angle_between(a.target, -pos), abs_angle)
+
+        a = sat._attitude(t, roll=angle, yaw=angle, rotation_order="zx")
+        assert isclose(angle_between(a.target, -pos), abs_angle)
+
+        # Rotations applied in reverse will result in different attitudes,
+        # but the angle between attitude and position vector should be equal.
+        axyz = sat._attitude(t, roll=angle, pitch=angle, yaw=angle)
+        azyx = sat._attitude(t, roll=angle, pitch=angle, yaw=angle,
+                             rotation_order="zyx")
+        assert isclose(
+            angle_between(axyz.target, -pos), angle_between(azyx.target, -pos)
+            )

skyfield/tests/test_topos.py

@@ -1,10 +1,11 @@
 from assay import assert_raises
+from collections import namedtuple
 from numpy import abs, arange, sqrt
 
 from skyfield import constants
 from skyfield.api import Distance, load, wgs84, wms
 from skyfield.functions import length_of
-from skyfield.positionlib import Apparent, Barycentric
+from skyfield.positionlib import Apparent, Barycentric, ICRF
 from skyfield.toposlib import ITRSPosition, iers2010
 
 angle = (-15, 15, 35, 45)
@@ -262,3 +263,88 @@ def test_deprecated_position_subpoint_method(ts, angle):
     error_degrees = abs(b.longitude.degrees - angle)
     error_mas = 60.0 * 60.0 * 1000.0 * error_degrees
     assert error_mas < 0.1
+
+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)
+    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
+    assert error_mas < 0.1
+    assert earth_point.elevation.m < 0.1
+
+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)
+    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
+    assert error_mas < 0.1
+
+    error_degrees = abs(earth_point.longitude.degrees - 0.0)
+    error_mas = 60.0 * 60.0 * 1000.0 * error_degrees
+    assert error_mas < 0.1
+    assert earth_point.elevation.m < 0.1
+
+def test_limb_intersection_points(ts):
+    t = ts.utc(2018, 1, 19, 14, 37, 55)
+    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)
+
+    # 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
+    zenith = Vector(pos, [-1.0, 0.0, 0.0], t)
+
+    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
+    zenith_lon = intersection_zenith.longitude.degrees
+    left_limb_lon = zenith_lon - 90.0
+    right_limb_lon = zenith_lon + 90.0
+
+    error_degrees = abs(intersection_bottom.latitude.degrees + 90.0)
+    assert error_degrees < 0.1
+    assert intersection_bottom.elevation.m < 0.1
+
+    error_degrees = abs(intersection_top.latitude.degrees - 90.0)
+    assert error_degrees < 0.1
+    assert intersection_top.elevation.m < 0.1
+
+    error_degrees = abs(intersection_left.latitude.degrees - 0.0)
+    assert error_degrees < 0.1
+
+    error_degrees = abs(intersection_left.longitude.degrees - left_limb_lon)
+    assert error_degrees < 0.1
+    assert intersection_left.elevation.m < 0.1
+
+    error_degrees = abs(intersection_right.latitude.degrees - 0.0)
+    assert error_degrees < 0.1
+
+    error_degrees = abs(intersection_right.longitude.degrees - right_limb_lon)
+    assert error_degrees < 0.1
+    assert intersection_right.elevation.m < 0.1
+