"""Figure for topocentric.qmd (img/fig-topo-pass.png): an ISS-like pass over Bengaluru.

Left: the sky plot (azimuth round the edge, zenith distance outwards) with the 5 degree mask.
Right: elevation, range and range rate (hence the Doppler shift) against time since acquisition of signal.
The orbit is circular at 420 km, i = 51.64 deg, Omega = 166 deg; theta_G = 100 deg and argument of latitude 0 at t = 0.

Run from the sfm-slides directory:   python3 codes/fig-topo-pass.py
"""
import numpy as np
import matplotlib.pyplot as plt
from scipy.optimize import brentq, minimize_scalar
from _manstyle import INK, MAROON, BLUE, GREEN, AMBER, PURPLE, despine
from _topo import circular_state, observe, WE, C

phi, lam, h = np.radians(12.97), np.radians(77.59), 0.92
a, inc, raan, tg0 = 6378.137 + 420, np.radians(51.64), np.radians(166.0), np.radians(100.0)


def obs(t):
    r, v = circular_state(a, inc, raan, 0.0, t)
    return observe(r, v, tg0 + WE * t, phi, lam, h)


el = lambda t: obs(t)[1]
tm = minimize_scalar(lambda t: -el(t), bounds=(-300, 600), method="bounded", options={"xatol": 1e-6}).x
mask = np.radians(5)
t1 = brentq(lambda t: el(t) - mask, tm - 600, tm)
t2 = brentq(lambda t: el(t) - mask, tm, tm + 600)
ts = np.linspace(t1, t2, 800)
O = np.array([obs(t) for t in ts])
az, ele, rng, daz, dele, drng = O.T
fig = plt.figure(figsize=(13.6, 5.6))
ax = fig.add_subplot(1, 3, 1, projection="polar")
ax.set_theta_zero_location("N"); ax.set_theta_direction(-1)
ax.plot(az, 90 - np.degrees(ele), color=MAROON, lw=2.6)
ax.scatter([az[0], az[-1]], [90 - np.degrees(ele[0]), 90 - np.degrees(ele[-1])], color=[GREEN, AMBER], s=50, zorder=5)
im = np.argmax(ele)
ax.scatter([az[im]], [90 - np.degrees(ele[im])], color=BLUE, s=50, zorder=5)
ax.text(az[0] + 0.05, 90 - np.degrees(ele[0]) - 12, "AOS", color=GREEN, fontsize=10, ha="left")
ax.text(az[-1] - 0.1, 90 - np.degrees(ele[-1]) + 2, "LOS", color=AMBER, fontsize=10, ha="right")
ax.text(az[im] + 0.25, 90 - np.degrees(ele[im]) - 2, f"max {np.degrees(ele[im]):.1f}$^\\circ$", color=BLUE, fontsize=10)
ax.set_ylim(0, 90); ax.set_yticks([0, 30, 60, 90]); ax.set_yticklabels(["90°", "60°", "30°", "horizon"], fontsize=8)
ax.set_title("Sky plot (elevation labelled)", fontsize=11.5, pad=14)
tt = ts - t1
b = fig.add_subplot(1, 3, 2)
b.plot(tt, np.degrees(ele), color=BLUE, lw=2.4); b.set_ylabel("elevation [deg]", color=BLUE); b.set_xlabel("time since AOS [s]")
b2 = b.twinx(); b2.plot(tt, rng, color=MAROON, lw=2.2); b2.set_ylabel("range [km]", color=MAROON)
b.set_title("Elevation and range", fontsize=11.5); despine(b)
c = fig.add_subplot(1, 3, 3)
c.plot(tt, drng, color=PURPLE, lw=2.4); c.axhline(0, color=INK, lw=0.8)
c.set_xlabel("time since AOS [s]"); c.set_ylabel(r"range rate $\dot\rho$ [km/s]", color=PURPLE)
c2 = c.twinx(); c2.plot(tt, -2.2e9 * drng / C / 1e3, color=AMBER, lw=0, alpha=0)
c2.set_ylabel("Doppler shift at 2.2 GHz [kHz]", color=AMBER)
lo, hi = c.get_ylim(); c2.set_ylim(-2.2e9 * lo / C / 1e3, -2.2e9 * hi / C / 1e3)
c.set_title("Range rate and Doppler (one-way)", fontsize=11.5); despine(c)
fig.tight_layout()
fig.savefig("img/fig-topo-pass.png", dpi=150, bbox_inches="tight")
print(f"AOS {t1:.2f} s, LOS {t2:.2f} s, duration {t2 - t1:.1f} s; max elevation {np.degrees(ele[im]):.2f} deg at {tm - t1:.1f} s after AOS")
print(f"AOS az {np.degrees(az[0]):.1f}, range {rng[0]:.1f}; LOS az {np.degrees(az[-1]):.1f}, range {rng[-1]:.1f}; min range {rng.min():.1f}; max |rdot| {abs(drng).max():.3f} km/s")