"""Q10: the Delta v budget of a GEO satellite, and what happens when the numbers are wrong.

Left: propellant used against year (chemical, Isp = 320 s; electric, Isp = 1800 s) with the
+/-3 % gauging band on the loaded mass.  Right: lifetime against the N-S station-keeping rate.

Run from the sfm-slides directory:   python3 codes/fig-mq-geo-budget.py
"""
import numpy as np
import matplotlib.pyplot as plt
from _manstyle import CREAM, INK, MAROON, BLUE, GREEN, AMBER, PURPLE, despine
from _mqcommon import g0

NS, EW, DISP, LIFE = 50.0, 3.0, 11.0, 15
m0 = 3000.0
Dv = LIFE*(NS + EW) + DISP                       # 806 m/s
ve_c, ve_e = 320*g0*1000, 1800*g0*1000
mp_c = m0*(1 - np.exp(-Dv/ve_c))

fig, (axL, axR) = plt.subplots(1, 2, figsize=(14.4, 5.6))

yrs = np.linspace(0, LIFE, 200)
dv_y = yrs*(NS + EW)
axL.plot(yrs, m0*(1 - np.exp(-dv_y/ve_c)), color=MAROON, lw=2.6, label=r"chemical, $I_{sp}=320$ s")
axL.plot(yrs, m0*(1 - np.exp(-dv_y/ve_e)), color=GREEN, lw=2.6, label=r"electric, $I_{sp}=1800$ s")
axL.fill_between([0, LIFE + 1.2], 0.97*mp_c, 1.03*mp_c, color=AMBER, alpha=0.25)
axL.text(0.4, 1.035*mp_c + 8, f"loaded propellant {mp_c:.0f} kg $\\pm$ 3 % gauging ($\\pm$20 kg)", color=AMBER, fontsize=10.5)
axL.plot(LIFE, mp_c, "o", color=MAROON, ms=8)
axL.annotate(f"806 m/s\n{mp_c:.0f} kg = 22.7 %", (LIFE, mp_c), (LIFE - 7.5, mp_c - 80), fontsize=10.5, color=MAROON,
             arrowprops=dict(arrowstyle="-", color=MAROON, lw=0.8))
axL.plot(LIFE, m0*(1 - np.exp(-Dv/ve_e)), "o", color=GREEN, ms=8)
axL.annotate("134 kg", (LIFE, m0*(1 - np.exp(-Dv/ve_e))), (LIFE - 3.2, 215), fontsize=10.5, color=GREEN,
             arrowprops=dict(arrowstyle="-", color=GREEN, lw=0.8))
axL.set_xlabel("years in orbit  (N-S 50 m/s/yr + E-W 3 m/s/yr; the 11 m/s disposal burn at the end)")
axL.set_ylabel(r"propellant used  [kg], 3000 kg satellite")
axL.set_xlim(0, LIFE + 1.2); axL.set_ylim(0, 800)
axL.legend(frameon=False, fontsize=10.5, loc="center left")
axL.set_title("A 15-year budget of 806 m/s", fontsize=11.5)
despine(axL)

rate = np.linspace(40, 62, 200)
life = (Dv - DISP)/(rate + EW)
axR.plot(rate, life, color=BLUE, lw=2.6, label="propellant for 806 m/s, Isp 320 s")
life_lowisp = (0.98*ve_c*np.log(m0/(m0 - mp_c)) - DISP)/(rate + EW)
axR.plot(rate, life_lowisp, color=PURPLE, lw=2.0, ls=(0, (5, 4)), label="same load, Isp 2 % below specification")
axR.axhline(15, color=INK, lw=1.0, ls=":")
for r_, col in ((50, INK), (55, MAROON)):
    l_ = (Dv - DISP)/(r_ + EW)
    axR.plot(r_, l_, "o", color=col, ms=9, zorder=6)
    axR.annotate(f"{r_} m/s/yr: {l_:.1f} yr", (r_, l_), (r_ + 0.8, l_ + 1.0), fontsize=10.5, color=col,
                 arrowprops=dict(arrowstyle="-", color=col, lw=0.8))
axR.set_xlabel("north-south station-keeping rate  [m/s per year]")
axR.set_ylabel("lifetime  [years]")
axR.set_xlim(40, 62); axR.set_ylim(10, 20)
axR.legend(frameon=False, fontsize=10.5, loc="upper right")
axR.set_title("A 10 % higher N-S rate costs about 1.3 years of life", fontsize=11.5)
despine(axR)

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