#!/usr/bin/env python3 """ elements-to-rv.py ================= The forward map as three coordinate transforms: (a,e,i,Omega,omega,f) -> (r,v). Part IV is two steps. The first is trivial -- the conic in its own plane, r_PQW = r (cos f, sin f, 0), v_PQW = sqrt(mu/p) (-sin f, e + cos f, 0), built with the perifocal frame {p, q, w} still sitting on top of the inertial frame {I, J, K}. The second step carries that frame to its true attitude, and it is the one worth animating, because it is three rotations and each one is about a DIFFERENT axis: R3(Omega) about K -- the 1-axis lands on the line of nodes R1(i) about the frame's own 1-axis (the node line) -- the 3-axis lands on h R3(omega) about the frame's own 3-axis (h) -- the 1-axis lands on e Q = R3(Omega) R1(i) R3(omega), r_IJK = Q r_PQW, v_IJK = Q v_PQW. Read left to right, those are INTRINSIC rotations: each one turns the frame about an axis the *previous* rotation put where it is, which is why the composition multiplies on the right. That is the reading animated here, and it is the same matrix as the deck's right-to-left reading, where the rotations are extrinsic -- about the fixed inertial axes -- and undo the sequence instead of building it. Same Q either way; only the story differs. The frame keeps its name throughout: it is always {p, q, w}. What changes is where it points, and at each stage one of its axes acquires a geometric meaning, which is what the axis labels record. Usage ----- python codes/elements-to-rv.py # step through by hand python codes/elements-to-rv.py --verify # the matrices, in text python codes/elements-to-rv.py --png img/fig-elements-to-rv.png python codes/elements-to-rv.py --storyboard img/fig-elements-to-rv-steps.png python codes/elements-to-rv.py --save img/elements-to-rv.gif python codes/elements-to-rv.py --example 5 # a hyperbolic flyby python codes/elements-to-rv.py --elements 12000 0.4 60 300 250 250 Units are km and km/s throughout, mu = 398600 km^3/s^2 (Earth). """ from __future__ import annotations import argparse import os import numpy as np MU_EARTH = 398600.0 # -- palette (shared with the other figures in this deck) --------------------- CREAM, INK, DIM = "#fcfbf7", "#2f2a24", "#8b857c" MAROON, BLUE, GREEN, AMBER, PURPLE = ("#6b1f1f", "#2f5ea8", "#2d6a2d", "#b45309", "#7e22ce") TEAL = "#0f766e" # velocity, and the 2-axis, which never earns a name AXIS = "#6b655c" # the inertial frame is backdrop, not an element # Element sets worked in Part V of the deck. Example 4 is the default so that # this animation and rv-to-elements.py are literally inverse to each other, on # the same orbit and from the same viewpoint. EXAMPLES = { "1": dict(a=8788.0951, e=0.171212, inc=153.2492, Om=255.2793, om=20.0683, f=28.4456, title="Examples 1-2 -- the round trip of the retrograde orbit"), "4": dict(a=12000.0, e=0.4, inc=60.0, Om=300.0, om=250.0, f=250.0, title="Example 4 -- the inverse of the " r"$(\mathbf{r},\mathbf{v})\to$ elements animation"), "5": dict(a=-16725.1964, e=1.4, inc=30.0, Om=40.0, om=60.0, f=30.0, title="Example 5 -- a hyperbolic flyby; the forward map never " "asks what conic it is drawing"), } # -- the three rotations ------------------------------------------------------ def R1(t): """Active rotation of a vector about the 1-axis by t (radians).""" c, s = np.cos(t), np.sin(t) return np.array([[1.0, 0.0, 0.0], [0.0, c, -s], [0.0, s, c]]) def R3(t): """Active rotation of a vector about the 3-axis by t (radians).""" c, s = np.cos(t), np.sin(t) return np.array([[c, -s, 0.0], [s, c, 0.0], [0.0, 0.0, 1.0]]) def attitude(Om, inc, om, stage_fracs): """The frame's attitude part way through the sequence. `stage_fracs` is (fO, fi, fw) in [0, 1]: how much of each of the three rotations has been applied so far. Because the product grows on the RIGHT, partial rotation k turns the frame about an axis that rotations 1..k-1 have already moved -- which is exactly what makes the sequence intrinsic. """ fO, fi, fw = stage_fracs return R3(fO * Om) @ R1(fi * inc) @ R3(fw * om) def rotation_axis(Om, inc, which): """The axis of each rotation, in inertial components. Rotation 1 turns about K. Rotation 2 turns about the frame's 1-axis as the first rotation left it, R3(Omega) e1 = n. Rotation 3 turns about the frame's 3-axis as the first two left it, R3(Omega) R1(i) e3 = h. """ if which == "Om": return np.array([0.0, 0.0, 1.0]) if which == "inc": return R3(Om) @ np.array([1.0, 0.0, 0.0]) return R3(Om) @ R1(inc) @ np.array([0.0, 0.0, 1.0]) def perifocal_state(a, e, f, mu=MU_EARTH): """Step 1: the conic in its own plane (Battin Eq. 3.36).""" p = a * (1.0 - e**2) r = p / (1.0 + e * np.cos(f)) r_pqw = r * np.array([np.cos(f), np.sin(f), 0.0]) v_pqw = np.sqrt(mu / p) * np.array([-np.sin(f), e + np.cos(f), 0.0]) return p, r, r_pqw, v_pqw def compose(prm, mu=MU_EARTH): """Everything the picture and the panel need, in one dictionary.""" a, e = prm["a"], prm["e"] inc, Om, om, f = (np.deg2rad(prm[k]) for k in ("inc", "Om", "om", "f")) p, r, r_pqw, v_pqw = perifocal_state(a, e, f, mu) Q = R3(Om) @ R1(inc) @ R3(om) return dict(mu=mu, a=a, e=e, inc=inc, Om=Om, om=om, f=f, p=p, r=r, r_pqw=r_pqw, v_pqw=v_pqw, Q=Q, r_ijk=Q @ r_pqw, v_ijk=Q @ v_pqw, h=np.sqrt(mu * p), r_p=p / (1.0 + e), r_a=(p / (1.0 - e) if e < 1.0 else np.inf), T=(2 * np.pi * np.sqrt(a**3 / mu) if a > 0 else np.inf)) def elements_from_state(r_vec, v_vec, mu=MU_EARTH): """The inverse map of Part III -- used only to check the round trip.""" r = np.linalg.norm(r_vec) h_vec = np.cross(r_vec, v_vec) n_vec = np.cross([0.0, 0.0, 1.0], h_vec) n = np.linalg.norm(n_vec) e_vec = np.cross(v_vec, h_vec) / mu - r_vec / r e = np.linalg.norm(e_vec) inc = np.arccos(np.clip(h_vec[2] / np.linalg.norm(h_vec), -1.0, 1.0)) Om = np.arctan2(n_vec[1], n_vec[0]) % (2 * np.pi) if n > 0 else 0.0 om = np.arccos(np.clip((n_vec @ e_vec) / (n * e), -1.0, 1.0)) if n > 0 else 0.0 if e_vec[2] < 0.0: om = 2 * np.pi - om f = np.arccos(np.clip((e_vec @ r_vec) / (e * r), -1.0, 1.0)) if r_vec @ v_vec < 0.0: f = 2 * np.pi - f a = 1.0 / (2.0 / r - (v_vec @ v_vec) / mu) return dict(a=a, e=e, inc=inc, Om=Om, om=om, f=f) # -- geometry helpers --------------------------------------------------------- def rodrigues(axis, ang): """Rotation matrix about a unit `axis` by `ang` (active rotation).""" k = np.asarray(axis, dtype=float) k = k / np.linalg.norm(k) K = np.array([[0.0, -k[2], k[1]], [k[2], 0.0, -k[0]], [-k[1], k[0], 0.0]]) return np.eye(3) + np.sin(ang) * K + (1.0 - np.cos(ang)) * (K @ K) def _basis(nhat, prefer=None): """An orthonormal pair spanning the plane through the origin normal nhat.""" nhat = np.asarray(nhat, dtype=float) u = np.array([1.0, 0.0, 0.0]) if prefer is None else np.asarray(prefer, float) if abs(u @ nhat) > 0.9: u = np.array([0.0, 1.0, 0.0]) u = u - (u @ nhat) * nhat u /= np.linalg.norm(u) return u, np.cross(nhat, u) def _square(nhat, half, prefer=None): """Corners of a square patch of the plane through the origin, normal nhat.""" u, w = _basis(nhat, prefer) return np.array([half * (su * u + sw * w) for su, sw in ((1, 1), (-1, 1), (-1, -1), (1, -1))]) def _grid(nhat, half, n=4, prefer=None): """Line segments of a square grid ruled on that plane (for depth cues).""" u, w = _basis(nhat, prefer) segs = [] for t in np.linspace(-half, half, n + 1): segs.append(np.array([t * u - half * w, t * u + half * w])) segs.append(np.array([-half * u + t * w, half * u + t * w])) return segs def _arc(v_from, axis, ang, radius, n=90): """Points of an arc: `v_from` swept about `axis` through `ang`, scaled.""" v = np.asarray(v_from, dtype=float) v = v / np.linalg.norm(v) ts = np.linspace(0.0, ang, max(int(n), 2)) return np.array([radius * (rodrigues(axis, t) @ v) for t in ts]) def _sector(v_from, axis, ang, radius, n=90): """The arc closed back on the vertex: the piece of plane the angle lives in.""" return np.vstack([np.zeros(3), _arc(v_from, axis, ang, radius, n=n)]) def ease(t): """Cosine ease in/out on [0, 1] -- no jerk at the start or end of a step.""" return 0.5 * (1.0 - np.cos(np.pi * np.clip(t, 0.0, 1.0))) # -- the steps ---------------------------------------------------------------- STEPS = [ ("given", "Given: the six elements. No inverse trigonometry anywhere " "in what follows."), ("pqw", r"**Step 1.** The conic in its own plane, with " r"$\{\hat{\mathbf{p}},\hat{\mathbf{q}},\hat{\mathbf{w}}\}$ " r"still lying on $\{\hat{\mathbf{I}},\hat{\mathbf{J}}," r"\hat{\mathbf{K}}\}$"), ("Om", r"**Step 2a.** $\mathbf{R}_3(\Omega)$, about $\hat{\mathbf{K}}$ " r"$\rightarrow$ the 1-axis lands on the line of nodes"), ("inc", r"**Step 2b.** $\mathbf{R}_1(i)$, about the frame's own 1-axis " r"$\rightarrow$ the 3-axis lands on $\mathbf{h}$"), ("om", r"**Step 2c.** $\mathbf{R}_3(\omega)$, about the frame's own " r"3-axis, $\mathbf{h}$ $\rightarrow$ the 1-axis lands on " r"$\mathbf{e}$"), ("state", r"**Done.** $\mathbf{r}=\mathbf{Q}\,\mathbf{r}_{PQW}$, " r"$\mathbf{v}=\mathbf{Q}\,\mathbf{v}_{PQW}$ -- a rotation, so " r"no length changed"), ] STEP_KEYS = [k for k, _ in STEPS] def _reveal(key, step, t): """How far along the reveal of `key` we are: 0 before, t during, 1 after.""" j = STEP_KEYS.index(key) if j < step: return 1.0 if j > step: return 0.0 return float(np.clip(t, 0.0, 1.0)) def _fracs(step, t): """How much of each of the three rotations has been applied.""" return (_reveal("Om", step, t), _reveal("inc", step, t), _reveal("om", step, t)) def triad_style(step): """(colour, label) for the 1-, 2- and 3-axis, as they acquire meanings. The frame is always {p, q, w}; the aliases are what the rotation just finished has made true of it. """ j = step one = ((INK, r"$\hat{\mathbf{p}}$") if j < 2 else (GREEN, r"$\hat{\mathbf{p}}=\hat{\mathbf{n}}$") if j < 4 else (PURPLE, r"$\hat{\mathbf{p}}=\hat{\mathbf{e}}$")) two = (TEAL, r"$\hat{\mathbf{q}}$") three = ((INK, r"$\hat{\mathbf{w}}$") if j < 3 else (BLUE, r"$\hat{\mathbf{w}}=\hat{\mathbf{h}}$")) return one, two, three def draw_scene(ax, S, step, t=1.0, zoom=1.68): """Redraw the 3D scene as it stands `t` of the way through `step`.""" from mpl_toolkits.mplot3d.art3d import Poly3DCollection elev, azim = ax.elev, ax.azim ax.clear() e, inc, Om, om, f = S["e"], S["inc"], S["Om"], S["om"], S["f"] zhat = np.array([0.0, 0.0, 1.0]) M = attitude(Om, inc, om, _fracs(step, t)) # the frame, right now p1, p2, p3 = (M @ np.eye(3)).T # its three axes L = (S["r_a"] if np.isfinite(S["r_a"]) else 2.4 * S["r_p"]) * 1.12 L = max(L, 1.25 * S["r"]) v_scale = 0.42 * L / np.linalg.norm(S["v_pqw"]) # -- reference plane, its grid, and the inertial axes --------------------- sq = L * np.array([[1, 1, 0], [-1, 1, 0], [-1, -1, 0], [1, -1, 0]], float) ax.add_collection3d(Poly3DCollection([sq], facecolor=DIM, alpha=0.09, edgecolor=DIM, linewidths=0.8)) for seg in _grid(zhat, L, n=6): ax.plot(*seg.T, color=DIM, lw=0.45, alpha=0.35) ax.text(-0.50 * L, -0.95 * L, 0.0, "reference plane (equator)", color=DIM, fontsize=9.5, ha="center") for vec, lab, ln, ratio in (([1.0, 0, 0], r"$\hat{\mathbf{I}}$", 1.02, 0.10), ([0, 1.0, 0], r"$\hat{\mathbf{J}}$", 0.95, 0.09), ([0, 0, 1.0], r"$\hat{\mathbf{K}}$", 0.80, 0.11)): ax.quiver(0, 0, 0, *(ln * L * np.asarray(vec)), color=AXIS, lw=1.2, arrow_length_ratio=ratio) ax.text(*((ln + 0.08) * L * np.asarray(vec)), lab, color=AXIS, fontsize=11, ha="center") ax.plot([0], [0], [0], "o", color=AMBER, ms=11, zorder=8) # -- the axis this step turns about, lit up before anything else --------- live = STEP_KEYS[step] if step < len(STEPS) else None if live in ("Om", "inc", "om"): ax_vec = rotation_axis(Om, inc, live) col = {"Om": AMBER, "inc": GREEN, "om": BLUE}[live] ax.plot(*np.vstack([-0.45 * L * ax_vec, 1.05 * L * ax_vec]).T, color=col, lw=7.0, alpha=0.20, solid_capstyle="round") # Only the first rotation gets the words. The other two turn about an # axis that IS one of the triad's own axes -- p for R1(i), w for # R3(omega) -- and that arrow already carries its label in this very # colour, so a second label would only collide with the first. if live == "Om": ax.text(*(1.02 * L * ax_vec + 0.16 * L * np.array([1.0, 0.0, 0.0])), "rotation axis", color=col, fontsize=9.5, ha="center", va="bottom", alpha=0.95) # -- the perifocal triad, in its current attitude ------------------------- g_frame = _reveal("pqw", step, t) if g_frame > 0: # each label is nudged along the NEXT axis of the triad, which both # separates the three and clears the inertial labels the triad sits # exactly on top of at step 1 nudge = (p2, p3, -p1) for vec, off, (col, lab) in zip((p1, p2, p3), nudge, triad_style(step)): ax.quiver(0, 0, 0, *(g_frame * 0.72 * L * vec), color=col, lw=2.6, arrow_length_ratio=0.12 / max(g_frame, 0.3)) ax.text(*(g_frame * (0.80 * L * vec + 0.10 * L * off)), lab, color=col, fontsize=13, ha="center", va="center") # -- the conic, carried along by the frame -------------------------------- if g_frame > 0: f_max = np.pi if e < 1.0 else 0.985 * np.arccos(-1.0 / e) fs = np.linspace(-f_max, f_max, 721) rs = S["p"] / (1.0 + e * np.cos(fs)) orb = (rs[:, None] * np.stack([np.cos(fs), np.sin(fs), np.zeros_like(fs)], 1)) @ M.T orb = orb[np.linalg.norm(orb, axis=1) <= 1.45 * L] if len(orb) > 1: ax.plot(*orb.T, color=INK, lw=2.4, alpha=g_frame) shadow = orb.copy() shadow[:, 2] = 0.0 ax.plot(*shadow.T, color=INK, lw=0.9, ls=(0, (3, 4)), alpha=0.42 * g_frame) # the orbital plane, once the tilt has put it somewhere of its own if _reveal("inc", step, t) > 0: gi = _reveal("inc", step, t) quad = _square(p3, 0.86 * L, prefer=p1) ax.add_collection3d(Poly3DCollection([quad], facecolor=BLUE, alpha=0.11 * gi, edgecolor=BLUE, linewidths=0.9)) # -- the state, carried along with it ------------------------------------- if g_frame > 0: r_vec = M @ S["r_pqw"] v_vec = M @ S["v_pqw"] ax.quiver(0, 0, 0, *(g_frame * r_vec), color=MAROON, lw=2.4, arrow_length_ratio=0.10 / max(g_frame, 0.3)) ax.plot(*(g_frame * r_vec).reshape(3, 1), "o", color=MAROON, ms=9, zorder=9) q_in = np.cross(p3, r_vec / np.linalg.norm(r_vec)) # in-plane normal ax.text(*(0.80 * g_frame * r_vec + 0.10 * L * q_in), r"$\mathbf{r}$", color=MAROON, fontsize=14, ha="center", va="center") ax.quiver(*(g_frame * r_vec), *(g_frame * v_scale * v_vec), color=TEAL, lw=2.4, arrow_length_ratio=0.12 / max(g_frame, 0.3)) ax.text(*(g_frame * (r_vec + v_scale * v_vec) + 0.07 * L * (v_vec / np.linalg.norm(v_vec))), r"$\mathbf{v}$", color=TEAL, fontsize=14, ha="center", va="center") # the true anomaly rides with the frame: it is a perifocal angle _angle(ax, p1, p3, f, 0.12 * L, MAROON, r"$f$", alpha=0.9 * g_frame, lab_r=1.34, lab_off=-0.55) # periapsis, the point the 1-axis is aiming for ax.plot(*(g_frame * S["r_p"] * p1).reshape(3, 1), "o", color=PURPLE, ms=7, zorder=9) # -- each angle stays once its rotation has drawn it ---------------------- nhat = rotation_axis(Om, inc, "inc") hhat = rotation_axis(Om, inc, "om") g = _reveal("Om", step, t) if g > 0: _angle(ax, np.array([1.0, 0.0, 0.0]), zhat, g * Om, 0.52 * L, AMBER, r"$\Omega$", alpha=1.0 if live == "Om" else 0.60, fill=True, lab_r=1.12, lab_off=-0.34, head=(live == "Om")) g = _reveal("inc", step, t) if g > 0: _angle(ax, zhat, nhat, g * inc, 0.36 * L, BLUE, r"$i$", alpha=1.0 if live == "inc" else 0.60, fill=True, lab_r=1.26, head=(live == "inc")) g = _reveal("om", step, t) if g > 0: _angle(ax, nhat, hhat, g * om, 0.22 * L, PURPLE, r"$\omega$", alpha=1.0 if live == "om" else 0.60, fill=True, lab_r=1.26, lab_off=+0.40, head=(live == "om")) # -- the node line, once there is one ------------------------------------- if _reveal("Om", step, t) >= 1.0: ax.plot(*np.vstack([-L * nhat, L * nhat]).T, color=GREEN, lw=1.6, ls=(0, (7, 4)), alpha=0.75) # -- the answer, with its components read off the inertial axes ----------- g = _reveal("state", step, t) if g > 0: r_vec = S["r_ijk"] ax.plot(*np.vstack([r_vec, [r_vec[0], r_vec[1], 0.0]]).T, color=MAROON, lw=1.0, ls=":", alpha=0.8 * g) ax.plot(*np.vstack([[0.0, 0.0, 0.0], [r_vec[0], 0.0, 0.0], [r_vec[0], r_vec[1], 0.0]]).T, color=MAROON, lw=1.0, ls=":", alpha=0.6 * g) # -- cosmetics ------------------------------------------------------------- zlo, zhi = -0.92 * L, 0.96 * L ax.set_xlim(-L, L) ax.set_ylim(-L, L) ax.set_zlim(zlo, zhi) ax.set_box_aspect((1.0, 1.0, (zhi - zlo) / (2.0 * L)), zoom=zoom) ax.set_facecolor(CREAM) ax.axis("off") ax.view_init(elev=elev, azim=azim) def _angle(ax, v_from, axis, ang, radius, color, label, alpha=1.0, fill=False, lab_r=1.20, lab_off=0.0, head=False): """Draw an angle where it is actually measured: an arc about `axis`. `lab_off` lifts the label off the arc along `axis`, i.e. along the normal of the plane the angle is swept in; three of these arcs are nearly concentric in projection, so a purely radial offset cannot separate them. `head` puts an arrowhead on the leading end, which is the only thing in the picture that shows the SENSE of the rotation. """ from mpl_toolkits.mplot3d.art3d import Poly3DCollection if abs(ang) < 1e-4: return n = max(8, int(90 * abs(ang) / (2 * np.pi)) + 8) pts = _arc(v_from, axis, ang, radius, n=n) if fill: ax.add_collection3d(Poly3DCollection( [_sector(v_from, axis, ang, radius, n=n)], facecolor=color, alpha=0.16 * alpha, edgecolor="none")) for end in (pts[0], pts[-1]): ax.plot(*np.vstack([np.zeros(3), end]).T, color=color, lw=0.9, alpha=0.75 * alpha) ax.plot(*pts.T, color=color, lw=2.0, alpha=alpha) if head and len(pts) > 2: tan = pts[-1] - pts[-2] nrm = np.linalg.norm(tan) if nrm > 0: ax.quiver(*pts[-1], *(0.28 * radius * tan / nrm), color=color, lw=2.0, arrow_length_ratio=0.9, alpha=alpha) nhat_ax = np.asarray(axis, dtype=float) nhat_ax = nhat_ax / np.linalg.norm(nhat_ax) ax.text(*(pts[len(pts) // 2] * lab_r + lab_off * radius * nhat_ax), label, color=color, fontsize=14, ha="center", va="center", alpha=alpha) # -- the running tally -------------------------------------------------------- def _vec(v, p=3): return "(" + ",".join(f"{x:9.{p}f}" for x in v) + ")" def panel_lines(S, step, t=1.0): """The lines of the side panel, as (text, colour, size, weight) tuples.""" deg = np.rad2deg out = [("GIVEN -- THE SIX ELEMENTS", DIM, 10.5, "bold"), (f"a = {S['a']:11.3f} km", INK, 10, "normal"), (f"e = {S['e']:11.6f}", PURPLE, 10, "normal"), (f"i = {deg(S['inc']):11.3f} deg", BLUE, 10, "normal"), (f"Ω = {deg(S['Om']):11.3f} deg", AMBER, 10, "normal"), (f"ω = {deg(S['om']):11.3f} deg", PURPLE, 10, "normal"), (f"f = {deg(S['f']):11.3f} deg", MAROON, 10, "normal")] out += [("", INK, 6, "normal"), ("STEP 1 -- THE CONIC IN ITS PLANE", DIM, 10.5, "bold")] if _reveal("pqw", step, t) > 0: out += [(f"p = a(1-e^2) = {S['p']:.3f} km", INK, 9.5, "normal"), (f"r = p/(1+e cos f) = {S['r']:.3f} km", INK, 9.5, "normal"), (f"r_PQW = {_vec(S['r_pqw'])} km", MAROON, 9.5, "normal"), (f"v_PQW = {_vec(S['v_pqw'], 6)} km/s", TEAL, 9.5, "normal")] out += [("", INK, 6, "normal"), ("STEP 2 -- Q = R3(Ω) R1(i) R3(ω)", DIM, 10.5, "bold")] fO, fi, fw = _fracs(step, t) done = [] if fO > 0: done.append((f"R3({deg(S['Om']) * fO:7.2f})", AMBER)) if fi > 0: done.append((f"R1({deg(S['inc']) * fi:7.2f})", BLUE)) if fw > 0: done.append((f"R3({deg(S['om']) * fw:7.2f})", PURPLE)) if done: out += [("applied so far, left to right:", DIM, 9.0, "normal")] out += [(" " + txt, col, 9.5, "normal") for txt, col in done] M = attitude(S["Om"], S["inc"], S["om"], (fO, fi, fw)) for k, row in enumerate(M): lead = "M = " if k == 1 else " " out += [(lead + "[" + " ".join(f"{x:+.4f}" for x in row) + "]", INK, 9.0, "normal")] if _reveal("state", step, t) > 0: out += [("", INK, 6, "normal"), ("RESULT", DIM, 10.5, "bold"), (f"r = {_vec(S['r_ijk'])} km", MAROON, 11, "bold"), (f"v = {_vec(S['v_ijk'], 6)} km/s", TEAL, 11, "bold")] if _reveal("state", step, t) >= 1.0: out += [(f"|r| = {np.linalg.norm(S['r_ijk']):.3f} km " f"(was {np.linalg.norm(S['r_pqw']):.3f})", DIM, 9.0, "normal"), (f"|v| = {np.linalg.norm(S['v_ijk']):.6f} km/s", DIM, 9.0, "normal"), ("a rotation cannot change a length", DIM, 9.0, "normal")] return out def draw_panel(ax, S, step, t=1.0): ax.clear() ax.set_xlim(0, 1) ax.set_ylim(0, 1) ax.axis("off") ax.set_facecolor(CREAM) y = 0.985 for text, color, size, weight in panel_lines(S, step, t): if text: ax.text(0.0, y, text, color=color, fontsize=size, weight=weight, family="DejaVu Sans Mono", ha="left", va="top", transform=ax.transAxes) y -= (size + 5.5) / 478.0 # -- frame schedule ----------------------------------------------------------- def schedule(draw=26, hold=16): """(step, t) for every frame: `draw` frames to build it, `hold` to read it.""" frames = [] for j in range(len(STEPS)): n_draw = 1 if STEP_KEYS[j] == "given" else draw for k in range(n_draw): frames.append((j, ease((k + 1) / n_draw))) frames += [(j, 1.0)] * hold frames += [(len(STEPS) - 1, 1.0)] * (2 * hold) return frames def _caption(step): """Captions carry mathtext, so the bold markers have to go.""" return STEPS[step][1].replace("**", "") def _figure(figsize=(12.6, 6.8), elev=24, azim=-12): import matplotlib.pyplot as plt plt.rcParams.update({"figure.facecolor": CREAM, "savefig.facecolor": CREAM, "axes.facecolor": CREAM, "text.color": INK, "font.size": 11}) fig = plt.figure(figsize=figsize) ax3 = fig.add_axes([-0.020, 0.050, 0.695, 0.905], projection="3d") ax3.view_init(elev=elev, azim=azim) axp = fig.add_axes([0.650, 0.050, 0.345, 0.905]) return fig, ax3, axp def render(path, S, title, fps=20, dpi=100, draw=26, hold=16, elev=24, azim=-12): import matplotlib matplotlib.use("Agg") from matplotlib import animation frames = schedule(draw=draw, hold=hold) fig, ax3, axp = _figure(elev=elev, azim=azim) cap = fig.text(0.5, 0.022, "", ha="center", fontsize=13.5, color=INK) fig.text(0.5, 0.972, title, ha="center", fontsize=12, color=DIM) def frame(k): step, t = frames[k] draw_scene(ax3, S, step, t) draw_panel(axp, S, step, t) cap.set_text(_caption(step)) return () ani = animation.FuncAnimation(fig, frame, frames=len(frames), blit=False) os.makedirs(os.path.dirname(os.path.abspath(path)), exist_ok=True) writer = "ffmpeg" if path.endswith(".mp4") else "pillow" print(f"writing {path} ({len(frames)} frames) ...") ani.save(path, writer=writer, fps=fps, dpi=dpi, savefig_kwargs=dict(facecolor=CREAM)) print("done:", path, f"{os.path.getsize(path) / 1e6:.2f} MB") def save_png(path, S, title, step=None, dpi=160, elev=24, azim=-12): import matplotlib matplotlib.use("Agg") step = len(STEPS) - 1 if step is None else step fig, ax3, axp = _figure(elev=elev, azim=azim) fig.text(0.5, 0.972, title, ha="center", fontsize=12, color=DIM) fig.text(0.5, 0.022, _caption(step), ha="center", fontsize=13.5, color=INK) draw_scene(ax3, S, step, 1.0) draw_panel(axp, S, step, 1.0) os.makedirs(os.path.dirname(os.path.abspath(path)), exist_ok=True) fig.savefig(path, dpi=dpi, facecolor=CREAM) print("wrote", path) def storyboard(path, S, title, dpi=150, elev=24, azim=-12, cols=3): """All six steps on one sheet -- the animation, for a printed handout.""" import matplotlib matplotlib.use("Agg") import matplotlib.pyplot as plt plt.rcParams.update({"figure.facecolor": CREAM, "savefig.facecolor": CREAM, "axes.facecolor": CREAM, "text.color": INK}) rows = int(np.ceil(len(STEPS) / cols)) fig = plt.figure(figsize=(5.0 * cols, 3.45 * rows)) fig.text(0.5, 0.996, title, ha="center", va="top", fontsize=13, color=DIM) for j in range(len(STEPS)): ax = fig.add_subplot(rows, cols, j + 1, projection="3d") ax.view_init(elev=elev, azim=azim) draw_scene(ax, S, j, 1.0, zoom=1.46) ax.set_title(_caption(j), color=INK, fontsize=9.0, pad=-4, wrap=True) fig.subplots_adjust(left=0.005, right=0.995, top=0.925, bottom=0.005, wspace=0.01, hspace=0.14) fig.savefig(path, dpi=dpi, facecolor=CREAM, bbox_inches="tight") print("wrote", path) def interactive(S, title, elev=24, azim=-12): """Step through it by hand: left/right arrows, or the two buttons.""" import matplotlib.pyplot as plt from matplotlib.widgets import Button fig, ax3, axp = _figure(elev=elev, azim=azim) fig.text(0.5, 0.972, title, ha="center", fontsize=12, color=DIM) cap = fig.text(0.5, 0.022, "", ha="center", fontsize=13.5, color=INK) state = {"step": 0} def show(): draw_scene(ax3, S, state["step"], 1.0) draw_panel(axp, S, state["step"], 1.0) cap.set_text(_caption(state["step"])) fig.canvas.draw_idle() def go(d): def _cb(_=None): state["step"] = int(np.clip(state["step"] + d, 0, len(STEPS) - 1)) show() return _cb bprev = Button(fig.add_axes([0.30, 0.955, 0.06, 0.038]), "< back", color="#eee9dd", hovercolor="#e2dbc9") bnext = Button(fig.add_axes([0.37, 0.955, 0.06, 0.038]), "next >", color="#eee9dd", hovercolor="#e2dbc9") bprev.on_clicked(go(-1)) bnext.on_clicked(go(+1)) fig.canvas.mpl_connect("key_press_event", lambda ev: ( go(+1)() if ev.key in ("right", " ") else go(-1)() if ev.key == "left" else None)) show() plt.show() # -- verification ------------------------------------------------------------- def verify(S): """The three matrices, their product, and the round trip back to Part III.""" deg = np.rad2deg np.set_printoptions(formatter={"float": lambda x: f"{x:+9.6f}"}) print("forward map (a, e, i, Omega, omega, f) -> (r, v)\n") print(f" given a = {S['a']:.4f} km e = {S['e']:.6f} " f"i = {deg(S['inc']):.4f}") print(f" Omega = {deg(S['Om']):.4f} omega = {deg(S['om']):.4f}" f" f = {deg(S['f']):.4f} (deg)\n") print(f" step 1 p = {S['p']:.4f} km, r = {S['r']:.4f} km") print(f" r_PQW = {_vec(S['r_pqw'], 4)} km") print(f" v_PQW = {_vec(S['v_pqw'], 6)} km/s\n") print(" step 2 the three rotations, applied left to right:\n") for nm, Rm in (("R3(Omega)", R3(S["Om"])), ("R1(i)", R1(S["inc"])), ("R3(omega)", R3(S["om"]))): for k, row in enumerate(Rm): print(f" {nm if k == 1 else '':>10} " f"[{' '.join(f'{x:+9.6f}' for x in row)}]") print() print(" Q = R3(Omega) R1(i) R3(omega):\n") for k, row in enumerate(S["Q"]): print(f" {'Q =' if k == 1 else '':>10} " f"[{' '.join(f'{x:+9.6f}' for x in row)}]") Q = S["Q"] print(f"\n orthogonality max|Q Q^T - 1| = " f"{np.abs(Q @ Q.T - np.eye(3)).max():.3e}") print(f" det Q = {np.linalg.det(Q):+.12f} (a rotation, not a " f"reflection)") print(f" third column = h-hat = {_vec(Q[:, 2], 6)} " f"(carries no omega, as it must not)") print(f"\n result r = {_vec(S['r_ijk'], 4)} km") print(f" v = {_vec(S['v_ijk'], 6)} km/s") print(f" |r| = {np.linalg.norm(S['r_ijk']):.6f} km " f"(perifocal: {np.linalg.norm(S['r_pqw']):.6f})") print(f" |v| = {np.linalg.norm(S['v_ijk']):.6f} km/s " f"(perifocal: {np.linalg.norm(S['v_pqw']):.6f})") el = elements_from_state(S["r_ijk"], S["v_ijk"], S["mu"]) err = max(abs(el["a"] - S["a"]), abs(el["e"] - S["e"]), *(abs((el[k] - S[k] + np.pi) % (2 * np.pi) - np.pi) for k in ("inc", "Om", "om", "f"))) print(f"\n round trip back through Part III: worst element error " f"{err:.3e}") # the intrinsic reading animated here, against the extrinsic reading of the # deck's sentence: the same Q, built the other way round zhat, xhat = np.array([0.0, 0.0, 1.0]), np.array([1.0, 0.0, 0.0]) n_ax = R3(S["Om"]) @ xhat h_ax = R3(S["Om"]) @ R1(S["inc"]) @ zhat Q_intrinsic = (rodrigues(h_ax, S["om"]) @ rodrigues(n_ax, S["inc"]) @ rodrigues(zhat, S["Om"])) print(f" intrinsic (about the frame's own moving axes) vs the product: " f"max diff {np.abs(Q_intrinsic - Q).max():.3e}") np.set_printoptions() def main(): ap = argparse.ArgumentParser(description=__doc__.split("\n")[3]) ap.add_argument("--example", choices=sorted(EXAMPLES), default="4", help="which worked example of the deck to compose") ap.add_argument("--elements", nargs=6, type=float, metavar=("A", "E", "I", "OM", "OMEGA", "F"), help="a, e in km and -, then i, Omega, omega, f in degrees") ap.add_argument("--mu", type=float, default=MU_EARTH) ap.add_argument("--verify", action="store_true", help="the three matrices, the product, and the round trip") ap.add_argument("--png", metavar="PATH", help="save the final frame") ap.add_argument("--storyboard", metavar="PATH", help="save all six steps on one sheet") ap.add_argument("--save", metavar="PATH", help="render the animation to PATH (.gif or .mp4)") ap.add_argument("--step", type=int, default=None, help="which step --png should show (0-5, default: last)") ap.add_argument("--fps", type=int, default=20) ap.add_argument("--dpi", type=int, default=100) ap.add_argument("--draw", type=int, default=26, help="frames per rotation") ap.add_argument("--hold", type=int, default=16, help="frames per pause") ap.add_argument("--elev", type=float, default=24.0) ap.add_argument("--azim", type=float, default=-12.0) args = ap.parse_args() if args.elements: a, e, inc, Om, om, f = args.elements prm = dict(a=a, e=e, inc=inc, Om=Om, om=om, f=f) title = (rf"$a={a:.1f}$ km, $e={e:.4f}$, $i={inc:.2f}^\circ$, " rf"$\Omega={Om:.2f}^\circ$, $\omega={om:.2f}^\circ$, " rf"$f={f:.2f}^\circ$") else: prm = dict(EXAMPLES[args.example]) title = prm.pop("title") S = compose(prm, mu=args.mu) if args.verify: verify(S) if args.png: save_png(args.png, S, title, step=args.step, elev=args.elev, azim=args.azim) if args.storyboard: storyboard(args.storyboard, S, title, elev=args.elev, azim=args.azim) if args.save: render(args.save, S, title, fps=args.fps, dpi=args.dpi, draw=args.draw, hold=args.hold, elev=args.elev, azim=args.azim) if not (args.verify or args.png or args.storyboard or args.save): interactive(S, title, elev=args.elev, azim=args.azim) if __name__ == "__main__": main()