"""Q2: a finite burn centred on periapsis gives up a little energy per m/s.

Left: the burn arc on an orbit of e = 0.6.  Right: the efficiency eta = <v>/v_p of a
burn that spans +/- theta of true anomaly, for three eccentricities (lines: the
small-arc formula 1 - e theta^2 / [6 (1+e)^2]; dots: exact, by integrating the orbit).
The two real burns of Q5 and Q6 are marked.

Run from the sfm-slides directory:   python3 codes/fig-mq-finite-burn.py
"""
import numpy as np
import matplotlib.pyplot as plt
from _manstyle import CREAM, INK, MAROON, BLUE, GREEN, AMBER, PURPLE, despine
from _mqcommon import mu, conic_xy, eta_exact, eta_formula

fig, (axL, axR) = plt.subplots(1, 2, figsize=(14.4, 5.6), gridspec_kw=dict(width_ratios=[1, 1.15]))

# ---- A : the arc -------------------------------------------------------------
e, rp = 0.6, 1.0
a = rp/(1 - e)
x, y = conic_xy(a, e)
axL.plot(x, y, color=BLUE, lw=2.2)
th = np.radians(40)
x2, y2 = conic_xy(a, e, fmax=th, n=200)
axL.plot(x2, y2, color=AMBER, lw=7, solid_capstyle="round", alpha=0.9, zorder=4)
for f in np.linspace(-th, th, 5):
    r = a*(1 - e*e)/(1 + e*np.cos(f))
    v = np.array([-np.sin(f)*(1 + e*np.cos(f)) + 0*np.cos(f), np.cos(f)*(1 + e*np.cos(f)) + e])  # direction of motion (up to scale)
    v = np.array([-np.sin(f), e + np.cos(f)])
    v = v/np.hypot(*v)
    axL.annotate("", xy=(r*np.cos(f) + 0.42*v[0], r*np.sin(f) + 0.42*v[1]), xytext=(r*np.cos(f), r*np.sin(f)),
                 arrowprops=dict(arrowstyle="-|>", color=MAROON, lw=1.6))
axL.plot(0, 0, "o", color=INK, ms=9)
axL.text(-0.15, -0.38, "Earth", fontsize=10, color=INK)
axL.plot(rp, 0, "o", color=AMBER, ms=8, zorder=6)
axL.text(rp + 0.1, -0.30, "periapsis", fontsize=10, color=AMBER)
axL.text(1.45, 0.9, r"burn arc" "\n" r"$-\theta \ldots +\theta$," "\n" r"thrust along $\mathbf{v}$", color=MAROON, fontsize=11)
axL.set_aspect("equal"); axL.axis("off")
axL.set_xlim(-4.4, 2.6); axL.set_ylim(-2.3, 2.3)
axL.set_title(r"$e=0.6$: the burn is centred in time on periapsis (arc drawn for $\theta=40^\circ$)", fontsize=11.5)

# ---- B : efficiency ----------------------------------------------------------
rpk = 7000.0
for e, col in ((0.3, GREEN), (0.6, BLUE), (0.9, MAROON)):
    thd = np.linspace(0, 75, 200)
    axR.plot(thd, 100*eta_formula(e, np.radians(thd)), color=col, lw=2.3, label=f"$e={e}$  (formula)")
    ra = rpk*(1 + e)/(1 - e)
    pts = []
    for thd_ in (20, 40, 60):
        a = (rpk + ra)/2; p = a*(1 - e*e); h = np.sqrt(mu*p)
        f, t = 0.0, 0.0
        while f < np.radians(thd_):
            f += h/(p/(1 + e*np.cos(f)))**2*0.5; t += 0.5
        pts.append((thd_, 100*eta_exact(rpk, ra, 2*t)))
    axR.plot(*zip(*pts), "o", color=col, ms=7, mfc=CREAM, mew=2)
axR.plot([], [], "o", color=INK, mfc=CREAM, mew=2, label="exact (integrated)")
# the two real burns
for lab, e, thd, dx, dy in (("Mangalyaan TMI  (e = 0.94, $\\pm$53$^\\circ$)", 0.936, 53.3, -43, -2.3),
                            ("Chandrayaan-3 EBN  (e = 0.73, $\\pm$36$^\\circ$)", 0.734, 36.5, 4, 1.6)):
    y = 100*eta_formula(e, np.radians(thd))
    axR.plot(thd, y, "s", color=AMBER, ms=9, zorder=6)
    axR.annotate(lab, (thd, y), (thd + dx, y + dy), fontsize=10, color=AMBER,
                 arrowprops=dict(arrowstyle="-", color=AMBER, lw=0.8))
axR.set_xlabel(r"half-arc $\theta$ of true anomaly spanned by the burn  [deg]")
axR.set_ylabel(r"energy gained per $\Delta v$, as % of the impulsive value ($v_p\Delta v$)")
axR.set_xlim(0, 75); axR.set_ylim(86, 100.5)
axR.legend(frameon=False, fontsize=10, loc="lower left")
axR.set_title(r"$\eta = 1-\dfrac{e\,\theta^{2}}{6\,(1+e)^{2}}$: the loss grows as the square of the burn duration", fontsize=11.5)
despine(axR)

plt.tight_layout()
fig.savefig("img/fig-mq-finite-burn.png", dpi=160, bbox_inches="tight")
plt.close(fig)
print("wrote img/fig-mq-finite-burn.png")