"""Crossing coaxial ellipses: one impulse at the intersection, or two at the apsides. Orbit 1 dips lower AND reaches higher than orbit 2, so the two conics cross. That opens a route the non-crossing cases do not have -- a single impulse at a crossing point -- but it has to rotate the velocity vector, and rotation is expensive. """ import numpy as np import matplotlib.pyplot as plt from _manstyle import CREAM, INK, MAROON, BLUE, GREEN, AMBER, PURPLE, despine mu = 398600.4418 rp1, ra1 = 8000.0, 30000.0 rp2, ra2 = 12000.0, 24000.0 a1, e1 = 0.5*(rp1 + ra1), (ra1 - rp1)/(ra1 + rp1) a2, e2 = 0.5*(rp2 + ra2), (ra2 - rp2)/(ra2 + rp2) p1, p2 = a1*(1 - e1**2), a2*(1 - e2**2) vis = lambda r, a: np.sqrt(mu*(2/r - 1/a)) # ---- the crossing ------------------------------------------------------- cosf = (p2 - p1)/(p1*e2 - p2*e1) fx = np.arccos(cosf) rx = p1/(1 + e1*cosf) def polar_v(p, e, f): h = np.sqrt(mu*p) return h*e*np.sin(f)/p, h*(1 + e*np.cos(f))/p # (v_r, v_theta) vr1, vt1 = polar_v(p1, e1, fx) vr2, vt2 = polar_v(p2, e2, fx) dv_single = np.hypot(vr2 - vr1, vt2 - vt1) # ---- the two-impulse apsidal alternatives ------------------------------- def apsidal(rA, aFrom, rB, aTo): at = 0.5*(rA + rB) d1 = abs(vis(rA, at) - vis(rA, aFrom)) d2 = abs(vis(rB, aTo) - vis(rB, at)) return at, d1, d2, d1 + d2, np.pi*np.sqrt(at**3/mu)/3600 atP1, e1P1, e2P1, totP1, tP1 = apsidal(rp1, a1, ra2, a2) atP2, d1P2, d2P2, totP2, tP2 = apsidal(ra1, a1, rp2, a2) # ---- figure ------------------------------------------------------------- fig, ax = plt.subplots(figsize=(9.4, 7.6)) f = np.linspace(0, 2*np.pi, 900) for p, e, col, lab in ((p1, e1, BLUE, f"orbit 1: $a={a1:.0f}$ km, $e={e1:.4f}$"), (p2, e2, GREEN, f"orbit 2: $a={a2:.0f}$ km, $e={e2:.4f}$")): r = p/(1 + e*np.cos(f)) ax.plot(r*np.cos(f), r*np.sin(f), color=col, lw=2.2, label=lab) # the P2 transfer: apoapsis 30000 on -x, periapsis 12000 on +x, flown lower half etP2 = (ra1 - rp2)/(ra1 + rp2) ft = np.linspace(np.pi, 2*np.pi, 400) rt = atP2*(1 - etP2**2)/(1 + etP2*np.cos(ft)) ax.plot(rt*np.cos(ft), rt*np.sin(ft), color=MAROON, lw=3.0, zorder=5, label=f"two-impulse: $a_t={atP2:.0f}$ km, {totP2:.3f} km/s, {tP2:.2f} h") ax.plot(0, 0, "o", color=INK, ms=8, zorder=9) ax.plot([-ra1, rp1], [0, 0], color=INK, lw=0.8, ls=":", zorder=2) # the two crossing points for sgn in (+1, -1): ax.plot(rx*np.cos(sgn*fx), rx*np.sin(sgn*fx), "o", color=PURPLE, ms=12, zorder=10) ax.annotate(f"the orbits cross at $r={rx:.0f}$ km,\n$f=\\pm{np.degrees(fx):.2f}^\\circ$ " "--- and both crossings cost the same", xy=(rx*np.cos(-fx), rx*np.sin(-fx)), xytext=(6000, -24000), color=PURPLE, fontsize=10.5, ha="center", arrowprops=dict(arrowstyle="-|>", color=PURPLE, lw=1.3)) # ---- velocity triangle, drawn as an inset clear of the conics ------------ X, Y = rx*np.cos(fx), rx*np.sin(fx) rhat = np.array([np.cos(fx), np.sin(fx)]) that = np.array([-np.sin(fx), np.cos(fx)]) VS = 2400.0 # km drawn per km/s V1 = (vr1*rhat + vt1*that)*VS V2 = (vr2*rhat + vt2*that)*VS OX, OY = 26500.0, 13000.0 # inset origin, clear of the conics ax.plot(OX, OY, "o", color=PURPLE, ms=7, zorder=9) for V, col, lab, off in ((V1, BLUE, r"$\mathbf{v}_1$" + f" {np.hypot(vr1,vt1):.3f}", (-1200, 2200)), (V2, GREEN, r"$\mathbf{v}_2$" + f" {np.hypot(vr2,vt2):.3f}", (-900, -2800))): ax.annotate("", xy=(OX + V[0], OY + V[1]), xytext=(OX, OY), arrowprops=dict(arrowstyle="-|>", color=col, lw=2.6), zorder=8) ax.text(OX + V[0] + off[0], OY + V[1] + off[1], lab, color=col, fontsize=12.5, ha="center") ax.annotate("", xy=(OX + V2[0], OY + V2[1]), xytext=(OX + V1[0], OY + V1[1]), arrowprops=dict(arrowstyle="-|>", color=AMBER, lw=3.0), zorder=9) ax.text(OX + 0.5*(V1[0] + V2[0]) + 4600, OY + 0.5*(V1[1] + V2[1]) - 700, f"$\\Delta v={dv_single:.3f}$", color=AMBER, fontsize=13, ha="left", va="center") ax.text(OX - 5200, OY + 6400, "velocity triangle at the crossing\n" f"turn {abs(np.degrees(np.arctan2(vr2,vt2))-np.degrees(np.arctan2(vr1,vt1))):.1f}" "$^\\circ$, speed " f"{np.hypot(vr2,vt2)-np.hypot(vr1,vt1):+.3f} km/s", color=INK, fontsize=10.5, ha="center", va="bottom") # the two apsidal burns: prograde at apoapsis(1), retrograde at periapsis(2) -- both -y for x, lab in ((-ra1, r"$\Delta v_1$" + f" {d1P2:.3f}"), (rp2, r"$\Delta v_2$" + f" {d2P2:.3f}")): ax.annotate("", xy=(x, -5600), xytext=(x, -900), arrowprops=dict(arrowstyle="-|>", color=MAROON, lw=2.6), zorder=9) ax.text(x + 1300, -4000, lab, color=MAROON, fontsize=11.5, va="center") ax.plot([-ra1, rp2], [0, 0], "o", color=MAROON, ms=9, zorder=9) ax.text(4000, -30500, f"one impulse at a crossing: {dv_single:.3f} km/s, immediate", color=AMBER, fontsize=12.5, ha="center") ax.text(4000, -34500, f"two impulses at the apsides: {totP2:.3f} km/s, " f"after a {tP2:.2f} h coast", color=MAROON, fontsize=12.5, ha="center") ax.text(4000, -38800, f"the single burn costs {dv_single/totP2:.2f}$\\times$ as much", color=INK, fontsize=13.5, ha="center", weight="bold") ax.set_aspect("equal"); ax.axis("off") ax.set_xlim(-36000, 44000); ax.set_ylim(-42000, 26000) ax.legend(frameon=False, fontsize=9.5, loc="upper left") ax.set_title("Crossing coaxial ellipses: one impulse, or two", fontsize=13) fig.tight_layout() fig.savefig("img/fig-man-coaxial-cross.png", dpi=160, bbox_inches="tight") plt.close(fig) print("wrote img/fig-man-coaxial-cross.png") print(f"\ncrossing: r={rx:.3f} km f=+/-{np.degrees(fx):.4f} deg cos f={cosf:.6f}") print(f" orbit 1 at crossing: vr={vr1:.5f} vt={vt1:.5f} |v|={np.hypot(vr1,vt1):.5f} " f"gamma={np.degrees(np.arctan2(vr1,vt1)):.3f} deg") print(f" orbit 2 at crossing: vr={vr2:.5f} vt={vt2:.5f} |v|={np.hypot(vr2,vt2):.5f} " f"gamma={np.degrees(np.arctan2(vr2,vt2)):.3f} deg") print(f" single impulse dv = {dv_single:.5f} km/s") print(f" apsidal P1 rp1->ra2: a_t={atP1:.0f} total={totP1:.5f} t={tP1:.3f} h") print(f" apsidal P2 ra1->rp2: a_t={atP2:.0f} dv1={d1P2:.5f} dv2={d2P2:.5f} " f"total={totP2:.5f} t={tP2:.3f} h") print(f" ratio single/two = {dv_single/totP2:.4f}")