diff --git a/skyfield/documentation/earth-satellites.rst b/skyfield/documentation/earth-satellites.rst index 0c6729450..1d4e703b0 100644 --- a/skyfield/documentation/earth-satellites.rst +++ b/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: + +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 ----------------------------- diff --git a/skyfield/framelib.py b/skyfield/framelib.py index 6d3f3d7b0..b7d40386b 100644 --- a/skyfield/framelib.py +++ b/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'], diff --git a/skyfield/sgp4lib.py b/skyfield/sgp4lib.py index b96a87e06..d500ae749 100644 --- a/skyfield/sgp4lib.py +++ b/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. diff --git a/skyfield/tests/test_earth_satellites.py b/skyfield/tests/test_earth_satellites.py index d524197ae..4016c1360 100644 --- a/skyfield/tests/test_earth_satellites.py +++ b/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) + ) \ No newline at end of file diff --git a/skyfield/tests/test_topos.py b/skyfield/tests/test_topos.py index 369dd4655..233bdfe95 100644 --- a/skyfield/tests/test_topos.py +++ b/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 + diff --git a/skyfield/toposlib.py b/skyfield/toposlib.py index 19d236d0d..ce37f7c95 100644 --- a/skyfield/toposlib.py +++ b/skyfield/toposlib.py @@ -5,8 +5,9 @@ from .earthlib import refract from .framelib import itrs from .functions import ( - _T, angular_velocity_matrix, mxm, mxv, rot_y, rot_z, + _T, angular_velocity_matrix, dots, length_of, mxm, mxv, rot_y, rot_z, ) +from .positionlib import Geocentric from .descriptorlib import reify from .units import Angle, Distance, _ltude from .vectorlib import VectorFunction @@ -276,6 +277,66 @@ def _compute_latitude(self, position): subpoint = geographic_position_of # deprecated method name + def _intersection_of(self, icrf_vector): + """Return the surface point intersected by a vector. + + This method calculates the point at which the input `icrf_vector` + intersects the surface of the Geoid. The input `icrf_vector` should be + an :class:`~skyfield.positionlib.ICRF` object (or subclass thereof) + and is implicitly assumed to extend to infinity. + 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 + + Returns a `GeographicPosition` with the geodetic ``latitude`` and + ``longitude`` at the point of intersection. + """ + if icrf_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`)' + ) + + # Define a distortion matrix to map the ellipsoid into a unit sphere + a = self.radius.au + c = a * (1.0 - 1.0 / self.inverse_flattening) + d_matrix = array( + [[1.0 / a, 0.0, 0.0], [0.0, 1.0 / a, 0.0], [0.0, 0.0, 1.0 / c]] + ) + 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 = icrf_vector.center.xyz.au + d_position_xyz = mxv(d_matrix, position_xyz) + d_direction = mxv(d_matrix, icrf_vector.target) + # Rescale the vector to unit length + d_direction = d_direction / length_of(d_direction) + + # The disciminant indicates the number of solutions (intersections) + discriminant = dots(d_direction, d_position_xyz) ** 2 - ( + length_of(d_position_xyz) ** 2 - 1.0 + ) + if discriminant < 0: + raise ValueError( + "The viewing angle does not intersect with the ellipsoid." + ) + # There will be another solution with `+ np.sqrt(discrimant)` (the + # intersection as the ray exits the opposite side of the ellipsoid), + # but we only need the solution on the near side. + distance = -dots(d_direction, d_position_xyz) - sqrt(discriminant) + + # Solve for the location of intersection + d_intersection = d_position_xyz + distance * d_direction + + # Apply inverse distortion and convert back to geocentric coords + intersection = Geocentric( + mxv(inv_d_matrix, d_intersection), t=icrf_vector.t + ) + + # Retrieve latitude and longitude + intersection_latlon = self.geographic_position_of(intersection) + + return intersection_latlon + wgs84 = Geoid('WGS84', 6378137.0, 298.257223563) iers2010 = Geoid('IERS2010', 6378136.6, 298.25642)