"""Q7: closing a 20 deg phase gap at ISS altitude: cost, perigee, and error against k.

Run from the sfm-slides directory:   python3 codes/fig-mq-iss-phasing.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, period

r = Re + 400.0
v = np.sqrt(mu/r); T = period(r)
dth = np.radians(20.0)

def phasing(k):
    Tph = T*(1 - dth/(2*np.pi*k)); aph = (mu*(Tph/(2*np.pi))**2)**(1/3)
    vph = np.sqrt(mu*(2/r - 1/aph))
    return 2*abs(vph - v)*1000, 2*aph - r - Re          # dv_total [m/s], perigee altitude [km]

ks = np.arange(1, 41)
dvs, alts = zip(*[phasing(k) for k in ks])
dvs, alts = np.array(dvs), np.array(alts)

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

axL.semilogy(ks, dvs, color=MAROON, lw=2.6, label=r"total $\Delta v$  [m/s]  (left axis)")
axL.set_xlabel("number of phasing revolutions $k$  (the mission takes about $1.5\\,k$ h)")
axL.set_ylabel(r"total $\Delta v$  [m/s]", color=MAROON)
axL.set_xlim(0, 40); axL.set_ylim(5, 400)
ax2 = axL.twinx()
ax2.plot(ks, alts, color=BLUE, lw=2.6, label="phasing-orbit perigee altitude (right axis)")
ax2.axhline(200, color=INK, lw=1.0, ls=(0, (5, 4)))
ax2.text(20, 168, "perigee floor 200 km (drag)", fontsize=10, color=INK)
ax2.set_ylabel("perigee altitude of the phasing orbit  [km]", color=BLUE)
ax2.set_ylim(-150, 420)
for k, ha in ((2, "left"), (3, "left"), (34, "right")):
    axL.plot(k, phasing(k)[0], "o", color=MAROON, ms=9, zorder=6)
    axL.annotate(f"$k={k}$:  {phasing(k)[0]:.1f} m/s,  perigee {phasing(k)[1]:.0f} km",
                 (k, phasing(k)[0]), (k + (2 if ha == 'left' else -2), phasing(k)[0]*({2:1.9, 3:0.62, 34:2.4}[k])),
                 fontsize=10, color=MAROON, ha=ha, arrowprops=dict(arrowstyle="-", color=MAROON, lw=0.8))
axL.set_title("Cost and safety: more revolutions, less $\\Delta v$ and a higher perigee", fontsize=11.5)
for ax in (axL, ax2):
    ax.spines["top"].set_visible(False)

ks2 = np.linspace(1, 40, 200)
rel = np.degrees(0.005*dth)*np.ones_like(ks2)
absn = np.degrees(6*np.pi*ks2*0.1e-3/v)
rss = np.hypot(rel, absn)
axR.plot(ks2, r*np.radians(rel), color=BLUE, lw=2.4, label=r"relative error 0.5 %  ($=\varepsilon\Delta\theta$, independent of $k$)")
axR.plot(ks2, r*np.radians(absn), color=AMBER, lw=2.4, label=r"absolute error $\sigma=0.1$ m/s  ($=6\pi k\sigma/v$)")
axR.plot(ks2, r*np.radians(rss), color=MAROON, lw=2.8, label="root-sum-square")
for k in (3, 34):
    y = r*np.radians(np.hypot(0.005*np.degrees(dth), np.degrees(6*np.pi*k*0.1e-3/v)))
    axR.plot(k, y, "o", color=MAROON, ms=9, zorder=6)
    axR.annotate(f"$k={k}$: {y:.0f} km", (k, y), (k + 1.5, y + 7), fontsize=10, color=MAROON,
                 arrowprops=dict(arrowstyle="-", color=MAROON, lw=0.8))
axR.set_xlabel("number of phasing revolutions $k$")
axR.set_ylabel("along-track error at the end of the phasing  [km]")
axR.set_xlim(0, 40); axR.set_ylim(0, 75)
axR.legend(frameon=False, fontsize=10, loc="upper left")
axR.set_title(r"Errors for $\Delta\theta=20^\circ$ at 400 km altitude", fontsize=11.5)
despine(axR)

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