"""Q8: GTO -> GEO.  Left: the velocity triangle at apogee for i = 28.5 deg.  Right: the cost against inclination.

Run from the sfm-slides directory:   python3 codes/fig-mq-gto-geo.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, Re, vis

rpg, rag = Re + 185.0, 42164.0
ag = (rpg + rag)/2
va, vg = vis(rag, ag), np.sqrt(mu/rag)

def dv_comb(inc): return np.sqrt(va**2 + vg**2 - 2*va*vg*np.cos(np.radians(inc)))

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

inc = 28.5
i = np.radians(inc)
vgv = vg*np.array([np.cos(i), np.sin(i)])
O = np.array([0.0, 0.0])
def arrow(ax, p, q, col, lw=2.8, ls="-"):
    ax.annotate("", xy=q, xytext=p, arrowprops=dict(arrowstyle="-|>", color=col, lw=lw, ls=ls, shrinkA=0, shrinkB=0))
arrow(axL, O, np.array([va, 0]), BLUE)
axL.text(va/2, -0.22, r"$v_{\rm apo}=1.596$ km/s (GTO)", color=BLUE, ha="center", va="top", fontsize=11)
arrow(axL, O, vgv, GREEN)
axL.text(vgv[0]*0.55 - 0.1, vgv[1]*0.55 + 0.1, r"$v_{\rm GEO}=3.075$ km/s," "\n" r"in the equatorial plane", color=GREEN, ha="right", fontsize=11)
arrow(axL, np.array([va, 0]), vgv, MAROON, lw=3.2)
mid = (np.array([va, 0]) + vgv)/2
axL.text(mid[0] + 0.45, mid[1] - 0.30, f"one burn:\n$\\Delta v$ = {dv_comb(inc):.3f} km/s", color=MAROON, fontsize=11.5, va="top")
beta = np.degrees(np.arctan2(vgv[1], vgv[0] - va))
th_ = np.linspace(0, np.radians(beta), 40)
axL.plot(va + 0.55*np.cos(th_), 0.55*np.sin(th_), color=MAROON, lw=1.4)
axL.text(va - 0.05, 0.42, f"$\\beta={beta:.0f}^\\circ$ from the velocity", color=MAROON, fontsize=10, ha="right")
ar = np.linspace(0, i, 40)
axL.plot(1.0*np.cos(ar), 1.0*np.sin(ar), color=INK, lw=1.2)
axL.text(1.12, 0.18, f"$i={inc}^\\circ$", color=INK, fontsize=11)
axL.set_xlim(-0.5, 5.2); axL.set_ylim(-0.9, 2.1)
axL.set_aspect("equal"); axL.axis("off")
axL.set_title("At apogee: circularise and change the plane in one burn", fontsize=11.5)

incs = np.linspace(0, 35, 200)
dvc = np.array([dv_comb(x) for x in incs])*1000
dsep = ((vg - va) + 2*vg*np.sin(np.radians(incs)/2))*1000
axR.plot(incs, dvc, color=MAROON, lw=2.8, label="plane change folded into the circularisation burn")
axR.plot(incs, dsep, color=INK, lw=2.0, ls=(0, (5, 4)), label="circularise, then change the plane (two burns)")
for lab, x, col in (("Kourou\n7$^\\circ$", 7.0, GREEN), ("Sriharikota\n19.2$^\\circ$", 19.2, AMBER), ("Cape Canaveral\n28.5$^\\circ$", 28.5, PURPLE)):
    axR.axvline(x, color=col, lw=1.0, ls=":")
    axR.plot(x, dv_comb(x)*1000, "o", color=col, ms=9, zorder=6)
    axR.text(x - 0.4, 780, lab, color=col, ha="right", fontsize=10)
    axR.text(x + 0.5, dv_comb(x)*1000 - 120, f"{dv_comb(x)*1000:.0f} m/s", color=col, fontsize=10, va="top")
axR.set_xlabel("inclination of the GTO  [deg]")
axR.set_ylabel(r"$\Delta v$ at apogee  [m/s]")
axR.set_xlim(0, 35); axR.set_ylim(700, 3400)
axR.legend(frameon=False, fontsize=10.5, loc="upper left")
axR.set_title("GTO 185 $\\times$ 35 786 km: a launch site fixes the plane-change bill", fontsize=11.5)
despine(axR)

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