"""The stopover: three impulses with r_b *short* of the target (R* < R). It looks like a bi-elliptic transfer and it is not. With r_b < r_2 the second ellipse runs from its *periapsis* r_b out to its *apoapsis* r_2, so * burn 2 must raise apoapsis, i.e. buy energy, out where v is smallest; * the spacecraft arrives at r_2 below circular speed, so burn 3 is prograde, dv3/vc1 = 1/sqrt(R) - sqrt(2R*/(R(R+R*))) > 0 (R* < R), which is the bi-elliptic dv3 with its sign reversed -- adding the published formulas as they stand therefore *understates* the cost by exactly 2 dv3. All three burns are prograde and the manoeuvre never beats Hohmann: at fixed R the cost is a hump over R* in (1, R) with both ends equal to the Hohmann cost, because both ends *are* the Hohmann transfer (R*=1 burns at r1 itself, R*=R circularises on arrival). Past R*=R the same curve becomes the genuine bi-elliptic family and drops below Hohmann whenever R > 15.58. Run from the sfm-slides directory: python3 codes/fig-man-stopover.py """ import numpy as np import matplotlib.pyplot as plt from _manstyle import CREAM, INK, MAROON, BLUE, GREEN, AMBER, PURPLE, despine BOX = dict(facecolor=CREAM, edgecolor="none", pad=1.6) R_CASE = 30.0 # radius ratio used on the slide and in the example RS_EX, RS_BI = 25.0, 40.0 # the worked stopover, and a genuine bi-elliptic def hoh(R): """Hohmann cost, normalised by v_c1.""" return np.sqrt(2*R/(1 + R)) - 1 + (1/np.sqrt(R))*(1 - np.sqrt(2/(1 + R))) def three(R, Rs): """Three-impulse cost via r_b = R* r_1, every magnitude taken positive. Valid on both sides of R* = R: for R* > R it is the bi-elliptic transfer, for R* < R the stopover. The two agree with Hohmann at R* = R. """ d1 = np.sqrt(2*Rs/(1 + Rs)) - 1 d2 = abs(np.sqrt(2*R/(Rs*(R + Rs))) - np.sqrt(2/(Rs*(1 + Rs)))) d3 = abs(1/np.sqrt(R) - np.sqrt(2*Rs/(R*(R + Rs)))) return d1 + d2 + d3 # -------------------------------------------------------------------------- # Figure # -------------------------------------------------------------------------- fig, (axA, axB) = plt.subplots(1, 2, figsize=(13.8, 5.9), gridspec_kw=dict(width_ratios=(1.0, 1.12))) # ---- (a) geometry, drawn for r1 : r_b : r2 = 1 : 3.5 : 5 ------------------ r1, rb, r2 = 1.0, 3.5, 5.0 th = np.linspace(0, 2*np.pi, 700) axA.plot(r1*np.cos(th), r1*np.sin(th), color=BLUE, lw=1.8) axA.plot(r2*np.cos(th), r2*np.sin(th), color=GREEN, lw=1.8) axA.plot(rb*np.cos(th), rb*np.sin(th), color=INK, lw=0.9, ls=(0, (4, 4)), alpha=0.45) def half(rp, ra, rot, **kw): """Half an ellipse, from periapsis at position angle `rot`, going CCW.""" a, e = 0.5*(rp + ra), (ra - rp)/(ra + rp) f = np.linspace(0.0, np.pi, 400) r = a*(1 - e*e)/(1 + e*np.cos(f)) axA.plot(r*np.cos(f + rot), r*np.sin(f + rot), **kw) half(r1, rb, 0.0, color=MAROON, lw=2.8, zorder=5) # ellipse 1: r1 -> r_b half(rb, r2, np.pi, color=PURPLE, lw=2.8, zorder=5) # ellipse 2: r_b -> r2 def burn(xy, up, lab, dx=0.28): """Prograde impulse; `up` is the sign of the motion along +y.""" s = 1.55 if up else -1.55 axA.annotate("", xy=(xy[0], xy[1] + s), xytext=(xy[0], xy[1] + 0.20*np.sign(s)), arrowprops=dict(arrowstyle="-|>", color=AMBER, lw=2.6), zorder=8) axA.text(xy[0] + dx, xy[1] + 1.05*np.sign(s), lab, color=AMBER, fontsize=13, zorder=8, bbox=BOX) burn((r1, 0.0), True, r"$\Delta v_1$") # at r1, moving +y burn((-rb, 0.0), False, r"$\Delta v_2$", dx=-1.30) # at r_b, moving -y burn((r2, 0.0), True, r"$\Delta v_3$") # at r2, moving +y axA.plot(0.0, 0.0, "o", color=INK, ms=8, zorder=9) for x, col in ((r1, BLUE), (r2, GREEN), (-rb, MAROON)): axA.plot([0, x], [0, 0], color=col, lw=0.9, ls=":") axA.text(0.45, -0.55, r"$r_1$", color=BLUE, fontsize=13, bbox=BOX, zorder=7) axA.text(3.1, -0.62, r"$r_2$", color=GREEN, fontsize=13, bbox=BOX, zorder=7) axA.text(-2.1, 0.30, r"$r_bR$ (wins: $R>15.58$)") axB.axhline((np.sqrt(2) - 1)*(1 + 1/np.sqrt(R_CASE)), color=INK, lw=1.4, ls=(0, (2, 3)), zorder=3) # the three places where the curve touches or is quoted on the slide imax = int(np.argmax(three(R_CASE, lo))) for x, y, col in ((1.0, hR, MAROON), (R_CASE, hR, MAROON), (lo[imax], three(R_CASE, lo[imax]), BLUE), (RS_EX, three(R_CASE, RS_EX), BLUE), (RS_BI, three(R_CASE, RS_BI), PURPLE)): axB.plot(x, y, "o", color=col, ms=7, mfc=CREAM, mew=1.8, zorder=7) axB.annotate(f"worst stopover\n$R^*={lo[imax]:.2f}$, $+{100*(three(R_CASE, lo[imax])/hR - 1):.0f}\\%$", xy=(lo[imax], three(R_CASE, lo[imax])), xytext=(lo[imax] + 8.5, three(R_CASE, lo[imax]) - 0.008), arrowprops=dict(arrowstyle="-|>", color=BLUE, lw=1.2), color=BLUE, fontsize=10, ha="left", va="center") axB.annotate(f"the worked example\n$R^*=25$: $+{100*(three(R_CASE, RS_EX)/hR - 1):.2f}\\%$", xy=(RS_EX, three(R_CASE, RS_EX)), xytext=(34.0, three(R_CASE, RS_EX) + 0.030), arrowprops=dict(arrowstyle="-|>", color=BLUE, lw=1.2), color=BLUE, fontsize=10, ha="left", va="center") axB.annotate(f"$R^*=40$: ${100*(three(R_CASE, RS_BI)/hR - 1):.2f}\\%$", xy=(RS_BI, three(R_CASE, RS_BI)), xytext=(RS_BI + 15.0, three(R_CASE, RS_BI) - 0.014), arrowprops=dict(arrowstyle="-|>", color=PURPLE, lw=1.2), color=PURPLE, fontsize=10, ha="left", va="center") # the two ringed points on the Hohmann line are the degenerate stopovers axB.text(3.0, 0.5455, "both ends of the hump\nare the Hohmann transfer", color=MAROON, fontsize=9.5, ha="left", va="center", bbox=BOX, zorder=8) axB.text(108, hR - 0.004, "Hohmann", color=MAROON, fontsize=10.5, ha="right", va="top") axB.text(108, (np.sqrt(2) - 1)*(1 + 1/np.sqrt(R_CASE)) + 0.004, r"bi-parabolic limit, $r_b\to\infty$", color=INK, fontsize=9.5, ha="right", va="bottom", alpha=0.8) axB.axvline(R_CASE, color=INK, lw=1.0, ls=":", alpha=0.6) axB.text(R_CASE - 0.8, 0.4955, r"$R^{*}=R$", color=INK, fontsize=10, ha="right", va="bottom", alpha=0.85) axB.set_xlim(0, 112); axB.set_ylim(0.485, 0.735) axB.set_xlabel(r"detour radius $R^{*}=r_b/r_1$", fontsize=11.5) axB.set_ylabel(r"$\Delta v_{\rm total}/v_{c1}$", fontsize=11.5) axB.set_title(fr"(b) the same three burns, swept in $r_b$ (fixed $R={R_CASE:.0f}$)", fontsize=12) axB.legend(frameon=False, fontsize=10, loc="upper right", bbox_to_anchor=(1.0, 1.0)) axB.grid(True, color=INK, alpha=0.16, lw=0.7) despine(axB) fig.subplots_adjust(left=0.015, right=0.975, top=0.90, bottom=0.115, wspace=0.13) fig.savefig("img/fig-man-stopover.png", dpi=160, bbox_inches="tight") plt.close(fig) print("wrote img/fig-man-stopover.png") # -------------------------------------------------------------------------- # Verification # -------------------------------------------------------------------------- if __name__ == "__main__": print("\nboth ends of the stopover branch are the Hohmann transfer exactly") for R in (5.0, 15.0, R_CASE, 60.0): print(f" R={R:5.1f}: three(R,1)={three(R, 1.0):.9f} " f"three(R,R)={three(R, R):.9f} hohmann={hoh(R):.9f}") print("\nthe stopover never beats Hohmann, over a wide sweep") worst = 1.0 for R in np.concatenate([np.linspace(1.05, 5, 80), np.linspace(5, 200, 400)]): Rs = np.linspace(1.0, R, 4000) worst = min(worst, float((three(R, Rs) - hoh(R)).min())) print(f" min over (R, R*) of Dv_stop - Dv_H = {worst:+.3e} (>= 0)") print(f"\nthe drawn case R={r2/r1:.0f}, R*={rb/r1:.1f}") print(f" stopover {three(r2/r1, rb/r1):.6f} Hohmann {hoh(r2/r1):.6f} " f"penalty {pen:.2f}%") print(f"\nfixed R={R_CASE:.0f}, the sweep in R* quoted on the slide") for Rs in (1.0, 4.0214, 10.0, RS_EX, 29.9, R_CASE, RS_BI, 60.0, 110.0): tag = "stopover " if Rs < R_CASE else ("Hohmann " if Rs == R_CASE else "bi-ellip.") print(f" R*={Rs:7.2f} {tag}: {three(R_CASE, Rs):.6f} " f"({100*(three(R_CASE, Rs)/hR - 1):+7.2f}% vs Hohmann)") print(f" bi-parabolic limit : " f"{(np.sqrt(2)-1)*(1+1/np.sqrt(R_CASE)):.6f} " f"({100*((np.sqrt(2)-1)*(1+1/np.sqrt(R_CASE))/hR - 1):+7.2f}%)") print("\nsigning dv3 the bi-elliptic way understates the stopover by 2 dv3") d3 = abs(1/np.sqrt(R_CASE) - np.sqrt(2*RS_EX/(R_CASE*(R_CASE + RS_EX)))) signed = three(R_CASE, RS_EX) - 2*d3 print(f" R={R_CASE:.0f}, R*={RS_EX:.0f}: true {three(R_CASE, RS_EX):.6f} " f"signed formula {signed:.6f} 2*dv3 = {2*d3:.6f}")