"""Rotating the apse line the long way round (section "Going round instead" of the manoeuvre notes). The orbit of the example, a = 12 000 km and e = 0.4, is to be turned through eta = 90 deg without changing its size or shape. The two-burn route: fly the initial orbit to an apse; burn 1 circularises it there; coast round the circle through eta; burn 2 puts the craft on the final orbit, which is the initial one turned through eta. Two versions, with the burns at apoapsis (the cheaper route) and at periapsis (the faster one). The craft moves at its true speed (fast near periapsis, slow near apoapsis, constant on the circle), only slowed down, and both animations run on the same clock; the burns are instantaneous and are shown as pauses. Run from the sfm-slides directory: python3 codes/anim-man-apse-detour.py Writes img/anim-man-apse-detour.gif (apoapsis) and img/anim-man-apse-detour-peri.gif (periapsis). Afterwards run python3 _filters/gif-frames.py (from the project root) to write the still frames that the PDF output uses. """ import numpy as np import _figfont # all figure text at one scale import matplotlib.pyplot as plt from matplotlib.patches import Circle from _manstyle import CREAM, INK, MAROON, BLUE, GREEN, AMBER from _gif import save_gif MU, RE = 398600.4418, 6378.137 A, E, ETA = 12000.0, 0.40, np.radians(90.0) B = A*np.sqrt(1 - E*E) RP, RA = A*(1 - E), A*(1 + E) vis = lambda r: np.sqrt(MU*(2/r - 1/A)) # speed on either ellipse (same a) vcirc = lambda r: np.sqrt(MU/r) T = 2*np.pi*np.sqrt(A**3/MU) # period of either ellipse TH1, TH2 = np.pi, np.pi + ETA # polar angles of burn 1 and burn 2 (the same in both versions) def kepler_E(M): Ea = M + E*np.sin(M) for _ in range(60): Ea = Ea - (Ea - E*np.sin(Ea) - M)/(1 - E*np.cos(Ea)) return Ea def ellipse(M, w): """Position on the orbit with periapsis in direction w, at mean anomaly M (array).""" Ea = kepler_E(np.asarray(M, dtype=float)) x, y = A*(np.cos(Ea) - E), B*np.sin(Ea) return np.stack([x*np.cos(w) - y*np.sin(w), x*np.sin(w) + y*np.cos(w)], axis=-1) def circle(theta, R): return R*np.stack([np.cos(theta), np.sin(theta)], axis=-1) def tangent(theta): # counter-clockwise direction of motion on a circle return np.array([-np.sin(theta), np.cos(theta)]) def smooth(u): return u*u*(3 - 2*u) # The two versions. Mb: mean anomaly of the burn point (pi = apoapsis, 0 = periapsis). The initial orbit's periapsis is # placed at w0 = TH1 - Mb, so that burn 1 is at polar angle TH1 in both: the burns sit at the left and at the bottom. # All text positions are in data units; "lead" puts a thin leader line from a label to its object. CASES = { "apo": dict( file="img/anim-man-apse-detour.gif", Mb=np.pi, R=RA, other="periapsis", rsym="a", title2="$a=12\\,000$ km, $e=0.4$: both burns at apoapsis, $r_a=" + f"{RA:,.0f}".replace(",", "\\,") + "$ km", xlim=(-27500, 17000), ylim=(-21500, 14500), status=(0.011, 0.986, "left", "top"), clock=(0.011, 0.933, "left", "top"), init=((-4800, 11500), "center", "bottom"), fin=((11500, -4800), "left", "center"), circ=((-13600, -13600), "right", "top", False), apse=(((RP + 400, 700), "left", "bottom"), ((0, RP + 600), "center", "bottom")), eta_r=(8000, 10200), b1=((-18300, -2140), "right", "center", False), b2=((-2140, -18500), "center", "top", False)), "peri": dict( file="img/anim-man-apse-detour-peri.gif", Mb=0.0, R=RP, other="apoapsis", rsym="p", title2="$a=12\\,000$ km, $e=0.4$: both burns at periapsis, $r_p=" + f"{RP:,.0f}".replace(",", "\\,") + "$ km", xlim=(-20500, 24500), ylim=(-18000, 18500), status=(0.99, 0.085, "right", "bottom"), clock=(0.99, 0.012, "right", "bottom"), init=((4800, 11500), "center", "bottom"), fin=((-11500, 4800), "right", "center"), circ=((-14500, -7500), "right", "top", True), apse=(((RA + 500, 700), "left", "bottom"), ((0, RA + 600), "center", "bottom")), eta_r=(11500, 13000), b1=((-12800, -500), "right", "center", True), b2=((0, -12900), "center", "top", True)), } def render(key): c = CASES[key] Mb, R = c["Mb"], c["R"] w0, w1 = TH1 - Mb, TH1 - Mb + ETA # periapsis directions of the initial and the final orbit V_O, V_C = vis(R), vcirc(R) # orbit speed and circular speed at the burns D1 = V_C - V_O # burn 1 (burn 2 is -D1); > 0 prograde, < 0 retrograde T_COAST = ETA*np.sqrt(R**3/MU) B1, B2 = circle(TH1, R), circle(TH2, R) # self-checks of the geometry the picture relies on assert np.allclose(ellipse(Mb, w0), B1) and np.allclose(ellipse(Mb, w1), B2) r2v, v2v = B2, V_O*tangent(TH2) # state just after burn 2: orbit speed, tangential ev = ((v2v @ v2v - MU/R)*r2v - (r2v @ v2v)*v2v)/MU # eccentricity vector of the orbit it defines assert np.isclose(np.hypot(*ev), E) and np.isclose(np.cos(np.arctan2(ev[1], ev[0]) - w1), 1.0) # ---- timeline: the same clock rate in both versions ------------------------------------------ N1 = 56 # frames for half a revolution of an ellipse N2 = max(8, int(round(N1*T_COAST/(T/2)))) # the coast, at the same clock rate NB, NE = 16, 26 # frames for a burn, and for the final hold timeline = ([("fly1", k/(N1 - 1)) for k in range(N1)] + [("burn1", k/(NB - 1)) for k in range(NB)] + [("coast", k/(N2 - 1)) for k in range(N2)] + [("burn2", k/(NB - 1)) for k in range(NB)] + [("fly2", k/(N1 - 1)) for k in range(N1)] + [("end", 1.0)]*NE) ORDER = ["fly1", "burn1", "coast", "burn2", "fly2", "end"] LABEL = {"fly1": "1. coasting on the initial orbit", "burn1": "2. burn 1: circularise", "coast": r"3. coasting on the circle, through $\eta$", "burn2": "4. burn 2: back to the orbit's shape", "fly2": "5. on the final orbit", "end": "5. on the final orbit"} # ---- the picture ------------------------------------------------------------------------------- fig, ax = plt.subplots(figsize=(9.6, 8.0)) ax.set_aspect("equal"); ax.axis("off") ax.set_xlim(*c["xlim"]); ax.set_ylim(*c["ylim"]) ax.add_patch(Circle((0, 0), RE, color=BLUE, alpha=0.30, zorder=1)) ax.plot(0, 0, "o", color=BLUE, ms=5, zorder=2) m = np.linspace(0, 2*np.pi, 500) th = np.linspace(TH1, TH2, 200) dash = dict(lw=1.3, ls=(0, (4, 4)), alpha=0.5, zorder=2) ax.plot(*ellipse(m, w0).T, color=BLUE, **dash) ax.plot(*ellipse(m, w1).T, color=GREEN, **dash) ax.plot(*circle(th, R).T, color=MAROON, **dash) for w, col in ((w0, BLUE), (w1, GREEN)): # the two apse lines, and the angle between them q1, q2 = RP*np.array([np.cos(w), np.sin(w)]), -RA*np.array([np.cos(w), np.sin(w)]) ax.plot([q1[0], q2[0]], [q1[1], q2[1]], color=col, lw=0.9, ls=":", alpha=0.7, zorder=2) ra_, rl_ = c["eta_r"] arc = np.linspace(TH1, TH2, 60) ax.plot(ra_*np.cos(arc), ra_*np.sin(arc), color=INK, lw=1.2, zorder=3) ax.text(rl_*np.cos((TH1 + TH2)/2), rl_*np.sin((TH1 + TH2)/2), r"$\eta=90^\circ$", fontsize=12, ha="center", va="center") (px, py), ha, va = c["init"]; ax.text(px, py, "initial orbit", color=BLUE, fontsize=12, ha=ha, va=va) (px, py), ha, va = c["fin"]; ax.text(px, py, "final orbit\n(turned by $\\eta$)", color=GREEN, fontsize=12, ha=ha, va=va) ctxt = f"circular orbit\n$r_{c['rsym']}=" + f"{R:,.0f}".replace(",", "\\,") + "$ km" (px, py), ha, va, lead = c["circ"] if lead: ax.annotate(ctxt, xy=circle((TH1 + TH2)/2, R), xytext=(px, py), color=MAROON, fontsize=11, ha=ha, va=va, arrowprops=dict(arrowstyle="-", color=MAROON, lw=0.8, alpha=0.6, shrinkA=3, shrinkB=0)) else: ax.text(px, py, ctxt, color=MAROON, fontsize=11, ha=ha, va=va) for (pos, ha, va), col in zip(c["apse"], (BLUE, GREEN)): # the other apse of each orbit ax.text(pos[0], pos[1], c["other"], color=col, fontsize=9, ha=ha, va=va, alpha=0.8) ax.text(0, -1500, "Earth", color=BLUE, fontsize=9, ha="center", va="top") trails = {k: ax.plot([], [], color=col, lw=3.2, zorder=5, solid_capstyle="round")[0] for k, col in (("fly1", BLUE), ("coast", MAROON), ("fly2", GREEN))} dot, = ax.plot([], [], "o", color=INK, ms=10, mec=CREAM, mew=1.5, zorder=9) # direction-of-travel arrowheads, one on each trail, shown once the craft has passed them heads = [] for k, pts, col in (("fly1", ellipse([Mb - np.pi/2 - 0.04, Mb - np.pi/2 + 0.04], w0), BLUE), ("coast", circle(np.array([(TH1 + TH2)/2 - 0.03, (TH1 + TH2)/2 + 0.03]), R), MAROON), ("fly2", ellipse([Mb + np.pi/2 - 0.04, Mb + np.pi/2 + 0.04], w1), GREEN)): ar = ax.annotate("", xy=pts[1], xytext=pts[0], zorder=6, annotation_clip=False, arrowprops=dict(arrowstyle="-|>", color=col, lw=2.4, mutation_scale=20, shrinkA=0, shrinkB=0)) ar.set_visible(False); heads.append((k, ar)) # the burns: a star on the path, a flash while it fires, the Delta-v arrow (along or against the velocity) LEN = 3900.0*abs(D1) burns = {} for k, P, th_, sg, nm, cfg in (("burn1", B1, TH1, np.sign(D1), "1", c["b1"]), ("burn2", B2, TH2, -np.sign(D1), "2", c["b2"])): txt = (f"burn {nm}\n$\\Delta v_{nm}={'+' if sg > 0 else '-'}{abs(D1):.3f}$ km/s\n" f"({'prograde' if sg > 0 else 'retrograde'})") tpos, ha, va, lead = cfg star, = ax.plot([], [], "*", color=AMBER, mec=MAROON, ms=19, zorder=8) ring = Circle(P, 1.0, fill=False, color=AMBER, lw=2.5, zorder=7); ring.set_visible(False); ax.add_patch(ring) d = sg*tangent(th_) arr = ax.annotate("", xy=P + 1e-3*d, xytext=P, zorder=8, arrowprops=dict(arrowstyle="-|>", color=AMBER, lw=3.4, mutation_scale=24, shrinkA=0, shrinkB=0)) arr.set_visible(False) if lead: label = ax.annotate(txt, xy=P, xytext=tpos, color=MAROON, fontsize=11, ha=ha, va=va, visible=False, arrowprops=dict(arrowstyle="-", color=MAROON, lw=0.8, alpha=0.6, shrinkA=3, shrinkB=9)) else: label = ax.text(*tpos, txt, color=MAROON, fontsize=11, ha=ha, va=va, visible=False) burns[k] = (P, d, star, ring, arr, label) sx, sy, sha, sva = c["status"]; cx, cy, cha, cva = c["clock"] status = ax.text(sx, sy, "", fontsize=13, va=sva, ha=sha, color=INK, weight="bold", transform=ax.transAxes) clock = ax.text(cx, cy, "", fontsize=12, va=cva, ha=cha, family="monospace", color=INK, transform=ax.transAxes) ax.set_title(r"Rotating the apse line by $\eta=90^\circ$ the long way round: circularise, coast, re-eccentricise" "\n" + c["title2"], fontsize=12.5) fig.tight_layout() def frame(i): ph, u = timeline[min(i, len(timeline) - 1)] k = ORDER.index(ph) prog = lambda name: 1.0 if k > ORDER.index(name) else (u if ph == name else 0.0) # position, speed, elapsed time, Delta-v spent so far if ph == "fly1": pos = ellipse(Mb - np.pi + np.pi*u, w0); v, t, spent = vis(np.linalg.norm(pos)), u*T/2, 0.0 elif ph == "burn1": pos, v, t, spent = B1, V_O + D1*smooth(u), T/2, abs(D1)*smooth(u) elif ph == "coast": pos, v, t, spent = circle(TH1 + ETA*u, R), V_C, T/2 + u*T_COAST, abs(D1) elif ph == "burn2": pos, v, t, spent = B2, V_C - D1*smooth(u), T/2 + T_COAST, abs(D1)*(1 + smooth(u)) else: uu = u if ph == "fly2" else 1.0 pos = ellipse(Mb + np.pi*uu, w1); v = vis(np.linalg.norm(pos)); t = T/2 + T_COAST + uu*T/2; spent = 2*abs(D1) dot.set_data([pos[0]], [pos[1]]) s = prog("fly1"); p = ellipse(Mb - np.pi + np.pi*np.linspace(0, s, max(2, int(200*s) + 2)), w0); trails["fly1"].set_data(*p.T) s = prog("coast"); p = circle(TH1 + ETA*np.linspace(0, s, max(2, int(100*s) + 2)), R); trails["coast"].set_data(*p.T) s = prog("fly2"); p = ellipse(Mb + np.pi*np.linspace(0, s, max(2, int(200*s) + 2)), w1); trails["fly2"].set_data(*p.T) for name, ar in heads: ar.set_visible(prog(name) >= 0.53) for name, (P, d, star, ring, arr, label) in burns.items(): kb = ORDER.index(name) star.set_data([P[0]] if k >= kb else [], [P[1]] if k >= kb else []) label.set_visible(k >= kb) uu = 1.0 if k > kb else (smooth(u) if ph == name else 0.0) arr.set_visible(uu > 0.02) if uu > 0.02: arr.xy = P + LEN*uu*d ring.set_visible(ph == name) if ph == name: ring.set_radius(500 + 2600*u); ring.set_alpha(1 - 0.8*u) status.set_text(LABEL[ph]) clock.set_text(f"t = {t/3600:5.2f} h v = {v:5.3f} km/s\n" f"Δv spent = {spent:5.3f} km/s") frame(len(timeline) - 1) return fig, frame, timeline, c["file"] if __name__ == "__main__": for key in CASES: fig, frame, timeline, path = render(key) save_gif(fig, frame, len(timeline), 20, path, dpi=90, colors=64)