"""Flight path angle across an elliptic orbit.

    tan(gamma) = v_r / v_perp = e sin f / (1 + e cos f)

Everything about the size of the orbit -- a, and mu with it -- cancels in that
ratio, so gamma is a function of (e, f) alone.  Its extrema sit at cos f = -e,
i.e. at r = a (the ends of the minor axis), where sin(gamma) = e.
"""

import os
import numpy as np
import matplotlib.pyplot as plt

os.makedirs("img", exist_ok=True)

CREAM, INK = "#fcfbf7", "#2f2a24"
MAROON, BLUE, GREEN, AMBER, PURPLE = "#6b1f1f", "#2f5ea8", "#2d6a2d", "#b45309", "#7e22ce"

plt.rcParams.update({
    "figure.facecolor": CREAM, "savefig.facecolor": CREAM, "axes.facecolor": CREAM,
    "axes.edgecolor": INK, "axes.labelcolor": INK, "text.color": INK,
    "xtick.color": INK, "ytick.color": INK, "font.size": 11,
    "axes.grid": True, "grid.alpha": 0.22, "grid.color": INK,
})

f = np.linspace(0.0, 360.0, 2001)
fr = np.deg2rad(f)

fig, (axE, axA, axM) = plt.subplots(1, 3, figsize=(15.0, 4.9))

# ── panel 1: the family in e ────────────────────────────────────────────────
ecc = [0.0, 0.2, 0.4, 0.6, 0.8]
colors = [INK, BLUE, GREEN, AMBER, PURPLE]
locus_f, locus_g = [], []

for e, col in zip(ecc, colors):
    g = np.rad2deg(np.arctan2(e*np.sin(fr), 1.0 + e*np.cos(fr)))
    axE.plot(f, g, color=col, lw=2.2 if e else 1.5,
             ls="-" if e else (0, (4, 3)), label=f"$e = {e:g}$")
    if e > 0.0:
        fm = np.rad2deg(np.arccos(-e))
        gm = np.rad2deg(np.arcsin(e))
        axE.plot([fm, 360.0 - fm], [gm, -gm], "o", color=col, ms=6, zorder=6)
        locus_f.append(fm); locus_g.append(gm)

axE.plot(locus_f, locus_g, color=INK, lw=1.0, ls=":", alpha=0.7, zorder=5)
axE.plot(360.0 - np.array(locus_f), -np.array(locus_g), color=INK, lw=1.0,
         ls=":", alpha=0.7, zorder=5)
axE.axhline(0.0, color=INK, lw=0.9, alpha=0.5)
axE.set_title(r"$\tan\gamma = \dfrac{e\sin f}{1 + e\cos f}$", fontsize=13, pad=10)
axE.set_xlabel("true anomaly $f$  [deg]")
axE.set_ylabel(r"flight path angle $\gamma$  [deg]")
axE.set_ylim(-70, 70)
axE.set_yticks(range(-60, 61, 30))
axE.legend(frameon=False, fontsize=9.5, ncol=2, loc="lower left",
           handlelength=1.6, labelspacing=0.3, columnspacing=1.0)
axE.text(268, 64, r"markers: $\gamma_{\max} = \arcsin e$ at $\cos f = -e$"
         "\n" r"(i.e. where $r = a$)", ha="center", va="top", fontsize=9.5,
         alpha=0.85)
axE.text(268, 30, "descending:\n" r"$180^\circ < f < 360^\circ$", ha="center",
         va="top", fontsize=10, alpha=0.85)
axE.text(172, -20, "climbing:\n" r"$0 < f < 180^\circ$", ha="right", va="top",
         fontsize=10, alpha=0.85)

# ── panel 2: a and mu do not enter ──────────────────────────────────────────
MU_EARTH = 398600.4418                                # km^3/s^2
e_fix = 0.6
cases = [(7000.0, r"$a = 7\,000$ km", BLUE, "-", 3.4),
         (10000.0, r"$a = 10\,000$ km", GREEN, (0, (6, 4)), 2.2),
         (26600.0, r"$a = 26\,600$ km", MAROON, (0, (1.5, 3)), 2.0)]

for a_km, lab, col, ls, lw in cases:
    p_km = a_km*(1.0 - e_fix**2)
    h = np.sqrt(MU_EARTH*p_km)                        # km^2/s
    vr = MU_EARTH/h*e_fix*np.sin(fr)                  # km/s
    vt = MU_EARTH/h*(1.0 + e_fix*np.cos(fr))          # km/s
    axA.plot(f, np.rad2deg(np.arctan2(vr, vt)), color=col, ls=ls, lw=lw,
             label=lab)

axA.axhline(0.0, color=INK, lw=0.9, alpha=0.5)
axA.set_title(rf"same $e = {e_fix}$, three sizes of orbit", fontsize=13, pad=10)
axA.set_xlabel("true anomaly $f$  [deg]")
axA.set_ylabel(r"$\gamma$  [deg]")
axA.set_ylim(-70, 70)
axA.set_yticks(range(-60, 61, 30))
axA.legend(frameon=False, fontsize=10, loc="lower left")
axA.text(180, 62, "the three curves coincide exactly:\n"
         r"$a$ (and $\mu$) cancel in $v_r/v_\perp$", ha="center", va="top",
         fontsize=10, alpha=0.9)

# ── panel 3: the envelope ───────────────────────────────────────────────────
es = np.linspace(0.0, 1.0, 400)
axM.plot(es, np.rad2deg(np.arcsin(es)), color=MAROON, lw=2.4)
for e, col in zip(ecc, colors):
    if e > 0.0:
        axM.plot([e], [np.rad2deg(np.arcsin(e))], "o", color=col, ms=7, zorder=6)
axM.set_title(r"$\gamma_{\max} = \arcsin e$", fontsize=13, pad=10)
axM.set_xlabel("eccentricity $e$")
axM.set_ylabel(r"$\gamma_{\max}$  [deg]")
axM.set_xlim(0.0, 1.0)
axM.set_ylim(0.0, 95.0)
axM.set_yticks(range(0, 91, 30))
axM.text(0.06, 84, "circle: $\\gamma \\equiv 0$\n"
         "parabola ($e = 1$): $\\gamma \\to 90^\\circ$,\n"
         "reached only at infinity", fontsize=9.5, va="top", alpha=0.85)
axM.plot([0.0167], [np.rad2deg(np.arcsin(0.0167))], "s", color=INK, ms=5)
axM.annotate("Earth's orbit, $e = 0.0167$:\n$\\gamma_{\\max} \\approx 0.96^\\circ$",
             xy=(0.0167, 0.96), xytext=(0.42, 6.0), fontsize=9.5, color=INK,
             arrowprops=dict(arrowstyle="->", color=INK, lw=1.0, alpha=0.8))

for a_ in (axE, axA, axM):
    for sp in ("top", "right"):
        a_.spines[sp].set_visible(False)
for a_ in (axE, axA):
    a_.set_xlim(0, 360)
    a_.set_xticks(range(0, 361, 90))

plt.tight_layout()
fig.savefig("img/fig-flight-path-angle.png", dpi=160, bbox_inches="tight")
plt.close(fig)