"""Elliptical phasing: same period change, two apsides, two very different bills.

A spacecraft already on an ellipse (a, e) must arrive dt late.  Kepler III fixes
the phasing semi-major axis,

    a_ph = a (1 + dt/T)^(2/3),

but *not* where the burn is made.  A prograde impulse at periapsis leaves r_p
alone and pushes apoapsis out; a prograde impulse at apoapsis leaves r_a alone
and pushes periapsis up.  Both ellipses have the same a_ph, hence the same
period and the same delay -- they differ only in cost.

Differentiating vis-viva at fixed r gives dv = mu da / (2 v a^2): for a
prescribed da the impulse is inversely proportional to the speed at the burn
point, so periapsis -- the fastest point -- is the cheapest.  The ratio of the
two apsidal costs is v_a/v_p = (1-e)/(1+e).

Numbers are Example 7 of the deck: a = 10000 km, e = 0.3, dt = 600 s.

Run from the sfm-slides directory:

    python3 codes/fig-man-ell-phasing.py
"""
import numpy as np
import matplotlib.pyplot as plt
from _manstyle import CREAM, INK, MAROON, BLUE, GREEN, AMBER, PURPLE, despine

mu, RE = 398600.4418, 6378.137

A, ECC, DT = 10000.0, 0.3, 600.0
RP, RA = A*(1 - ECC), A*(1 + ECC)
T = 2*np.pi*np.sqrt(A**3/mu)
APH = A*(1 + DT/T)**(2/3)
DA = APH - A


def speed(r, a):
    return np.sqrt(mu*(2.0/r - 1.0/a))


def total_dv(r):
    """Both impulses: out onto the phasing ellipse at r, and back again."""
    return 2.0*abs(speed(r, APH) - speed(r, A))


def ellipse(rp, ra, n=800):
    """Conic with the given apsides, periapsis on +x, focus at the origin."""
    a = 0.5*(rp + ra)
    e = (ra - rp)/(ra + rp)
    f = np.linspace(0.0, 2*np.pi, n)
    r = a*(1 - e**2)/(1 + e*np.cos(f))
    return r*np.cos(f), r*np.sin(f)


# the two phasing ellipses: same a_ph, one apsis held fixed by each burn
RA_P = 2*APH - RP                     # periapsis burn: r_p held, apoapsis rises
RP_A = 2*APH - RA                     # apoapsis burn: r_a held, periapsis rises
DV_P, DV_A = total_dv(RP), total_dv(RA)

fig, (ax, axc) = plt.subplots(1, 2, figsize=(13.8, 5.6),
                              gridspec_kw=dict(width_ratios=[1.22, 1.0]))

# --------------------------------------------------------------------------
#  Left: the two candidate phasing ellipses (geometry exaggerated)
# --------------------------------------------------------------------------
# Example 7's own numbers (e = 0.3, da/a = 4%) draw as three indistinguishable
# near-circles, so the *shapes* here use a fatter ellipse and a bigger period
# change.  Every quoted number stays with the right-hand panel, where it is
# Example 7 exactly.
AG, EG, AGP = 1.0, 0.45, 1.18
RPG, RAG = AG*(1 - EG), AG*(1 + EG)
RAG_P = 2*AGP - RPG                   # periapsis burn: r_p held, apoapsis rises
RPG_A = 2*AGP - RAG                   # apoapsis burn: r_a held, periapsis rises

curves = (((RPG, RAG), BLUE, 2.4, "-", 4, r"original ellipse,  $a$"),
          ((RPG, RAG_P), MAROON, 2.6, "-", 5,
           r"phasing ellipse, burn at periapsis"),
          ((RPG_A, RAG), PURPLE, 2.0, (0, (6, 4)), 3,
           r"phasing ellipse, burn at apoapsis"))
for (rp, ra), col, lw, ls, z, lab in curves:
    x, y = ellipse(rp, ra)
    ax.plot(x, y, color=col, lw=lw, ls=ls, zorder=z, label=lab)

ax.plot(0, 0, "o", color=INK, ms=7, zorder=8)
ax.plot([-RAG_P, RPG_A], [0, 0], color=INK, lw=0.8, ls=":", alpha=0.5, zorder=1)

# the two burn points and their impulses (prograde; motion counter-clockwise)
L = 0.42
ax.plot(RPG, 0, "o", color=MAROON, ms=10, zorder=9)
ax.annotate("", xy=(RPG, L), xytext=(RPG, 0.03),
            arrowprops=dict(arrowstyle="-|>", color=AMBER, lw=2.6))
ax.text(RPG + 0.28, 0.74, "prograde at periapsis:\n"
        r"$r_p$ held, apoapsis out" "\n"
        r"$v$ largest $\Rightarrow$ cheapest",
        color=MAROON, fontsize=10, ha="left", va="center")

ax.plot(-RAG, 0, "o", color=PURPLE, ms=10, zorder=9)
ax.annotate("", xy=(-RAG, -L), xytext=(-RAG, -0.03),
            arrowprops=dict(arrowstyle="-|>", color=AMBER, lw=2.6))
ax.text(-RAG - 0.12, -0.95, "prograde at apoapsis:\n"
        r"$r_a$ held, periapsis lifted" "\n"
        r"$v$ smallest $\Rightarrow$ dearest",
        color=PURPLE, fontsize=10, ha="right", va="center")

ax.text(0.0, -1.30,
        r"both phasing ellipses share the same $a_{ph}=a\,(1+\Delta t/T)^{2/3}$,"
        "\n" r"so both return $\Delta t$ late --- they differ only in the bill",
        color=INK, fontsize=11, ha="center", va="top")
ax.text(-3.06, 1.44, "shapes exaggerated", color=INK, fontsize=8.5, alpha=0.55,
        ha="left", va="top")

ax.legend(frameon=False, fontsize=9.5, loc="upper right",
          bbox_to_anchor=(1.03, 1.04), handlelength=2.6)
ax.set_aspect("equal"); ax.axis("off")
ax.set_xlim(-3.08, 2.42)
ax.set_ylim(-1.78, 1.48)
ax.set_title("Same period change, two places to buy it", fontsize=12.5)

# --------------------------------------------------------------------------
#  Right: the cost as a function of where the burn is made
# --------------------------------------------------------------------------
f = np.linspace(0.0, 2*np.pi, 721)
r = A*(1 - ECC**2)/(1 + ECC*np.cos(f))
v = speed(r, A)
exact = 2.0*(speed(r, APH) - v)*1000.0          # m/s
approx = mu*DA/(v*A**2)*1000.0                  # m/s, from dv = mu da / (2 v a^2)
fd = np.degrees(f)

axc.plot(fd, exact, color=MAROON, lw=2.6, zorder=4, label=r"exact  $2|v_{ph}-v|$")
axc.plot(fd, approx, color=MAROON, lw=1.4, ls=(0, (5, 4)), alpha=0.8, zorder=3,
         label=r"$\mu\,\Delta a/(v\,a^{2})$  (infinitesimal law)")
axc.plot([0.0, 360.0], [1000*DV_P]*2, color=INK, lw=0.8, ls=":", alpha=0.5)
axc.plot([180.0], [1000*DV_A], "o", color=PURPLE, ms=9, zorder=5)
axc.plot([0.0], [1000*DV_P], "o", color=MAROON, ms=9, zorder=5)
axc.plot([360.0], [1000*DV_P], "o", color=MAROON, ms=9, zorder=5)

axc.annotate(f"periapsis: {1000*DV_P:.1f} m/s", xy=(0.0, 1000*DV_P),
             xytext=(18.0, 1000*DV_P - 34.0), color=MAROON, fontsize=10.5,
             ha="left", va="top",
             arrowprops=dict(arrowstyle="-", color=MAROON, lw=0.9, alpha=0.7))
axc.annotate(f"apoapsis: {1000*DV_A:.1f} m/s", xy=(180.0, 1000*DV_A),
             xytext=(180.0, 1000*DV_A + 26.0), color=PURPLE, fontsize=10.5,
             ha="center", va="bottom",
             arrowprops=dict(arrowstyle="-", color=PURPLE, lw=0.9, alpha=0.7))

axc.annotate("", xy=(316.0, 1000*DV_A), xytext=(316.0, 1000*DV_P),
             arrowprops=dict(arrowstyle="<|-|>", color=INK, lw=1.4, alpha=0.8))
axc.text(311.0, 500*(DV_P + DV_A), r"$\dfrac{1+e}{1-e}=$"
         f"{(1 + ECC)/(1 - ECC):.2f}", color=INK, fontsize=12.5, ha="right",
         va="center",
         bbox=dict(facecolor=CREAM, edgecolor="none", pad=1.5))

# the speed itself, to show that the cost curve is simply its mirror image
axv = axc.twinx()
axv.plot(fd, v, color=GREEN, lw=1.8, alpha=0.75)
axv.set_ylabel("speed at the burn point  $v$   [km/s]", color=GREEN)
axv.tick_params(axis="y", colors=GREEN)
axv.spines["top"].set_visible(False)
axv.spines["right"].set_color(GREEN)
axv.set_ylim(3.6, 9.6)
axv.text(96.0, speed(A*(1 - ECC**2)/(1 + ECC*np.cos(np.radians(96))), A) + 0.22,
         r"$v(f)$", color=GREEN, fontsize=13, ha="center", va="bottom")

axc.set_xlabel(r"true anomaly of the burn point  $f$   [deg]")
axc.set_ylabel(r"total $\Delta v$   [m/s]", color=MAROON)
axc.tick_params(axis="y", colors=MAROON)
axc.set_xlim(0, 360)
axc.set_xticks(np.arange(0, 361, 60))
axc.set_ylim(120, 400)
axc.set_zorder(axv.get_zorder() + 1)
axc.patch.set_visible(False)
axc.legend(frameon=False, fontsize=10, loc="upper left")
axc.set_title(r"Cost $\propto 1/v$: burn where you are fastest", fontsize=12.5)
despine(axc)
axc.spines["left"].set_color(MAROON)

fig.subplots_adjust(left=0.02, right=0.93, top=0.90, bottom=0.12, wspace=0.06)
fig.savefig("img/fig-man-ell-phasing.png", dpi=160, bbox_inches="tight")
plt.close(fig)
print("wrote img/fig-man-ell-phasing.png")

# --------------------------------------------------------------------------
#  Verification against Example 7 of the deck
# --------------------------------------------------------------------------
if __name__ == "__main__":
    print(f"\nT      = {T:.5f} s   (deck: 9952.01)")
    print(f"a_ph   = {APH:.5f} km (deck: 10397.99)")
    for nm, rr in (("periapsis", RP), ("apoapsis", RA)):
        print(f"  {nm:>9}: v = {speed(rr, A):.6f} -> {speed(rr, APH):.6f}, "
              f"dv each = {speed(rr, APH) - speed(rr, A):.6f}, "
              f"total = {total_dv(rr):.6f} km/s")
    print(f"  cost ratio apo/peri = {DV_A/DV_P:.4f}  vs  "
          f"(1+e)/(1-e) = {(1 + ECC)/(1 - ECC):.4f}")
    print(f"  min of the exact curve at f = {fd[np.argmin(exact)]:.1f} deg, "
          f"max at f = {fd[np.argmax(exact)]:.1f} deg")