Quaternions¶

Quaternions are the standard way to represent 3D rotations in astrodynamics. Compared to Euler angles they avoid gimbal lock, and compared to rotation matrices they are more compact (4 numbers vs 9) and cheaper to compose. A unit quaternion $q = w + x\mathbf{i} + y\mathbf{j} + z\mathbf{k}$ encodes a rotation of angle $\theta$ about a unit axis $\hat{u}$ as:

$$q = \cos\frac{\theta}{2} + \sin\frac{\theta}{2}\,(u_x\mathbf{i} + u_y\mathbf{j} + u_z\mathbf{k})$$

In satkit, quaternions rotate vectors via the * operator: q * v applies the rotation encoded by q to the vector v.

Construction¶

Quaternions can be created from axis-angle, principal axis rotations, rotation matrices, or by finding the rotation between two vectors.

In [1]:

import satkit as sk
import numpy as np
import math as m

# Identity quaternion (no rotation)
q_id = sk.quaternion()
print(f"Identity: w={q_id.w}, x={q_id.x}, y={q_id.y}, z={q_id.z}")

# Rotation about principal axes
q_x = sk.quaternion.rotx(m.radians(90))   # 90 deg about x-axis
q_y = sk.quaternion.roty(m.radians(45))   # 45 deg about y-axis
q_z = sk.quaternion.rotz(m.radians(30))   # 30 deg about z-axis

# Rotation about arbitrary axis
axis = np.array([1, 1, 1]) / np.sqrt(3)  # diagonal axis
q_arb = sk.quaternion.from_axis_angle(axis, m.radians(120))
print(f"\n120 deg about [1,1,1]: angle = {m.degrees(q_arb.angle):.1f} deg")
print(f"  axis = {q_arb.axis}")

# Rotation that takes one vector to another
v1 = np.array([1.0, 0.0, 0.0])
v2 = np.array([0.0, 1.0, 0.0])
q_between = sk.quaternion.rotation_between(v1, v2)
print(f"\nRotation from x to y: angle = {m.degrees(q_between.angle):.1f} deg")
print(f"  axis = {q_between.axis}")
print(f"  check: q * [1,0,0] = {q_between * v1}")

import satkit as sk import numpy as np import math as m # Identity quaternion (no rotation) q_id = sk.quaternion() print(f"Identity: w={q_id.w}, x={q_id.x}, y={q_id.y}, z={q_id.z}") # Rotation about principal axes q_x = sk.quaternion.rotx(m.radians(90)) # 90 deg about x-axis q_y = sk.quaternion.roty(m.radians(45)) # 45 deg about y-axis q_z = sk.quaternion.rotz(m.radians(30)) # 30 deg about z-axis # Rotation about arbitrary axis axis = np.array([1, 1, 1]) / np.sqrt(3) # diagonal axis q_arb = sk.quaternion.from_axis_angle(axis, m.radians(120)) print(f"\n120 deg about [1,1,1]: angle = {m.degrees(q_arb.angle):.1f} deg") print(f" axis = {q_arb.axis}") # Rotation that takes one vector to another v1 = np.array([1.0, 0.0, 0.0]) v2 = np.array([0.0, 1.0, 0.0]) q_between = sk.quaternion.rotation_between(v1, v2) print(f"\nRotation from x to y: angle = {m.degrees(q_between.angle):.1f} deg") print(f" axis = {q_between.axis}") print(f" check: q * [1,0,0] = {q_between * v1}")

Identity: w=1.0, x=0.0, y=0.0, z=0.0

120 deg about [1,1,1]: angle = 120.0 deg
  axis = [0.57735027 0.57735027 0.57735027]

Rotation from x to y: angle = 90.0 deg
  axis = [0. 0. 1.]
  check: q * [1,0,0] = [2.22044605e-16 1.00000000e+00 0.00000000e+00]

Applying Rotations¶

The * operator rotates vectors. When two quaternions are multiplied, the rightmost rotation is applied first (same convention as rotation matrices).

In [2]:

import satkit as sk
import numpy as np
import math as m

v = np.array([1.0, 0.0, 0.0])

# Rotate x-hat by 90 deg about z: should give y-hat
q = sk.quaternion.rotz(m.radians(90))
print(f"Rz(90) * [1,0,0] = {q * v}")

# Composing rotations: first rotate 90 about z, then 90 about x
qx = sk.quaternion.rotx(m.radians(90))
qz = sk.quaternion.rotz(m.radians(90))
q_composed = qx * qz  # qz applied first, then qx
print(f"Rx(90) * Rz(90) * [1,0,0] = {q_composed * v}")

# Verify: step by step gives same result
v_step = qz * v      # [0, 1, 0]
v_step = qx * v_step  # [0, 0, 1]
print(f"Step by step:              {v_step}")

# Inverse rotation via conjugate
q = sk.quaternion.roty(m.radians(45))
v_rotated = q * v
v_back = q.conj * v_rotated
print(f"\nRound-trip: {v} -> {v_rotated} -> {v_back}")

import satkit as sk import numpy as np import math as m v = np.array([1.0, 0.0, 0.0]) # Rotate x-hat by 90 deg about z: should give y-hat q = sk.quaternion.rotz(m.radians(90)) print(f"Rz(90) * [1,0,0] = {q * v}") # Composing rotations: first rotate 90 about z, then 90 about x qx = sk.quaternion.rotx(m.radians(90)) qz = sk.quaternion.rotz(m.radians(90)) q_composed = qx * qz # qz applied first, then qx print(f"Rx(90) * Rz(90) * [1,0,0] = {q_composed * v}") # Verify: step by step gives same result v_step = qz * v # [0, 1, 0] v_step = qx * v_step # [0, 0, 1] print(f"Step by step: {v_step}") # Inverse rotation via conjugate q = sk.quaternion.roty(m.radians(45)) v_rotated = q * v v_back = q.conj * v_rotated print(f"\nRound-trip: {v} -> {v_rotated} -> {v_back}")

Rz(90) * [1,0,0] = [2.22044605e-16 1.00000000e+00 0.00000000e+00]
Rx(90) * Rz(90) * [1,0,0] = [2.22044605e-16 2.22044605e-16 1.00000000e+00]
Step by step:              [2.22044605e-16 2.22044605e-16 1.00000000e+00]

Round-trip: [1. 0. 0.] -> [ 0.70710678  0.         -0.70710678] -> [ 1.00000000e+00  0.00000000e+00 -1.38777878e-16]

Conversion To/From Other Representations¶

In [3]:

import satkit as sk
import numpy as np
import math as m

q = sk.quaternion.rotz(m.radians(30)) * sk.quaternion.roty(m.radians(20)) * sk.quaternion.rotx(m.radians(10))

# To rotation matrix (DCM)
R = q.as_rotation_matrix()
print("Rotation matrix:")
print(np.array2string(R, precision=4, suppress_small=True))

# Back from rotation matrix
q_from_R = sk.quaternion.from_rotation_matrix(R)
print(f"\nRound-trip angle error: {abs(q.angle - q_from_R.angle):.1e} rad")

# To Euler angles (roll, pitch, yaw) - intrinsic ZYX convention
roll, pitch, yaw = q.as_euler()
print(f"\nEuler angles (ZYX): roll={m.degrees(roll):.1f}, pitch={m.degrees(pitch):.1f}, yaw={m.degrees(yaw):.1f} deg")

# Axis-angle decomposition
print(f"\nAxis-angle: {m.degrees(q.angle):.2f} deg about {q.axis}")

import satkit as sk import numpy as np import math as m q = sk.quaternion.rotz(m.radians(30)) * sk.quaternion.roty(m.radians(20)) * sk.quaternion.rotx(m.radians(10)) # To rotation matrix (DCM) R = q.as_rotation_matrix() print("Rotation matrix:") print(np.array2string(R, precision=4, suppress_small=True)) # Back from rotation matrix q_from_R = sk.quaternion.from_rotation_matrix(R) print(f"\nRound-trip angle error: {abs(q.angle - q_from_R.angle):.1e} rad") # To Euler angles (roll, pitch, yaw) - intrinsic ZYX convention roll, pitch, yaw = q.as_euler() print(f"\nEuler angles (ZYX): roll={m.degrees(roll):.1f}, pitch={m.degrees(pitch):.1f}, yaw={m.degrees(yaw):.1f} deg") # Axis-angle decomposition print(f"\nAxis-angle: {m.degrees(q.angle):.2f} deg about {q.axis}")

Rotation matrix:
[[ 0.8138 -0.441   0.3785]
 [ 0.4698  0.8826  0.018 ]
 [-0.342   0.1632  0.9254]]

Round-trip angle error: 1.1e-16 rad

Euler angles (ZYX): roll=10.0, pitch=20.0, yaw=30.0 deg

Axis-angle: 35.82 deg about [0.12401544 0.61563806 0.77820945]

SLERP: Smooth Interpolation¶

Spherical Linear Interpolation (SLERP) produces smooth, constant-angular-velocity rotation between two orientations. This is essential for attitude interpolation — linear interpolation of Euler angles produces uneven motion and can fail near gimbal lock.

In [4]:

import satkit as sk
import numpy as np
import math as m
import matplotlib.pyplot as plt
import scienceplots  # noqa: F401
plt.style.use(["science", "no-latex", "../satkit.mplstyle"])
%config InlineBackend.figure_formats = ['svg']

# Interpolate between two attitudes
q_start = sk.quaternion()  # identity
q_end = sk.quaternion.rotz(m.radians(120)) * sk.quaternion.rotx(m.radians(60))

fracs = np.linspace(0, 1, 50)
q_interp = [q_start.slerp(q_end, f) for f in fracs]

# Track where the x-axis points during the rotation
x_hat = np.array([1.0, 0.0, 0.0])
trajectories = np.array([q * x_hat for q in q_interp])

# Also track the rotation angle over time — should be linear for SLERP
angles = [q.angle for q in q_interp]

fig = plt.figure(figsize=(12, 5))

ax = fig.add_subplot(121, projection='3d')
ax.plot(trajectories[:, 0], trajectories[:, 1], trajectories[:, 2], 'o-', markersize=3)
ax.scatter(*trajectories[0], color='C2', s=80, zorder=5, label='Start')
ax.scatter(*trajectories[-1], color='C3', s=80, zorder=5, label='End')
# Draw unit sphere wireframe
u, v = np.mgrid[0:2*np.pi:20j, 0:np.pi:10j]
ax.plot_wireframe(np.cos(u)*np.sin(v), np.sin(u)*np.sin(v), np.cos(v),
                  color='gray', alpha=0.1, linewidth=0.5)
ax.set_xlabel('X')
ax.set_ylabel('Y')
ax.set_zlabel('Z')
ax.set_title('SLERP Path of $\\hat{x}$ on Unit Sphere')
ax.legend()

ax2 = fig.add_subplot(122)
ax2.plot(fracs, np.degrees(angles))
ax2.set_xlabel('Interpolation Fraction')
ax2.set_ylabel('Rotation Angle [deg]')
ax2.set_title('SLERP Produces Constant Angular Rate')

plt.tight_layout()
plt.show()

import satkit as sk import numpy as np import math as m import matplotlib.pyplot as plt import scienceplots # noqa: F401 plt.style.use(["science", "no-latex", "../satkit.mplstyle"]) %config InlineBackend.figure_formats = ['svg'] # Interpolate between two attitudes q_start = sk.quaternion() # identity q_end = sk.quaternion.rotz(m.radians(120)) * sk.quaternion.rotx(m.radians(60)) fracs = np.linspace(0, 1, 50) q_interp = [q_start.slerp(q_end, f) for f in fracs] # Track where the x-axis points during the rotation x_hat = np.array([1.0, 0.0, 0.0]) trajectories = np.array([q * x_hat for q in q_interp]) # Also track the rotation angle over time — should be linear for SLERP angles = [q.angle for q in q_interp] fig = plt.figure(figsize=(12, 5)) ax = fig.add_subplot(121, projection='3d') ax.plot(trajectories[:, 0], trajectories[:, 1], trajectories[:, 2], 'o-', markersize=3) ax.scatter(*trajectories[0], color='C2', s=80, zorder=5, label='Start') ax.scatter(*trajectories[-1], color='C3', s=80, zorder=5, label='End') # Draw unit sphere wireframe u, v = np.mgrid[0:2*np.pi:20j, 0:np.pi:10j] ax.plot_wireframe(np.cos(u)*np.sin(v), np.sin(u)*np.sin(v), np.cos(v), color='gray', alpha=0.1, linewidth=0.5) ax.set_xlabel('X') ax.set_ylabel('Y') ax.set_zlabel('Z') ax.set_title('SLERP Path of $\\hat{x}$ on Unit Sphere') ax.legend() ax2 = fig.add_subplot(122) ax2.plot(fracs, np.degrees(angles)) ax2.set_xlabel('Interpolation Fraction') ax2.set_ylabel('Rotation Angle [deg]') ax2.set_title('SLERP Produces Constant Angular Rate') plt.tight_layout() plt.show()

Quaternions for Frame Rotations¶

In satellite astrodynamics, the most common use of quaternions is rotating between coordinate frames. satkit.frametransform provides quaternions for all standard frame pairs.

In [5]:

import satkit as sk
import numpy as np

t = sk.time(2024, 6, 15, 12, 0, 0)

# A ground station in Boston
boston = sk.itrfcoord(latitude_deg=42.36, longitude_deg=-71.06, altitude=0)

# Position in ITRF (Earth-fixed, meters)
pos_itrf = boston.vector
print(f"Boston ITRF: [{pos_itrf[0]/1e3:.1f}, {pos_itrf[1]/1e3:.1f}, {pos_itrf[2]/1e3:.1f}] km")

# Rotate to GCRF (inertial) — the Earth has rotated, so the
# inertial position depends on the time
q = sk.frametransform.qitrf2gcrf(t)
pos_gcrf = q * pos_itrf
print(f"Boston GCRF: [{pos_gcrf[0]/1e3:.1f}, {pos_gcrf[1]/1e3:.1f}, {pos_gcrf[2]/1e3:.1f}] km")

# The magnitude is preserved (rigid-body rotation)
print(f"\n|ITRF| = {np.linalg.norm(pos_itrf)/1e3:.3f} km")
print(f"|GCRF| = {np.linalg.norm(pos_gcrf)/1e3:.3f} km")

# Round-trip: ITRF -> GCRF -> ITRF
pos_back = q.conj * pos_gcrf
print(f"\nRound-trip error: {np.linalg.norm(pos_back - pos_itrf):.2e} m")

# Local NED frame at Boston
q_ned = boston.qned2itrf
north_itrf = q_ned * np.array([1.0, 0.0, 0.0])  # north in ITRF
print(f"\nNorth direction at Boston (ITRF): {north_itrf}")

import satkit as sk import numpy as np t = sk.time(2024, 6, 15, 12, 0, 0) # A ground station in Boston boston = sk.itrfcoord(latitude_deg=42.36, longitude_deg=-71.06, altitude=0) # Position in ITRF (Earth-fixed, meters) pos_itrf = boston.vector print(f"Boston ITRF: [{pos_itrf[0]/1e3:.1f}, {pos_itrf[1]/1e3:.1f}, {pos_itrf[2]/1e3:.1f}] km") # Rotate to GCRF (inertial) — the Earth has rotated, so the # inertial position depends on the time q = sk.frametransform.qitrf2gcrf(t) pos_gcrf = q * pos_itrf print(f"Boston GCRF: [{pos_gcrf[0]/1e3:.1f}, {pos_gcrf[1]/1e3:.1f}, {pos_gcrf[2]/1e3:.1f}] km") # The magnitude is preserved (rigid-body rotation) print(f"\n|ITRF| = {np.linalg.norm(pos_itrf)/1e3:.3f} km") print(f"|GCRF| = {np.linalg.norm(pos_gcrf)/1e3:.3f} km") # Round-trip: ITRF -> GCRF -> ITRF pos_back = q.conj * pos_gcrf print(f"\nRound-trip error: {np.linalg.norm(pos_back - pos_itrf):.2e} m") # Local NED frame at Boston q_ned = boston.qned2itrf north_itrf = q_ned * np.array([1.0, 0.0, 0.0]) # north in ITRF print(f"\nNorth direction at Boston (ITRF): {north_itrf}")

Boston ITRF: [1532.1, -4464.6, 4275.2] km
Boston GCRF: [4611.3, 1053.1, 4264.3] km

|ITRF| = 6368.473 km
|GCRF| = 6368.473 km

Round-trip error: 3.49e-09 m

North direction at Boston (ITRF): [-0.21869622  0.63730719  0.73892591]