From 38971e4b4616acbff02fe7215becb23c80254046 Mon Sep 17 00:00:00 2001 From: James Grimmett Date: Sat, 25 Mar 2023 16:02:47 +1100 Subject: [PATCH 01/14] add vector intersection methods --- skyfield/tests/test_topos.py | 73 +++++++++++++++++++++++++++++++++++- skyfield/toposlib.py | 73 +++++++++++++++++++++++++++++++++++- 2 files changed, 144 insertions(+), 2 deletions(-) diff --git a/skyfield/tests/test_topos.py b/skyfield/tests/test_topos.py index 369dd4655..e223c31eb 100644 --- a/skyfield/tests/test_topos.py +++ b/skyfield/tests/test_topos.py @@ -4,7 +4,7 @@ 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 +262,74 @@ 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) + attitude = -p.xyz.au / length_of(p.xyz.au) + earth_point = wgs84.intersection_of(p, attitude) + + error_degrees = abs(earth_point.latitude.degrees - 90.0) + error_mas = 60.0 * 60.0 * 1000.0 * error_degrees + assert error_mas < 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) + attitude = -p.xyz.au / length_of(p.xyz.au) + earth_point = wgs84.intersection_of(p, attitude) + + 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 + +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) + + # 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] / np.sqrt(d**2 + c**2) + # 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) + + 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) + + # 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 + + error_degrees = abs(intersection_top.latitude.degrees - 90.0) + assert error_degrees < 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 + + 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 + diff --git a/skyfield/toposlib.py b/skyfield/toposlib.py index 19d236d0d..472d1b4d8 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 ICRF, Geocentric from .descriptorlib import reify from .units import Angle, Distance, _ltude from .vectorlib import VectorFunction @@ -276,6 +277,76 @@ def _compute_latitude(self, position): subpoint = geographic_position_of # deprecated method name + def intersection_of(self, position, direction): + """Return the surface point intersected by a vector. + + The position must be provided in ICRF coordinates with a ``.center`` of + 399, the center of the Earth. 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: + 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 = position.xyz.au + d_position_xyz = mxv(d_matrix, position_xyz) + d_direction = mxv(d_matrix, direction) + # 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=position.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) From d0b59d8257e5aa1b51d343d1f14032550c0c642f Mon Sep 17 00:00:00 2001 From: James Grimmett Date: Sat, 25 Mar 2023 16:05:09 +1100 Subject: [PATCH 02/14] add lvlh inertial frame --- skyfield/framelib.py | 55 +++++++++++++++++++++++++++++++++++++++++--- skyfield/toposlib.py | 12 ++++++---- 2 files changed, 59 insertions(+), 8 deletions(-) diff --git a/skyfield/framelib.py b/skyfield/framelib.py index 6d3f3d7b0..3e2a43e33 100644 --- a/skyfield/framelib.py +++ b/skyfield/framelib.py @@ -1,10 +1,13 @@ # -*- coding: utf-8 -*- """Raw transforms between coordinate frames, as NumPy matrices.""" -from numpy import array -from .constants import ANGVEL, ASEC2RAD, DAY_S, tau +from numpy import array, cos, sin + +from .constants import ANGVEL, ASEC2RAD, DAY_S, pi, tau from .data.spice import inertial_frames as _inertial_frames -from .functions import mxm, rot_x, rot_z +from .elementslib import osculating_elements_of +from .functions import mxm, mxmxm, mxv, rot_x, rot_y, rot_z + def build_matrix(): # 'xi0', 'eta0', and 'da0' are ICRS frame biases in arcseconds taken @@ -149,6 +152,52 @@ def __init__(self, doc, matrix): def rotation_at(self, t): return self._matrix +class LVLH(InertialFrame): + """The Local Vertical Local Horizontal (LVLH) reference frame. + + An inertial 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): + elements = osculating_elements_of(position) + i = elements.inclination.radians + u = elements.argument_of_latitude.radians + om = elements.longitude_of_ascending_node.radians + + 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 + equatorial_B1950_frame = InertialFrame( 'Reference frame of the Earth’s mean equator and equinox at B1950.', _inertial_frames['B1950'], diff --git a/skyfield/toposlib.py b/skyfield/toposlib.py index 472d1b4d8..efd4f6202 100644 --- a/skyfield/toposlib.py +++ b/skyfield/toposlib.py @@ -280,11 +280,13 @@ def _compute_latitude(self, position): def intersection_of(self, position, direction): """Return the surface point intersected by a vector. - The position must be provided in ICRF coordinates with a ``.center`` of - 399, the center of the Earth. Returns a `GeographicPosition` giving the - geodetic ``latitude`` and ``longitude`` at the point that the ray - intersects the surface of the Earth. - + 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. + 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 """ From fcf15ae0a6e16b8adedd9d8ef190210fdeae2f73 Mon Sep 17 00:00:00 2001 From: James Grimmett Date: Mon, 27 Mar 2023 20:18:39 +1100 Subject: [PATCH 03/14] add geoid intersection example to docs --- skyfield/documentation/earth-satellites.rst | 57 +++++++++++++++++++++ 1 file changed, 57 insertions(+) diff --git a/skyfield/documentation/earth-satellites.rst b/skyfield/documentation/earth-satellites.rst index 0c6729450..a267fa4e6 100644 --- a/skyfield/documentation/earth-satellites.rst +++ b/skyfield/documentation/earth-satellites.rst @@ -540,6 +540,63 @@ 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, first contruct an inertial :data:`~skyfield.framelib.LVLH` frame +for the satellite position, then use the frame to define an attitude vector: + +.. testcode:: + + from skyfield.framelib import LVLH + + satellite_lvlh = LVLH(geocentric) + attitude = satellite_lvlh.icrf_looking_vector(pitch=0.0, roll=0.0) + print('Frame:', satellite_lvlh) + print('Attitude:', attitude) + +.. testoutput:: + + Frame: + Attitude: [ 0.57745914 0.2781511 -0.76757599] + +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:: + + neg_position = -geocentric.position.au / length_of(geocentric.position.au) + print(negative_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(geocentric, attitude) + print('Latitude:', intersection.latitude) + print('Longitude:', intersection.longitude) + +.. testoutput:: + + Longitude: -86deg 23' 23.3" + Latitude: 50deg 15' 19.6" + + Find a satellite’s range rate ----------------------------- From 96316b3addb299c9b538354c73b112199371642f Mon Sep 17 00:00:00 2001 From: James Grimmett Date: Tue, 28 Mar 2023 09:38:26 +1100 Subject: [PATCH 04/14] fix typo in example docs --- skyfield/documentation/earth-satellites.rst | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/skyfield/documentation/earth-satellites.rst b/skyfield/documentation/earth-satellites.rst index a267fa4e6..fbacb75c7 100644 --- a/skyfield/documentation/earth-satellites.rst +++ b/skyfield/documentation/earth-satellites.rst @@ -576,11 +576,11 @@ anti-aligned with the position vector: .. testcode:: neg_position = -geocentric.position.au / length_of(geocentric.position.au) - print(negative_position) + print(neg_position) .. testoutput:: - [-0.57745914 -0.2781511 0.76757599] + [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: @@ -593,8 +593,8 @@ use the :meth:`~skyfield.toposlib.Geoid.intersection_of()` method: .. testoutput:: - Longitude: -86deg 23' 23.3" Latitude: 50deg 15' 19.6" + Longitude: -86deg 23' 23.3" Find a satellite’s range rate From eef729f71915be9a05a3a9d66dd1ebbeea748c2d Mon Sep 17 00:00:00 2001 From: James Grimmett Date: Sun, 2 Apr 2023 15:11:28 +1000 Subject: [PATCH 05/14] add elevation check to intersection tests --- skyfield/tests/test_topos.py | 8 +++++++- 1 file changed, 7 insertions(+), 1 deletion(-) diff --git a/skyfield/tests/test_topos.py b/skyfield/tests/test_topos.py index e223c31eb..870e176cb 100644 --- a/skyfield/tests/test_topos.py +++ b/skyfield/tests/test_topos.py @@ -272,6 +272,7 @@ def test_intersection_from_pole(ts): 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) @@ -286,6 +287,7 @@ def test_intersection_from_equator(ts): 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) @@ -297,7 +299,7 @@ def test_limb_intersection_points(ts): 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] / np.sqrt(d**2 + c**2) + attitude_right_tangent = [-d, a, 0.0] / sqrt(d**2 + c**2) # Attitude vector pointing straight down attitude_zenith = [-1.0, 0.0, 0.0] @@ -317,19 +319,23 @@ def test_limb_intersection_points(ts): 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 From 836a366951c1bc3d4edf15f7276ae042296232b8 Mon Sep 17 00:00:00 2001 From: James Date: Sat, 8 Apr 2023 13:15:33 +1000 Subject: [PATCH 06/14] add a custom __repr__ method for inertial frame objects --- skyfield/documentation/earth-satellites.rst | 6 ++++-- skyfield/framelib.py | 5 ++++- skyfield/toposlib.py | 2 +- 3 files changed, 9 insertions(+), 4 deletions(-) diff --git a/skyfield/documentation/earth-satellites.rst b/skyfield/documentation/earth-satellites.rst index fbacb75c7..9e1cab703 100644 --- a/skyfield/documentation/earth-satellites.rst +++ b/skyfield/documentation/earth-satellites.rst @@ -566,7 +566,7 @@ for the satellite position, then use the frame to define an attitude vector: .. testoutput:: - Frame: + Frame: Center (399) Pointing Local Vertical Local Horizontal reference frame. Attitude: [ 0.57745914 0.2781511 -0.76757599] In this case we have set the pitch and roll to zero, so the attitude vector @@ -575,12 +575,14 @@ 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] + [ 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: diff --git a/skyfield/framelib.py b/skyfield/framelib.py index 3e2a43e33..fe1f6189f 100644 --- a/skyfield/framelib.py +++ b/skyfield/framelib.py @@ -149,6 +149,9 @@ def __init__(self, doc, matrix): self.__doc__ = doc self._matrix = matrix + def __repr__(self): + return '<{0}> {1}'.format(self.__class__.__name__, self.__doc__) + def rotation_at(self, t): return self._matrix @@ -170,7 +173,7 @@ def __init__(self, position): matrix = mxmxm(matrix, rot_z(pi - i), rot_y(-u)) doc = ( - "Center {0} Pointing Local Vertical Local Horizontal" + "Center ({0}) Pointing Local Vertical Local Horizontal" " reference frame." ).format(position.center) super(LVLH, self).__init__(doc, matrix) diff --git a/skyfield/toposlib.py b/skyfield/toposlib.py index efd4f6202..9a5042332 100644 --- a/skyfield/toposlib.py +++ b/skyfield/toposlib.py @@ -286,7 +286,7 @@ def intersection_of(self, position, direction): 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. - + 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 """ From d1465e9604cae9c4efdb92aa19ad17caf715a147 Mon Sep 17 00:00:00 2001 From: James Date: Sun, 20 Aug 2023 17:32:29 +1000 Subject: [PATCH 07/14] refactor attitude from framelib to sgp4lib --- skyfield/framelib.py | 47 ++++++++----------------- skyfield/sgp4lib.py | 26 ++++++++++++++ skyfield/tests/test_earth_satellites.py | 43 ++++++++++++++++++++-- skyfield/tests/test_topos.py | 43 +++++++++++++--------- skyfield/toposlib.py | 32 +++++------------ 5 files changed, 116 insertions(+), 75 deletions(-) diff --git a/skyfield/framelib.py b/skyfield/framelib.py index fe1f6189f..38ddf0b07 100644 --- a/skyfield/framelib.py +++ b/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.', diff --git a/skyfield/sgp4lib.py b/skyfield/sgp4lib.py index b96a87e06..954915aff 100644 --- a/skyfield/sgp4lib.py +++ b/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. diff --git a/skyfield/tests/test_earth_satellites.py b/skyfield/tests/test_earth_satellites.py index d524197ae..4810cc551 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,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) \ No newline at end of file diff --git a/skyfield/tests/test_topos.py b/skyfield/tests/test_topos.py index 870e176cb..58534649a 100644 --- a/skyfield/tests/test_topos.py +++ b/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 diff --git a/skyfield/toposlib.py b/skyfield/toposlib.py index 9a5042332..695abee91 100644 --- a/skyfield/toposlib.py +++ b/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 From 21adbf8955f56701a7749503dfb54b1a9937e24a Mon Sep 17 00:00:00 2001 From: James Date: Tue, 14 Nov 2023 15:48:27 +1100 Subject: [PATCH 08/14] simplify calculation of lvlh frame rotation --- skyfield/framelib.py | 23 +++++++++++------------ 1 file changed, 11 insertions(+), 12 deletions(-) diff --git a/skyfield/framelib.py b/skyfield/framelib.py index 38ddf0b07..4d879682d 100644 --- a/skyfield/framelib.py +++ b/skyfield/framelib.py @@ -1,12 +1,10 @@ # -*- coding: utf-8 -*- """Raw transforms between coordinate frames, as NumPy matrices.""" -from numpy import array - -from .constants import ANGVEL, ASEC2RAD, DAY_S, pi, tau +from numpy import array, cross +from .constants import ANGVEL, ASEC2RAD, DAY_S, tau from .data.spice import inertial_frames as _inertial_frames -from .elementslib import osculating_elements_of -from .functions import mxm, mxmxm, rot_x, rot_y, rot_z +from .functions import length_of, mxm, rot_x, rot_z def build_matrix(): @@ -172,13 +170,14 @@ def __init__(self, satellite): 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 - om = elements.longitude_of_ascending_node.radians - - 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)) + 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 From 2c3adf997571cfbb687f2608680fdaf5f5f6a5b8 Mon Sep 17 00:00:00 2001 From: James Date: Tue, 14 Nov 2023 18:09:58 +1100 Subject: [PATCH 09/14] reduce calculations in attitude calculation and represent result using Geometric object --- skyfield/sgp4lib.py | 34 ++++++++++++++++++++-------------- 1 file changed, 20 insertions(+), 14 deletions(-) diff --git a/skyfield/sgp4lib.py b/skyfield/sgp4lib.py index 954915aff..7b0ce8e33 100644 --- a/skyfield/sgp4lib.py +++ b/skyfield/sgp4lib.py @@ -4,12 +4,12 @@ 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 .positionlib import Geometric from .searchlib import _find_discrete, find_maxima from .timelib import compute_calendar_date from .vectorlib import VectorFunction @@ -284,26 +284,32 @@ def below_horizon_at(t): 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. + """Return a unit vector in the attitude direction in ICRF 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 are applied 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": 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, + "x": (roll, rot_x), + "y": (pitch, rot_y), + "z": (yaw, rot_z), } 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) + 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 = Geometric( + 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 From 4bdc3bc7947cbea2e1dbbb8f6ea70c3e1c359c1e Mon Sep 17 00:00:00 2001 From: James Date: Tue, 14 Nov 2023 18:10:38 +1100 Subject: [PATCH 10/14] add docstring and use explicit variable names --- skyfield/toposlib.py | 24 +++++++++++++----------- 1 file changed, 13 insertions(+), 11 deletions(-) diff --git a/skyfield/toposlib.py b/skyfield/toposlib.py index 695abee91..5a0b027b9 100644 --- a/skyfield/toposlib.py +++ b/skyfield/toposlib.py @@ -7,7 +7,7 @@ from .functions import ( _T, angular_velocity_matrix, dots, length_of, mxm, mxv, rot_y, rot_z, ) -from .positionlib import ICRF, Geocentric +from .positionlib import Geocentric from .descriptorlib import reify from .units import Angle, Distance, _ltude from .vectorlib import VectorFunction @@ -277,18 +277,20 @@ def _compute_latitude(self, position): subpoint = geographic_position_of # deprecated method name - def intersection_of(self, vector): + def intersection_of(self, icrf_vector): """Return the surface point intersected by a vector. - TODO ... - Returns a `GeographicPosition` giving the geodetic ``latitude`` and - ``longitude`` at the point that the ray intersects the surface of the - Earth. - + 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 vector.t is None: + 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`)' @@ -303,9 +305,9 @@ def intersection_of(self, vector): 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 = vector.center.xyz.au + position_xyz = icrf_vector.center.xyz.au d_position_xyz = mxv(d_matrix, position_xyz) - d_direction = mxv(d_matrix, vector.target) + d_direction = mxv(d_matrix, icrf_vector.target) # Rescale the vector to unit length d_direction = d_direction / length_of(d_direction) @@ -327,7 +329,7 @@ def intersection_of(self, vector): # Apply inverse distortion and convert back to geocentric coords intersection = Geocentric( - mxv(inv_d_matrix, d_intersection), t=vector.t + mxv(inv_d_matrix, d_intersection), t=icrf_vector.t ) # Retrieve latitude and longitude From 5bd963d8849b04a8a152ce86cf004953cd9e0488 Mon Sep 17 00:00:00 2001 From: James Date: Tue, 14 Nov 2023 18:11:34 +1100 Subject: [PATCH 11/14] update unit tests and documentation --- skyfield/documentation/earth-satellites.rst | 22 +++++++-------- skyfield/tests/test_earth_satellites.py | 30 +++++++++++++++++---- 2 files changed, 36 insertions(+), 16 deletions(-) diff --git a/skyfield/documentation/earth-satellites.rst b/skyfield/documentation/earth-satellites.rst index 9e1cab703..b5eda3d2c 100644 --- a/skyfield/documentation/earth-satellites.rst +++ b/skyfield/documentation/earth-satellites.rst @@ -552,22 +552,22 @@ position in the line-of-sight: * :meth:`~skyfield.toposlib.Geoid.intersection_of()` -For example, first contruct an inertial :data:`~skyfield.framelib.LVLH` frame -for the satellite position, then use the frame to define an attitude vector: +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:: - from skyfield.framelib import LVLH - - satellite_lvlh = LVLH(geocentric) - attitude = satellite_lvlh.icrf_looking_vector(pitch=0.0, roll=0.0) - print('Frame:', satellite_lvlh) - print('Attitude:', attitude) + 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:: - Frame: Center (399) Pointing Local Vertical Local Horizontal reference frame. - Attitude: [ 0.57745914 0.2781511 -0.76757599] + 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 @@ -589,7 +589,7 @@ use the :meth:`~skyfield.toposlib.Geoid.intersection_of()` method: .. testcode:: - intersection = wgs84.intersection_of(geocentric, attitude) + intersection = wgs84.intersection_of(attitude) print('Latitude:', intersection.latitude) print('Longitude:', intersection.longitude) diff --git a/skyfield/tests/test_earth_satellites.py b/skyfield/tests/test_earth_satellites.py index 4810cc551..4016c1360 100644 --- a/skyfield/tests/test_earth_satellites.py +++ b/skyfield/tests/test_earth_satellites.py @@ -220,7 +220,7 @@ def test_neutral_attitude(): 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.center.xyz.au, p) assert allclose(a.target, neg_p) def test_attitude(): @@ -232,17 +232,37 @@ def test_attitude(): 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) - a = sat._attitude(t, pitch=angle, yaw=angle) - assert allclose(angle_between(a.target, -pos), abs_angle) + # 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) - a = sat._attitude(t, roll=angle, yaw=angle) - assert allclose(angle_between(a.target, -pos), abs_angle) \ No newline at end of file + # 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 From 97d29e6c390e3c0b2c88ccb12986096e2505d1b4 Mon Sep 17 00:00:00 2001 From: James Date: Thu, 16 Nov 2023 20:19:03 +1100 Subject: [PATCH 12/14] revert unnecessary changes to framelib.IntertialFrame --- skyfield/framelib.py | 4 ---- 1 file changed, 4 deletions(-) diff --git a/skyfield/framelib.py b/skyfield/framelib.py index 4d879682d..b7d40386b 100644 --- a/skyfield/framelib.py +++ b/skyfield/framelib.py @@ -6,7 +6,6 @@ from .data.spice import inertial_frames as _inertial_frames from .functions import length_of, mxm, rot_x, rot_z - def build_matrix(): # 'xi0', 'eta0', and 'da0' are ICRS frame biases in arcseconds taken # from IERS (2003) Conventions, Chapter 5. @@ -147,9 +146,6 @@ def __init__(self, doc, matrix): self.__doc__ = doc self._matrix = matrix - def __repr__(self): - return '<{0}> {1}'.format(self.__class__.__name__, self.__doc__) - def rotation_at(self, t): return self._matrix From c76de23838eee11d2eee31c665f466b096cff63b Mon Sep 17 00:00:00 2001 From: James Date: Thu, 16 Nov 2023 20:38:35 +1100 Subject: [PATCH 13/14] add leading underscore to intersection_of --- skyfield/documentation/earth-satellites.rst | 6 +++--- skyfield/tests/test_topos.py | 14 +++++++------- skyfield/toposlib.py | 2 +- 3 files changed, 11 insertions(+), 11 deletions(-) diff --git a/skyfield/documentation/earth-satellites.rst b/skyfield/documentation/earth-satellites.rst index b5eda3d2c..1d4e703b0 100644 --- a/skyfield/documentation/earth-satellites.rst +++ b/skyfield/documentation/earth-satellites.rst @@ -550,7 +550,7 @@ 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()` +* :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 @@ -585,11 +585,11 @@ anti-aligned with the position vector: [ 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: +use the :meth:`~skyfield.toposlib.Geoid._intersection_of()` method: .. testcode:: - intersection = wgs84.intersection_of(attitude) + intersection = wgs84._intersection_of(attitude) print('Latitude:', intersection.latitude) print('Longitude:', intersection.longitude) diff --git a/skyfield/tests/test_topos.py b/skyfield/tests/test_topos.py index 58534649a..233bdfe95 100644 --- a/skyfield/tests/test_topos.py +++ b/skyfield/tests/test_topos.py @@ -270,7 +270,7 @@ def test_intersection_from_pole(ts): 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) + 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 @@ -283,7 +283,7 @@ def test_intersection_from_equator(ts): 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) + 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 @@ -314,11 +314,11 @@ def test_limb_intersection_points(ts): # 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) + 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 diff --git a/skyfield/toposlib.py b/skyfield/toposlib.py index 5a0b027b9..ce37f7c95 100644 --- a/skyfield/toposlib.py +++ b/skyfield/toposlib.py @@ -277,7 +277,7 @@ def _compute_latitude(self, position): subpoint = geographic_position_of # deprecated method name - def intersection_of(self, icrf_vector): + 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` From 4188a22475645593762766e55d9494613c45dede Mon Sep 17 00:00:00 2001 From: James Date: Sun, 3 Dec 2023 16:47:03 +1100 Subject: [PATCH 14/14] avoid using deprecated Geometric class and improve docstring --- skyfield/sgp4lib.py | 12 ++++++------ 1 file changed, 6 insertions(+), 6 deletions(-) diff --git a/skyfield/sgp4lib.py b/skyfield/sgp4lib.py index 7b0ce8e33..d500ae749 100644 --- a/skyfield/sgp4lib.py +++ b/skyfield/sgp4lib.py @@ -9,7 +9,7 @@ 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 Geometric +from .positionlib import ICRF from .searchlib import _find_discrete, find_maxima from .timelib import compute_calendar_date from .vectorlib import VectorFunction @@ -286,10 +286,10 @@ def below_horizon_at(t): 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 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 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), @@ -305,7 +305,7 @@ def _attitude(self, t, roll=0.0, pitch=0.0, yaw=0.0, rotation_order="xyz"): R = LVLH(self).rotation_at(t) direction = mxv(R, lvlh_attitude) - attitude = Geometric( + 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