"""Figures for patched-conics.qmd (img/fig-pc-*.png).

  fig-pc-geometry     the three bodies of the sphere-of-influence derivation
  fig-pc-ratios       Laplace's two perturbation ratios against distance from the Earth
  fig-pc-overview     the three conics of an Earth-Mars patched-conic transfer
  fig-pc-departure    the departure hyperbola from a 300 km parking orbit (Example 1)
  fig-pc-mars-arrival the arrival excess velocity at Mars, and the approach from in front (Example 1)
  fig-pc-moon         the lunar patch on the Moon's sphere of influence (Example 2)
  fig-pc-flyby        the Jupiter flyby and its velocity triangle (Example 3)
  fig-pc-dla          declination of a Hohmann departure asymptote over the year (Example 4)
  fig-pc-plane        out-of-plane arrival: tilted transfer against a broken-plane transfer (Example 5)
  fig-pc-arrival      the family of capture planes through an arrival asymptote (Example 6)

Run from the sfm-slides directory:   python3 codes/fig-pc-patched-conics.py
"""
import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import solve_ivp
from _manstyle import INK, MAROON, BLUE, GREEN, AMBER, PURPLE, despine

muS, AU = 1.32712440018e11, 1.495978707e8
mu, Re = 398600.4418, 6378.137
muM, RM, D = 4902.800066, 1737.4, 384400.0
muMa, RMa, aMa = 42828.37, 3396.2, 1.52368 * AU
muJ, RJ, aJ = 1.26686534e8, 71492.0, 5.2026 * AU
EPS = np.radians(23.4392911)
t = np.linspace(0, 2 * np.pi, 400)


def save(fig, name):
    fig.savefig(f"img/{name}.png", dpi=150, bbox_inches="tight")
    plt.close(fig)


def arrow(ax, p, q, c, lw=2.0, txt=None, off=(0, 0), fs=13):
    ax.annotate("", xy=q, xytext=p, arrowprops=dict(arrowstyle="-|>", color=c, lw=lw, shrinkA=0, shrinkB=0))
    if txt:
        m = (np.array(p) + np.array(q)) / 2 + np.array(off)
        ax.text(*m, txt, color=c, fontsize=fs, ha="center", va="center")


# ------------------------------------------------------------------ geometry of the derivation
fig, ax = plt.subplots(figsize=(8.4, 4.6))
S, Pl, s = np.array([0.0, 0.0]), np.array([7.0, 0.0]), np.array([8.3, 1.7])
ax.scatter(*S, s=700, color=AMBER, zorder=3)
ax.text(S[0], S[1] - 0.75, "Sun, $M$", ha="center", fontsize=12)
ax.scatter(*Pl, s=220, color=BLUE, zorder=3)
ax.text(Pl[0], Pl[1] - 0.6, "planet, $m$", ha="center", fontsize=12)
ax.scatter(*s, s=40, color=INK, zorder=3)
ax.text(s[0] + 0.15, s[1] + 0.15, "spacecraft", fontsize=12)
arrow(ax, S, Pl, BLUE, txt=r"$\mathbf{R}$", off=(0, -0.3))
arrow(ax, Pl, s, MAROON, txt=r"$\mathbf{r}$", off=(0.35, -0.1))
arrow(ax, S, s, GREEN, txt=r"$\boldsymbol{\rho}=\mathbf{R}+\mathbf{r}$", off=(-0.6, 0.45))
a = np.linspace(0, np.arctan2(s[1] - Pl[1], s[0] - Pl[0]), 30)
ax.plot(Pl[0] + 0.55 * np.cos(a), Pl[1] + 0.55 * np.sin(a), color=MAROON, lw=1.2)
ax.plot([Pl[0], Pl[0] + 1.6], [0, 0], color=INK, lw=0.7, ls=":")
ax.text(Pl[0] + 0.7, 0.18, r"$\theta$", color=MAROON, fontsize=13)
ax.set_aspect("equal")
ax.axis("off")
ax.set_xlim(-0.8, 10.2)
ax.set_ylim(-1.2, 2.5)
save(fig, "fig-pc-geometry")

# ------------------------------------------------------------------ the two ratios
mE = mu / muS
r = np.logspace(4.3, 6.6, 500)                    # km from the Earth
R = AU
fig, ax = plt.subplots(figsize=(8.6, 5.0))
for th, ls, lab in ((0.0, "-", r"along the Sun-Earth line ($\theta=0$)"), (np.pi / 2, "--", r"perpendicular ($\theta=90^\circ$)")):
    Qp = (r / R) ** 3 / mE * np.sqrt(1 + 3 * np.cos(th) ** 2)       # planet-centred: Sun's tide / Earth's pull
    Qs = mE * (R / r) ** 2                                          # Sun-centred: Earth's pull / Sun's pull
    ax.loglog(r, Qp, color=BLUE, ls=ls, lw=2, label=f"planet-centred ratio, {lab}")
    ax.loglog(r, Qs, color=MAROON, ls="-", lw=2, label="Sun-centred ratio (any direction)" if th == 0 else None)
rs = R * mE ** 0.4
for x, c, lab in ((R * np.sqrt(mE) / (1 + np.sqrt(mE)), INK, "equal pulls\n260 000 km"),
                  (rs, GREEN, "Laplace radius\n925 000 km"), (R * (mE / 3) ** (1 / 3), PURPLE, "Hill radius\n1.50 million km")):
    ax.axvline(x, color=c, lw=1.2, ls=":")
    ax.text(x * 1.04, 3e-6 if c != GREEN else 3e-4, lab, color=c, fontsize=9)
ax.set_xlabel("distance from the Earth, $r$ [km]")
ax.set_ylabel("perturbing / main acceleration")
ax.set_ylim(1e-6, 1e2)
ax.legend(frameon=False, fontsize=8.5, loc="upper left")
ax.set_title("Which centre makes the perturbation smaller? (Earth and Sun)", fontsize=12)
despine(ax)
save(fig, "fig-pc-ratios")

# ------------------------------------------------------------------ overview of a patched-conic transfer
fig = plt.figure(figsize=(14, 5.2))
gs = fig.add_gridspec(1, 3, width_ratios=[1, 1.6, 1])
a1, a2, a3 = [fig.add_subplot(gs[0, k]) for k in range(3)]
at = (AU + aMa) / 2
eT = (aMa - AU) / (aMa + AU)
f = np.linspace(0, np.pi, 300)
rT = at * (1 - eT ** 2) / (1 + eT * np.cos(f))
a2.plot(np.cos(t), np.sin(t), color=BLUE, lw=1, ls=(0, (4, 3)))
a2.plot(aMa / AU * np.cos(t), aMa / AU * np.sin(t), color=MAROON, lw=1, ls=(0, (4, 3)))
a2.plot(rT * np.cos(f) / AU, rT * np.sin(f) / AU, color=GREEN, lw=2.2)
a2.scatter([0], [0], s=150, color=AMBER, zorder=3)
a2.scatter([1], [0], s=50, color=BLUE, zorder=3)
a2.scatter([-aMa / AU], [0], s=50, color=MAROON, zorder=3)
for c0, rr, lab in ((np.array([1.0, 0]), 0.0062, "Earth's\nSOI"), (np.array([-aMa / AU, 0]), 0.0039, "Mars's\nSOI")):
    a2.add_patch(plt.Circle(c0, 0.09, fill=False, color=INK, lw=1, ls=":"))
    a2.text(c0[0], c0[1] - 0.2, lab + "\n(drawn 15x)", ha="center", va="top", fontsize=8.5)
a2.text(0, 1.1, "2. heliocentric ellipse\n(Sun only)", ha="center", fontsize=10, color=GREEN)
a2.set_aspect("equal")
a2.axis("off")
a2.set_xlim(-1.8, 1.4)
a2.set_ylim(-1.7, 1.7)
for ax, mup, Rp, vinf, rp, lab, col, sgn in ((a1, mu, Re, 2.9447, Re + 300, "1. departure hyperbola\n(Earth only)", BLUE, 1),
                                               (a3, muMa, RMa, 2.6489, RMa + 300, "3. arrival hyperbola\n(Mars only)", MAROON, -1)):
    e = 1 + rp * vinf ** 2 / mup
    p = rp * (1 + e)
    ff = np.linspace(-np.arccos(-1 / e) + 0.05, np.arccos(-1 / e) - 0.05, 400)
    rh = p / (1 + e * np.cos(ff))
    ok = rh < 12 * Rp
    ax.plot((rh * np.cos(ff))[ok] / Rp, (rh * np.sin(ff))[ok] / Rp, color=col, lw=2)
    ax.add_patch(plt.Circle((0, 0), 1, color=col, alpha=0.5))
    ax.plot(rp / Rp * np.cos(t), rp / Rp * np.sin(t), color=INK, lw=0.8, ls=":")
    ax.set_aspect("equal")
    ax.axis("off")
    ax.set_xlim(-6, 6)
    ax.set_ylim(-7, 7)
    ax.text(0, 6.2, lab, ha="center", fontsize=10, color=col)
    ax.text(0, -6.6, "planet radii", ha="center", fontsize=8.5)
save(fig, "fig-pc-overview")

# ------------------------------------------------------------------ departure hyperbola, Example 1
vinf = 2.9446935
rp = Re + 300
e = 1 + rp * vinf ** 2 / mu
finf = np.arccos(-1 / e)
w = np.pi / 2 - finf                      # periapsis direction so that the asymptote points along +y
p = rp * (1 + e)
ff = np.linspace(-0.3, finf - 0.02, 600)
rh = p / (1 + e * np.cos(ff))
ok = rh < 9 * Re
fig, ax = plt.subplots(figsize=(7.4, 7.2))
ax.add_patch(plt.Circle((0, 0), 1, color=BLUE, alpha=0.5))
ax.plot(rp / Re * np.cos(t), rp / Re * np.sin(t), color=INK, lw=1, ls=":", label="300 km parking orbit")
ax.plot((rh * np.cos(ff + w))[ok] / Re, (rh * np.sin(ff + w))[ok] / Re, color=GREEN, lw=2.2, label="departure hyperbola")
c = np.array([np.cos(w), np.sin(w)]) * (p / (e * e - 1)) * e / Re          # centre: |a| e from the focus, beyond periapsis
ax.plot([c[0], c[0]], [c[1], c[1] + 30], color=GREEN, lw=1, ls="--", label="asymptote, parallel to the Earth's velocity")
pp = rp / Re * np.array([np.cos(w), np.sin(w)])
ax.scatter(*pp, s=60, color=MAROON, zorder=4)
ax.text(pp[0] + 0.25, pp[1] - 0.35, "burn here\n(periapsis)", color=MAROON, fontsize=10)
aa = np.linspace(w, w + finf, 60)
ax.plot(0.55 * np.cos(aa), 0.55 * np.sin(aa), color=MAROON, lw=1.4)
ax.text(-2.9, -1.3, r"$f_\infty=150.8^\circ$", color=MAROON, fontsize=11)
arrow(ax, (3.6, -3.6), (3.6, -1.6), INK, txt=r"Earth's velocity $\mathbf{V}_\oplus$", off=(0.0, -1.35), fs=10)
arrow(ax, (-1.6, -4.4), (-3.6, -4.4), AMBER, txt="to the Sun", off=(0, 0.35), fs=10)
ax.set_aspect("equal")
ax.set_xlim(-5, 5)
ax.set_ylim(-5, 9)
ax.set_xlabel("Earth radii")
ax.legend(frameon=False, fontsize=9, loc="upper right")
ax.set_title("Example 1: leaving the Earth for Mars", fontsize=12)
save(fig, "fig-pc-departure")

# ------------------------------------------------------------------ arrival at Mars, Example 1
at = (AU + aMa) / 2
VM = np.sqrt(muS / aMa)
Vt2 = np.sqrt(muS * (2 / aMa - 1 / at))
vi2 = VM - Vt2
rp2 = RMa + 300
e2 = 1 + rp2 * vi2 ** 2 / muMa
fig, (aL, aR) = plt.subplots(1, 2, figsize=(14, 5.0), gridspec_kw=dict(width_ratios=[1.0, 1.15]))
# left: the subtraction at the patch, all three vectors along Mars's direction of motion
arrow(aL, (0, 2.0), (VM, 2.0), MAROON, lw=2.4, txt=r"$\mathbf{V}_{\rm M}=24.129$ (Mars)", off=(0, 0.35), fs=12)
arrow(aL, (0, 0.0), (Vt2, 0.0), GREEN, lw=2.4, txt=r"$\mathbf{V}_{t2}=21.480$ (spacecraft, at aphelion)", off=(0, -0.35), fs=12)
arrow(aL, (VM, 1.0), (Vt2, 1.0), BLUE, lw=2.4)
aL.text((VM + Vt2) / 2, 1.38, r"$\mathbf{v}_{\infty2}=\mathbf{V}_{t2}-\mathbf{V}_{\rm M}$", color=BLUE, fontsize=12, ha="center")
aL.text((VM + Vt2) / 2, 0.62, "2.649, backwards", color=BLUE, fontsize=11, ha="center")
for x in (Vt2, VM):
    aL.plot([x, x], [-0.2, 2.2], color=INK, lw=0.8, ls=":")
aL.set_xlim(-0.5, 27)
aL.set_ylim(-0.7, 2.7)
aL.set_yticks([])
aL.set_xlabel("speed along Mars's direction of motion [km/s]")
aL.set_title("At the patch: the spacecraft is slower than Mars", fontsize=12)
despine(aL)
aL.spines["left"].set_visible(False)
# right: Mars-centred view; Mars moves towards +x, the Sun is towards +y
fi = np.arccos(-1 / e2)
w = fi                                    # incoming asymptote direction (where it comes from) is +x
pp2 = rp2 * (1 + e2)
ff = np.linspace(-fi + 0.02, 0.0, 600)
rh = pp2 / (1 + e2 * np.cos(ff))
ok = rh < 14 * RMa
aR.plot((rh * np.cos(ff + w))[ok] / RMa, (rh * np.sin(ff + w))[ok] / RMa, color=GREEN, lw=2.2, label="arrival hyperbola")
ff2 = np.linspace(0.0, fi - 0.02, 600)
rh2 = pp2 / (1 + e2 * np.cos(ff2))
ok2 = rh2 < 14 * RMa
aR.plot((rh2 * np.cos(ff2 + w))[ok2] / RMa, (rh2 * np.sin(ff2 + w))[ok2] / RMa, color=GREEN, lw=1, ls="--",
        label="(if no burn)")
aR.plot(rp2 / RMa * np.cos(t), rp2 / RMa * np.sin(t), color=INK, lw=1, ls=":", label="300 km capture orbit")
aR.add_patch(plt.Circle((0, 0), 1, color=MAROON, alpha=0.5))
bA = rp2 * np.sqrt(1 + 2 * muMa / (rp2 * vi2 ** 2)) / RMa
aR.plot([0, 14], [bA, bA], color=BLUE, lw=0.8, ls=(0, (2, 2)))
aR.plot([0, 14], [0, 0], color=BLUE, lw=0.8, ls=(0, (2, 2)))
aR.annotate("", xy=(7.5, bA), xytext=(7.5, 0), arrowprops=dict(arrowstyle="<->", color=BLUE, lw=1, shrinkA=0, shrinkB=0))
aR.text(7.8, bA / 2, r"$b=7667$ km", color=BLUE, fontsize=10, va="center")
arrow(aR, (12.5, bA + 0.9), (9.0, bA + 0.9), BLUE, lw=2.0, txt=r"$\mathbf{v}_{\infty2}$", off=(0, 0.6), fs=12)
pk = rp2 / RMa * np.array([np.cos(w), np.sin(w)])
aR.scatter(*pk, s=50, color=MAROON, zorder=4)
aR.text(pk[0] - 0.3, pk[1] + 0.5, "capture burn", color=MAROON, fontsize=10, ha="right")
arrow(aR, (-12.5, -6.0), (-8.5, -6.0), INK, txt=r"Mars's velocity $\mathbf{V}_{\rm M}$", off=(0, -0.8), fs=10)
arrow(aR, (8.0, 5.0), (8.0, 7.5), AMBER, txt="to the Sun", off=(2.0, 0), fs=10)
aR.axvline(0, color=INK, lw=0.6, ls="-.", alpha=0.5)
aR.text(7.0, -6.6, "ahead of Mars\n(leading side)", fontsize=10, ha="center")
aR.text(-7.0, -2.6, "behind Mars\n(trailing side)", fontsize=10, ha="center")
aR.set_aspect("equal")
aR.set_xlim(-14, 14)
aR.set_ylim(-8, 8)
aR.set_xlabel("Mars radii")
aR.legend(frameon=False, fontsize=9, loc="upper left")
aR.set_title("Seen from Mars: the spacecraft comes in from in front", fontsize=12)
save(fig, "fig-pc-mars-arrival")

# ------------------------------------------------------------------ lunar patch, Example 2
Rs = D * (muM / mu) ** 0.4
vM = np.sqrt((mu + muM) / D)
nM = vM / D
r0, v0, lam = Re + 300, 10.84, np.radians(60.0)
h0 = r0 * v0
pp0 = h0 ** 2 / mu
e0 = pp0 / r0 - 1
r1 = np.sqrt(D ** 2 + Rs ** 2 - 2 * D * Rs * np.cos(lam))
gam = np.arcsin(Rs / r1 * np.sin(lam))
f1 = np.arccos((pp0 / r1 - 1) / e0)
thM = f1 - gam
a0 = pp0 / (1 - e0 ** 2)
E1 = 2 * np.arctan(np.sqrt((1 - e0) / (1 + e0)) * np.tan(f1 / 2))
t1 = (E1 - e0 * np.sin(E1)) * np.sqrt(a0 ** 3 / mu)
M0 = thM - nM * t1


def rhs(tt, y):
    rr = y[:2]
    th = M0 + nM * tt
    rm = D * np.array([np.cos(th), np.sin(th)])
    acc = -mu * rr / np.linalg.norm(rr) ** 3 - muM * ((rr - rm) / np.linalg.norm(rr - rm) ** 3 + rm / D ** 3)
    return [y[2], y[3], acc[0], acc[1]]


num = solve_ivp(rhs, [0, 4.2 * 86400], [r0, 0, 0, v0], rtol=1e-11, atol=1e-8, dense_output=True, max_step=600)
fig, (aL, aR) = plt.subplots(1, 2, figsize=(14, 6.2))
ff = np.linspace(0, f1, 300)
rg = pp0 / (1 + e0 * np.cos(ff))
aL.plot(rg * np.cos(ff) / 1e3, rg * np.sin(ff) / 1e3, color=GREEN, lw=2, label="geocentric ellipse (Earth only)")
aL.plot(D * np.cos(t) / 1e3, D * np.sin(t) / 1e3, color=INK, lw=0.8, ls=(0, (4, 3)))
Mp = D * np.array([np.cos(thM), np.sin(thM)])
aL.add_patch(plt.Circle(Mp / 1e3, Rs / 1e3, fill=False, color=MAROON, lw=1.2, ls="--"))
aL.scatter(*Mp / 1e3, color="#9a9388", s=40, zorder=4)
P1 = r1 * np.array([np.cos(f1), np.sin(f1)])
aL.scatter(*P1 / 1e3, color=MAROON, s=40, zorder=4)
aL.plot([0, Mp[0] / 1e3], [0, Mp[1] / 1e3], color=INK, lw=0.7, ls=":")
aL.plot([0, P1[0] / 1e3], [0, P1[1] / 1e3], color=INK, lw=0.7, ls=":")
aL.plot([Mp[0] / 1e3, P1[0] / 1e3], [Mp[1] / 1e3, P1[1] / 1e3], color=MAROON, lw=1)
aL.text(P1[0] / 1e3 - 40, P1[1] / 1e3 + 18, "patch point", color=MAROON, fontsize=10)
aL.text(Mp[0] / 1e3 - 10, Mp[1] / 1e3 - 30, "Moon", fontsize=10)
aL.text(Mp[0] / 1e3 + 50, Mp[1] / 1e3 + 40, "Moon's SOI\n66 200 km", color=MAROON, fontsize=9)
aL.scatter([0], [0], s=80, color=BLUE, zorder=4)
aL.text(10, -25, "Earth", fontsize=10, color=BLUE)
aL.set_aspect("equal")
aL.set_xlabel("[1000 km]")
aL.legend(frameon=False, fontsize=9, loc="lower left")
aL.set_title(r"Example 2: the patch at $\lambda_1=60^\circ$ on the Moon's SOI", fontsize=12)
aL.set_xlim(-420, 120)
aL.set_ylim(-100, 420)
# right: Moon-centred view, patched conic against numerical
ts = np.linspace(0, 4.2 * 86400, 3000)
Y = num.sol(ts)
th = M0 + nM * ts
Xm = np.cos(th) * (Y[0] - D * np.cos(th)) + np.sin(th) * (Y[1] - D * np.sin(th))
Ym = -np.sin(th) * (Y[0] - D * np.cos(th)) + np.cos(th) * (Y[1] - D * np.sin(th))
near = np.hypot(Xm, Ym) < 1.3 * Rs
aR.plot(Xm[near] / 1e3, Ym[near] / 1e3, color=INK, lw=1.6, label="numerical (Earth and Moon), perilune 1030 km")
rho = P1 - Mp
VMv = vM * np.array([-np.sin(thM), np.cos(thM)])
V1 = mu / h0 * np.array([-np.sin(f1), e0 + np.cos(f1)])
v2 = V1 - VMv
sel = solve_ivp(lambda tt, y: [y[2], y[3], *(-muM * y[:2] / np.linalg.norm(y[:2]) ** 3)], [0, 1.6e5], [*rho, *v2],
                rtol=1e-11, atol=1e-8, dense_output=True)
Z = sel.sol(np.linspace(0, 1.6e5, 2000))
Rot = np.array([[np.cos(thM), np.sin(thM)], [-np.sin(thM), np.cos(thM)]])       # into the Earth-Moon frame at patch time
Zr = Rot @ Z[:2]
ok = np.hypot(*Zr) < 1.3 * Rs
aR.plot(Zr[0][ok] / 1e3, Zr[1][ok] / 1e3, color=MAROON, lw=1.8, ls="--", label="patched conic (Moon only), perilune 1627 km")
aR.add_patch(plt.Circle((0, 0), Rs / 1e3, fill=False, color=MAROON, lw=1, ls=":"))
aR.add_patch(plt.Circle((0, 0), RM / 1e3, color="#9a9388", alpha=0.7))
aR.arrow(-75, -70, -15, 0, head_width=3, color=BLUE)
aR.text(-75, -64, "to the Earth", fontsize=9, color=BLUE)
aR.set_aspect("equal")
aR.set_xlim(-90, 90)
aR.set_ylim(-80, 80)
aR.set_xlabel("[1000 km], Moon-centred, $x$ away from the Earth")
aR.legend(frameon=False, fontsize=9, loc="upper left")
aR.set_title("Inside the Moon's SOI (frame turning with the Moon)", fontsize=12)
save(fig, "fig-pc-moon")

# ------------------------------------------------------------------ Jupiter flyby, Example 3
VJ = np.sqrt(muS / aJ)
atJ = (AU + aJ) / 2
Va = np.sqrt(muS * (2 / aJ - 1 / atJ))
vi = VJ - Va
rpJ = 5 * RJ
eJ = 1 + rpJ * vi ** 2 / muJ
dl = 2 * np.arcsin(1 / eJ)
fig, (aL, aR) = plt.subplots(1, 2, figsize=(14, 6.2))
bJ = rpJ * np.sqrt(1 + 2 * muJ / (rpJ * vi ** 2))           # impact parameter
best = None
for side in (1, -1):                                          # pass on either side; keep the one turning towards +x
    tr = solve_ivp(lambda tt, y: [y[2], y[3], *(-muJ * y[:2] / np.linalg.norm(y[:2]) ** 3)], [0, 6e6],
                   [side * bJ, 60 * RJ, 0.0, -vi], rtol=1e-11, atol=1e-6, dense_output=True, max_step=2000)
    if tr.y[2, -1] > 0:
        best = tr
Zj = best.sol(np.linspace(0, best.t[-1], 4000))
ok = np.hypot(Zj[0], Zj[1]) < 40 * RJ
aL.plot(Zj[0][ok] / RJ, Zj[1][ok] / RJ, color=GREEN, lw=2.2)
aL.annotate("", xy=(Zj[0][ok][40] / RJ, Zj[1][ok][40] / RJ), xytext=(Zj[0][ok][0] / RJ, Zj[1][ok][0] / RJ),
            arrowprops=dict(arrowstyle="-|>", color=GREEN))
aL.add_patch(plt.Circle((0, 0), 1, color=AMBER, alpha=0.6))
arrow(aL, (14, -14), (14, -6), INK, txt=r"$\mathbf{V}_J$", off=(1.6, 0), fs=12)
aL.text(-15, 15, "passes behind Jupiter:\nthe trailing side", fontsize=10, color=GREEN)
aL.set_aspect("equal")
aL.set_xlim(-20, 20)
aL.set_ylim(-20, 20)
aL.set_xlabel("Jupiter radii")
aL.set_title(r"Example 3: flyby at $r_p=5R_J$ (Jupiter-centred)", fontsize=12)
O = np.array([0.0, 0.0])
Vj = np.array([0.0, VJ])
vin = np.array([0.0, -vi])
vout = vi * np.array([np.sin(dl), -np.cos(dl)])
arrow(aR, O, Vj, INK, txt=r"$\mathbf{V}_J=13.06$", off=(-1.6, 0), fs=11)
arrow(aR, Vj, Vj + vin, BLUE, txt=r"$\mathbf{v}_{\infty}^-$", off=(-0.9, 0), fs=12)
arrow(aR, Vj, Vj + vout, MAROON, txt=r"$\mathbf{v}_{\infty}^+$", off=(0.4, 0.6), fs=12)
arrow(aR, O, Vj + vin, BLUE, lw=1.2, txt=r"$\mathbf{V}^-=7.41$", off=(1.4, 0), fs=11)
arrow(aR, O, Vj + vout, MAROON, lw=1.2, txt=r"$\mathbf{V}^+=17.41$", off=(1.6, -0.5), fs=11)
aa = np.linspace(-np.pi / 2, -np.pi / 2 + dl, 60)
aR.plot(Vj[0] + 1.8 * np.cos(aa), Vj[1] + 1.8 * np.sin(aa), color=GREEN, lw=1.4)
aR.text(Vj[0] + 2.0, Vj[1] - 1.0, r"$\delta=133.1^\circ$", color=GREEN, fontsize=11)
aR.set_aspect("equal")
aR.set_xlim(-6, 9)
aR.set_ylim(-1, 19)
aR.set_xlabel("km/s")
aR.set_title("Velocity triangle (heliocentric)", fontsize=12)
despine(aR)
save(fig, "fig-pc-flyby")

# ------------------------------------------------------------------ declination of the departure asymptote, Example 4
ls = np.linspace(0, 360, 361)
dla = np.degrees(np.arcsin(np.sin(EPS) * np.sin(np.radians(ls - 90))))
fig, ax = plt.subplots(figsize=(9, 4.8))
ax.plot(ls, dla, color=BLUE, lw=2, label="Hohmann departure to Mars: asymptote along the Earth's velocity")
ax.axhspan(-13.72, 13.72, color=GREEN, alpha=0.12, label="reachable from a due-east parking orbit at Sriharikota ($i=13.7^\\circ$)")
for x, y, lab in ((0, -23.44, "near 20 Mar"), (248, 8.57, "near 30 Nov")):
    ax.scatter([x], [y], color=MAROON, zorder=4)
    ax.text(x + 6, y + (-2.5 if y < 0 else 1.5), lab, color=MAROON, fontsize=10)
ax.set_xlabel(r"the Sun's ecliptic longitude on the departure day, $\lambda_\odot$ [deg]")
ax.set_ylabel(r"declination of the asymptote, $\delta_\infty$ [deg]")
ax.set_xticks(range(0, 361, 45))
ax.set_ylim(-30, 30)
ax.legend(frameon=False, fontsize=9, loc="upper left")
ax.set_title("Example 4: where in the sky the departure asymptote points", fontsize=12)
despine(ax)
save(fig, "fig-pc-dla")

# ------------------------------------------------------------------ out-of-plane arrival, Example 5
at = (AU + aMa) / 2
Vt = np.sqrt(muS * (2 / AU - 1 / at))
VE = np.sqrt(muS / AU)
beta = np.radians(1.85)
thd = np.linspace(120, 176, 300)
it = np.arcsin(np.sin(beta) / np.sin(np.radians(thd)))
vinf_t = np.sqrt(Vt ** 2 + VE ** 2 - 2 * Vt * VE * np.cos(it))
dv_leo = np.sqrt(vinf_t ** 2 + 2 * mu / (Re + 300)) - np.sqrt(mu / (Re + 300))
fig, ax = plt.subplots(figsize=(9, 4.8))
ax.plot(thd, dv_leo, color=BLUE, lw=2, label="tilt the whole transfer at departure")
ax.axhline(3.590, color=INK, lw=1, ls=":", label="coplanar Hohmann, 3.590 km/s")
ax.axhline(3.590 + 0.875, color=MAROON, lw=1.6, ls="--", label="broken plane: Hohmann + 0.875 km/s at 90 deg before arrival")
ax.set_xlabel("transfer angle [deg]")
ax.set_ylabel(r"$\Delta v$ from a 300 km orbit [km/s]")
ax.set_ylim(3.4, 7.0)
ax.legend(frameon=False, fontsize=9, loc="upper left")
ax.set_title(r"Example 5: reaching Mars $1.85^\circ$ out of the ecliptic", fontsize=12)
despine(ax)
save(fig, "fig-pc-plane")

# ------------------------------------------------------------------ the capture planes, Example 6
dinf = np.radians(30.25)
thb = np.linspace(0, 360, 361)
inc = np.degrees(np.arccos(np.cos(dinf) * np.sin(np.radians(thb))))
fig, ax = plt.subplots(figsize=(9, 4.8))
ax.plot(thb, inc, color=BLUE, lw=2)
ax.axhline(150.0, color=MAROON, lw=1.2, ls="--", label="MOM's Mars orbit, $i=150.0^\\circ$")
ax.axhline(180 - 30.25, color=GREEN, lw=1.0, ls=":", label=r"$180^\circ-|\delta_\infty|=149.75^\circ$")
ax.axhline(30.25, color=GREEN, lw=1.0, ls=":")
ax.set_xlabel(r"angle $\theta$ of the orbit's normal about the arrival asymptote [deg]")
ax.set_ylabel("inclination to Mars's equator [deg]")
ax.set_xticks(range(0, 361, 45))
ax.set_ylim(0, 180)
ax.legend(frameon=False, fontsize=9, loc="center right")
ax.set_title(r"Example 6: every plane through an asymptote with $\delta_\infty=-30.25^\circ$", fontsize=12)
despine(ax)
save(fig, "fig-pc-arrival")
print("wrote nine figures")