{
"cells": [
{
"cell_type": "markdown",
"id": "cell-0",
"metadata": {},
"source": [
"# Quaternions\n",
"\n",
"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:\n",
"\n",
"$$q = \\cos\\frac{\\theta}{2} + \\sin\\frac{\\theta}{2}\\,(u_x\\mathbf{i} + u_y\\mathbf{j} + u_z\\mathbf{k})$$\n",
"\n",
"In `satkit`, quaternions rotate vectors via the `*` operator: `q * v` applies the rotation encoded by `q` to the vector `v`."
]
},
{
"cell_type": "markdown",
"id": "cell-1",
"metadata": {},
"source": [
"## Construction\n",
"\n",
"Quaternions can be created from axis-angle, principal axis rotations, rotation matrices, or by finding the rotation between two vectors."
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "cell-2",
"metadata": {},
"outputs": [],
"source": [
"import satkit as sk\n",
"import numpy as np\n",
"import math as m\n",
"\n",
"# Identity quaternion (no rotation)\n",
"q_id = sk.quaternion()\n",
"print(f\"Identity: w={q_id.w}, x={q_id.x}, y={q_id.y}, z={q_id.z}\")\n",
"\n",
"# Rotation about principal axes\n",
"q_x = sk.quaternion.rotx(m.radians(90)) # 90 deg about x-axis\n",
"q_y = sk.quaternion.roty(m.radians(45)) # 45 deg about y-axis\n",
"q_z = sk.quaternion.rotz(m.radians(30)) # 30 deg about z-axis\n",
"\n",
"# Rotation about arbitrary axis\n",
"axis = np.array([1, 1, 1]) / np.sqrt(3) # diagonal axis\n",
"q_arb = sk.quaternion.from_axis_angle(axis, m.radians(120))\n",
"print(f\"\\n120 deg about [1,1,1]: angle = {m.degrees(q_arb.angle):.1f} deg\")\n",
"print(f\" axis = {q_arb.axis}\")\n",
"\n",
"# Rotation that takes one vector to another\n",
"v1 = np.array([1.0, 0.0, 0.0])\n",
"v2 = np.array([0.0, 1.0, 0.0])\n",
"q_between = sk.quaternion.rotation_between(v1, v2)\n",
"print(f\"\\nRotation from x to y: angle = {m.degrees(q_between.angle):.1f} deg\")\n",
"print(f\" axis = {q_between.axis}\")\n",
"print(f\" check: q * [1,0,0] = {q_between * v1}\")"
]
},
{
"cell_type": "markdown",
"id": "cell-3",
"metadata": {},
"source": [
"## Applying Rotations\n",
"\n",
"The `*` operator rotates vectors. When two quaternions are multiplied, the rightmost rotation is applied first (same convention as rotation matrices)."
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "cell-4",
"metadata": {},
"outputs": [],
"source": [
"import satkit as sk\n",
"import numpy as np\n",
"import math as m\n",
"\n",
"v = np.array([1.0, 0.0, 0.0])\n",
"\n",
"# Rotate x-hat by 90 deg about z: should give y-hat\n",
"q = sk.quaternion.rotz(m.radians(90))\n",
"print(f\"Rz(90) * [1,0,0] = {q * v}\")\n",
"\n",
"# Composing rotations: first rotate 90 about z, then 90 about x\n",
"qx = sk.quaternion.rotx(m.radians(90))\n",
"qz = sk.quaternion.rotz(m.radians(90))\n",
"q_composed = qx * qz # qz applied first, then qx\n",
"print(f\"Rx(90) * Rz(90) * [1,0,0] = {q_composed * v}\")\n",
"\n",
"# Verify: step by step gives same result\n",
"v_step = qz * v # [0, 1, 0]\n",
"v_step = qx * v_step # [0, 0, 1]\n",
"print(f\"Step by step: {v_step}\")\n",
"\n",
"# Inverse rotation via conjugate\n",
"q = sk.quaternion.roty(m.radians(45))\n",
"v_rotated = q * v\n",
"v_back = q.conj * v_rotated\n",
"print(f\"\\nRound-trip: {v} -> {v_rotated} -> {v_back}\")"
]
},
{
"cell_type": "markdown",
"id": "cell-5",
"metadata": {},
"source": [
"## Conversion To/From Other Representations"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "cell-6",
"metadata": {},
"outputs": [],
"source": [
"import satkit as sk\n",
"import numpy as np\n",
"import math as m\n",
"\n",
"q = sk.quaternion.rotz(m.radians(30)) * sk.quaternion.roty(m.radians(20)) * sk.quaternion.rotx(m.radians(10))\n",
"\n",
"# To rotation matrix (DCM)\n",
"R = q.as_rotation_matrix()\n",
"print(\"Rotation matrix:\")\n",
"print(np.array2string(R, precision=4, suppress_small=True))\n",
"\n",
"# Back from rotation matrix\n",
"q_from_R = sk.quaternion.from_rotation_matrix(R)\n",
"print(f\"\\nRound-trip angle error: {abs(q.angle - q_from_R.angle):.1e} rad\")\n",
"\n",
"# To Euler angles (roll, pitch, yaw) - intrinsic ZYX convention\n",
"roll, pitch, yaw = q.as_euler()\n",
"print(f\"\\nEuler angles (ZYX): roll={m.degrees(roll):.1f}, pitch={m.degrees(pitch):.1f}, yaw={m.degrees(yaw):.1f} deg\")\n",
"\n",
"# Axis-angle decomposition\n",
"print(f\"\\nAxis-angle: {m.degrees(q.angle):.2f} deg about {q.axis}\")"
]
},
{
"cell_type": "markdown",
"id": "cell-7",
"metadata": {},
"source": [
"## SLERP: Smooth Interpolation\n",
"\n",
"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."
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "cell-8",
"metadata": {},
"outputs": [],
"source": [
"import satkit as sk\n",
"import numpy as np\n",
"import math as m\n",
"import matplotlib.pyplot as plt\n",
"import scienceplots # noqa: F401\n",
"plt.style.use([\"science\", \"no-latex\", \"../satkit.mplstyle\"])\n",
"%config InlineBackend.figure_formats = ['svg']\n",
"\n",
"# Interpolate between two attitudes\n",
"q_start = sk.quaternion() # identity\n",
"q_end = sk.quaternion.rotz(m.radians(120)) * sk.quaternion.rotx(m.radians(60))\n",
"\n",
"fracs = np.linspace(0, 1, 50)\n",
"q_interp = [q_start.slerp(q_end, f) for f in fracs]\n",
"\n",
"# Track where the x-axis points during the rotation\n",
"x_hat = np.array([1.0, 0.0, 0.0])\n",
"trajectories = np.array([q * x_hat for q in q_interp])\n",
"\n",
"# Also track the rotation angle over time — should be linear for SLERP\n",
"angles = [q.angle for q in q_interp]\n",
"\n",
"fig = plt.figure(figsize=(12, 5))\n",
"\n",
"ax = fig.add_subplot(121, projection='3d')\n",
"ax.plot(trajectories[:, 0], trajectories[:, 1], trajectories[:, 2], 'o-', markersize=3)\n",
"ax.scatter(*trajectories[0], color='C2', s=80, zorder=5, label='Start')\n",
"ax.scatter(*trajectories[-1], color='C3', s=80, zorder=5, label='End')\n",
"# Draw unit sphere wireframe\n",
"u, v = np.mgrid[0:2*np.pi:20j, 0:np.pi:10j]\n",
"ax.plot_wireframe(np.cos(u)*np.sin(v), np.sin(u)*np.sin(v), np.cos(v),\n",
" color='gray', alpha=0.1, linewidth=0.5)\n",
"ax.set_xlabel('X')\n",
"ax.set_ylabel('Y')\n",
"ax.set_zlabel('Z')\n",
"ax.set_title('SLERP Path of $\\\\hat{x}$ on Unit Sphere')\n",
"ax.legend()\n",
"\n",
"ax2 = fig.add_subplot(122)\n",
"ax2.plot(fracs, np.degrees(angles))\n",
"ax2.set_xlabel('Interpolation Fraction')\n",
"ax2.set_ylabel('Rotation Angle [deg]')\n",
"ax2.set_title('SLERP Produces Constant Angular Rate')\n",
"\n",
"plt.tight_layout()\n",
"plt.show()"
]
},
{
"cell_type": "markdown",
"id": "cell-9",
"metadata": {},
"source": [
"## Quaternions for Frame Rotations\n",
"\n",
"In satellite astrodynamics, the most common use of quaternions is rotating between coordinate frames. `satkit.frametransform` provides quaternions for all standard frame pairs."
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "cell-10",
"metadata": {},
"outputs": [],
"source": [
"import satkit as sk\n",
"import numpy as np\n",
"\n",
"t = sk.time(2024, 6, 15, 12, 0, 0)\n",
"\n",
"# A ground station in Boston\n",
"boston = sk.itrfcoord(latitude_deg=42.36, longitude_deg=-71.06, altitude=0)\n",
"\n",
"# Position in ITRF (Earth-fixed, meters)\n",
"pos_itrf = boston.vector\n",
"print(f\"Boston ITRF: [{pos_itrf[0]/1e3:.1f}, {pos_itrf[1]/1e3:.1f}, {pos_itrf[2]/1e3:.1f}] km\")\n",
"\n",
"# Rotate to GCRF (inertial) — the Earth has rotated, so the\n",
"# inertial position depends on the time\n",
"q = sk.frametransform.qitrf2gcrf(t)\n",
"pos_gcrf = q * pos_itrf\n",
"print(f\"Boston GCRF: [{pos_gcrf[0]/1e3:.1f}, {pos_gcrf[1]/1e3:.1f}, {pos_gcrf[2]/1e3:.1f}] km\")\n",
"\n",
"# The magnitude is preserved (rigid-body rotation)\n",
"print(f\"\\n|ITRF| = {np.linalg.norm(pos_itrf)/1e3:.3f} km\")\n",
"print(f\"|GCRF| = {np.linalg.norm(pos_gcrf)/1e3:.3f} km\")\n",
"\n",
"# Round-trip: ITRF -> GCRF -> ITRF\n",
"pos_back = q.conj * pos_gcrf\n",
"print(f\"\\nRound-trip error: {np.linalg.norm(pos_back - pos_itrf):.2e} m\")\n",
"\n",
"# Local NED frame at Boston\n",
"q_ned = boston.qned2itrf\n",
"north_itrf = q_ned * np.array([1.0, 0.0, 0.0]) # north in ITRF\n",
"print(f\"\\nNorth direction at Boston (ITRF): {north_itrf}\")"
]
}
],
"metadata": {
"kernelspec": {
"display_name": "Python 3 (ipykernel)",
"language": "python",
"name": "python3"
},
"language_info": {
"name": "python",
"version": "3.12.0"
}
},
"nbformat": 4,
"nbformat_minor": 5
}