Request to add attitude to positionlib.ICRF

[craigim]

I have a feature request possibly related to #492.

I'm trying to do calculations for satellite pointing where I know the attitude quaternion in ECI and would like to know the attitude in ECEF. Currently, I'm doing this by converting the quaternion to a line-of-sight vector by applying the quaternion to a unit vector along the positive Z axis, picking an arbitrary dummy distance unit, doing the coordinate transform, and then removing the units.

What would be preferable would be to define my satellite position at time T with position, velocity, and attitude in my given frame, and then have all three simultaneously transformed when I go from ecliptic_J2000_frame to itrs. I can then project my LOS vector to the Earth and figure out what exactly I'm looking at.

Failing that, is there a way to directly get the coordinate transformation quaternion between frames for a given time so I can just multiply the quaternions to get a spacecraft body local to ECEF rotation?

[brandon-rhodes]

Coordinate frames are, at the moment, functions of time that return a 3×3 rotation matrix from their .rotation_at(t) method. The matrix (if I’m remembering the direction correctly) rotates from the ITRS and to that coordinate frame.

The 3×3 rotation matrix between frames would, I suspect, be produced by multiplying the matrix of the target frame by the transpose of the origin frame matrix. Something like:

from skyfield.functions import mxm  # "matrix times matrix"
R1 = ecliptic_J2000_frame.rotation_at(t)
R2 = itrs.rotation_at(t)
R = mxm(R2, R1.T)  # Or the other way round? I always have to try it out.

Adding attitude to position and velocity is an interesting idea — maybe a subclass of the position classes could be added with that addition, to avoid making the main class more complicated? What does your code currently look like that does the transform for you with position and velocity — is it from_time_and_frame_vectors() followed by frame_xyz_and_velocity()?

[craigim]

I have a knack for always starting with the wrong rotation direction, no matter how hard I try to outwit myself 😄.

Matrices are fine by me. I prefer quaternions because that's what the spacecraft uses, but I can use scipy to convert the matrix.

Currently, I'm using astropy to do the coordinate transformations in this case.

from scipy.spatial.transform import Rotation

from astropy import coordinates as coord
from astropy import units as u
from astropy.time import Time
from astropy.utils import iers

# astrotime is the time for a time slice extracted from satellite data expressed as an astropy.time.Time object
itrs_frame = coord.ITRS(obstime=astrotime)

# eci_p and eci_v are the position and velocity Cartesian vectors for time `astrotime`
velocity = coord.CartesianDifferential(eci_v, xyz_axis=1)
position = coord.CartesianRepresentation(eci_p, xyz_axis=1, differentials=velocity)
gcrs = coord.GCRS(position, obstime=astrotime)
itrs = gcrs.transform_to(itrs_frame)

# Satellite position and velocity in ECEF (ITRS) frame
latitude = itrs.earth_location.lat.degree
longitude = itrs.earth_location.lon.degree
altitude_m = itrs.earth_location.height.value

GNC_ECEF_meters = itrs.earth_location.value.view((np.float64, 3))
GNC_ECEF_velocity = itrs.velocity.d_xyz.to(u.m / u.s).value.T

# Calculate the line-of-sight vectors in each frame
xyz_w = Rotation.from_quat(GNC_ATT_XYZW)
Z = np.array([0, 0, 1])
GNC_ECI_LOS = xyz_w.apply(Z)

attitude_astro = coord.CartesianRepresentation(GNC_ECI_LOS * u.m, xyz_axis=1)
attitude_gcrs = coord.GCRS(attitude_astro, obstime=astrotime)
attitude_itrs = attitude_gcrs.transform_to(itrs_frame)
GNC_ECEF_LOS = attitude_itrs.earth_location.value.view((np.float64, 3))

There is a lot of tacking on units, doing the calculation, and then stripping them back off, when at the end of the day it's just a rotation. When I looked at replacing the code with SkyField since I use it elsewhere and find it generally more user friendly, I saw that maybe I could abuse the code and similarly enter the values as meters, or even AU, since it's a normalized unitless vector after all and the transform is scale invariant. But I'd rather either do the rotation directly with the matrix/quaternion or else see if it would be more generally useful.

In the meantime, I think the raw rotation matrices will suit my needs