"""Ground-track note: what the ground track is used for.

    img/fig-gt-footprint.png   coverage circle geometry, coverage radius vs altitude, ground-station passes vs latitude
    img/fig-gt-launch.png      launch azimuth vs inclination, and a polar launch from Sriharikota

Run from the deck directory:   python codes/fig-gt-apps.py
"""
import numpy as np
import matplotlib.pyplot as plt
from matplotlib.patches import Circle, Wedge, Arc
from _gt import *
from _timeutil import RE, OMEGA_E, MU_E, D2R, DAY
from _manstyle import CREAM, INK, MAROON, BLUE, GREEN, AMBER, PURPLE, despine

# ================================================================== footprint / visibility
fig = plt.figure(figsize=(16.2, 5.0))
gs = fig.add_gridspec(1, 3, width_ratios=[1, 1.05, 1.15], wspace=0.28)

# (a) geometry, drawn with an exaggerated altitude so that it is legible
ax = fig.add_subplot(gs[0]); ax.set_aspect("equal"); ax.axis("off")
hR, eps = 0.45, np.radians(20)
eta = np.arcsin(np.cos(eps)/(1 + hR)); lam = np.pi/2 - eps - eta
ax.add_patch(Circle((0, 0), 1, color="#cfe3f5", ec=INK, lw=1.4))
G = np.array([0, 1.0]); Sx = (1 + hR)*np.array([np.sin(lam), np.cos(lam)])
ax.plot([0, G[0]], [0, G[1]], color=INK, lw=1.0); ax.plot([0, Sx[0]], [0, Sx[1]], color=INK, lw=1.0)
ax.plot([G[0], Sx[0]], [G[1], Sx[1]], color=MAROON, lw=2.2)
ax.plot([-0.6, 0.9], [1, 1], color=INK, lw=0.8, ls="--")
nadir_pt = Sx/np.linalg.norm(Sx)
ax.plot([Sx[0], nadir_pt[0]], [Sx[1], nadir_pt[1]], color=BLUE, lw=1.6, ls=":")
ax.plot(*G, "o", color=MAROON, ms=8); ax.text(-0.05, 1.05, "ground\nstation", ha="right", fontsize=10, color=MAROON)
ax.plot(*Sx, "s", color=BLUE, ms=10); ax.text(Sx[0] + 0.05, Sx[1] + 0.05, "satellite", color=BLUE, fontsize=10.5)
ax.plot(0, 0, "o", color=INK, ms=4); ax.text(0.03, -0.12, "O", fontsize=10)
ax.add_patch(Arc(G, 0.55, 0.55, theta1=0, theta2=np.degrees(np.arctan2(Sx[1] - G[1], Sx[0] - G[0])), color=MAROON, lw=1.4))
ax.text(0.30, 1.13, r"$\varepsilon$", color=MAROON, fontsize=14)
ax.add_patch(Arc((0, 0), 0.6, 0.6, theta1=90 - np.degrees(lam), theta2=90, color=GREEN, lw=1.6))
ax.text(0.09, 0.36, r"$\lambda$", color=GREEN, fontsize=14)
ax.add_patch(Arc(Sx, 0.5, 0.5, theta1=np.degrees(np.arctan2(-Sx[1], -Sx[0])) - 8 - np.degrees(eta) + 8, theta2=np.degrees(np.arctan2(G[1] - Sx[1], G[0] - Sx[0])) + 0.0, color=BLUE, lw=1.4))
ax.text(Sx[0] - 0.10, Sx[1] - 0.45, r"$\eta$", color=BLUE, fontsize=14)
ax.text(-0.08, 0.5, r"$R$", fontsize=11, ha="right", va="center"); ax.text(0.5*Sx[0] + 0.07, 0.5*Sx[1] - 0.10, r"$R+h$", fontsize=11)
ax.text(0.0, -1.35, r"$\sin\eta=\dfrac{R}{R+h}\cos\varepsilon,\qquad \lambda=90^\circ-\varepsilon-\eta$", ha="center", fontsize=12)
ax.set_xlim(-1.3, 1.9); ax.set_ylim(-1.7, 1.75)
ax.set_title("(a) What a station (or a sensor) can see", fontsize=11.5)

# (b) coverage half-angle vs altitude
b = fig.add_subplot(gs[1])
h = np.logspace(np.log10(300), np.log10(60000), 300)
for e_, col in ((0, INK), (5, BLUE), (10, GREEN), (30, AMBER), (60, MAROON)):
    b.semilogx(h, earth_central_angle(h, e_), color=col, lw=2.2, label="$\\varepsilon_{\\min}=%d^\\circ$" % e_)
for name, hh, dy in (("LEO 800 km", 800, 4), ("GPS 20 200 km", 20200, 4), ("GEO 35 786 km", 35786, -12)):
    b.axvline(hh, color=INK, lw=0.6, ls=":"); b.text(hh*1.06, 62 + dy, name, rotation=90, fontsize=8.8, va="top")
b.set_xlabel("altitude $h$ (km)"); b.set_ylabel("half-angle $\\lambda$ of the coverage circle (deg)")
b.set_ylim(0, 90); b.legend(frameon=False, fontsize=9, loc="upper left"); despine(b)
b.set_title("(b) Higher orbit, bigger footprint", fontsize=11.5)
print("coverage half-angle: 800 km, eps 10 deg = %.2f deg (radius %.0f km); 800 km, eps 0 = %.2f deg; GEO eps 0 = %.2f deg; GEO eps 5 = %.2f deg; 420 km eps 5 = %.2f deg" % (
    earth_central_angle(800, 10), earth_central_angle(800, 10)*np.pi/180*RE, earth_central_angle(800, 0), earth_central_angle(35786, 0), earth_central_angle(35786, 5), earth_central_angle(420, 5)))
print("GEO covers %.1f %% of the sphere's area at eps=0 (1 - cos(lambda))/2" % (100*(1 - np.cos(np.radians(earth_central_angle(35786, 0))))/2))

# (c) passes per day vs station latitude
def passes(a, inc, lat_st, lon_st, emin=5.0, days=10, dt=10.0, raan=0.0):
    t = np.arange(0, days*DAY, dt)
    tr = track(a, 0.0, inc, raan, 0.0, 0.0, t, j2=True)
    el, _ = elevation(tr["ecef"], lat_st, lon_st)
    vis = el > emin
    edges = np.diff(vis.astype(int))
    starts, ends = np.where(edges == 1)[0], np.where(edges == -1)[0]
    if vis[0]: starts = np.r_[0, starts]
    if vis[-1]: ends = np.r_[ends, len(vis) - 1]
    n = min(len(starts), len(ends))
    dur = (ends[:n] - starts[:n])*dt
    return n/days, vis.sum()*dt/days/60, (dur.mean()/60 if n else 0.0)
c = fig.add_subplot(gs[2])
lats = np.arange(0, 86, 5.0)
aS = RE + 800; iS = sso_inclination(aS)
res = {}
for lab, a_, i_, col in (("SSO, 800 km, $i=98.6^\\circ$", aS, iS, MAROON), ("ISS-like, 420 km, $i=51.6^\\circ$", RE + 420, 51.64, BLUE)):
    pp = np.array([passes(a_, i_, la, 20.0)[0] for la in lats])
    res[lab] = pp
    c.plot(lats, pp, "o-", color=col, lw=2.2, ms=5, label=lab)
c.set_xlabel("latitude of the ground station (deg)"); c.set_ylabel("passes per day above $5^\\circ$ elevation"); despine(c)
c.legend(frameon=False, fontsize=9.3, loc="upper left")
c.set_title("(c) How often a station sees the satellite", fontsize=11.5)
fig.savefig("img/fig-gt-footprint.png", dpi=160, bbox_inches="tight"); plt.close(fig)

print("\nSTATION PASSES (10-day average, elevation > 5 deg), SSO 800 km / ISS 420 km")
for name, la, lo in (("Thiruvananthapuram", 8.5, 76.9), ("Bengaluru", 13.0, 77.6), ("Delhi", 28.6, 77.2), ("Svalbard", 78.2, 15.4), ("Equator", 0.0, 20.0)):
    p1 = passes(aS, iS, la, lo); p2 = passes(RE + 420, 51.64, la, lo)
    print("  %-20s SSO: %.2f passes/day, %.1f min/day, %.1f min/pass   |   ISS: %.2f passes/day, %.1f min/day, %.1f min/pass" % ((name,) + p1 + p2))

# ================================================================== launch
fig = plt.figure(figsize=(14.6, 5.4))
gs = fig.add_gridspec(1, 2, width_ratios=[1, 1.15], wspace=0.22)
ax = fig.add_subplot(gs[0])
sites = (("Kourou 5.2°N", 5.24, GREEN), ("Sriharikota 13.7°N", 13.72, MAROON), ("Cape Canaveral 28.5°N", 28.5, BLUE), ("Baikonur 45.9°N", 45.9, AMBER))
inc = np.linspace(0.5, 179.5, 700)
for name, phi, col in sites:
    sinaz = np.cos(np.radians(inc))/np.cos(np.radians(phi))
    ok = np.abs(sinaz) <= 1
    az = np.full_like(inc, np.nan)
    az[ok] = np.degrees(np.arcsin(sinaz[ok]))            # north-going branch; the south-going one is 180 deg minus this
    ax.plot(inc[ok], az[ok], color=col, lw=2.2, label=name)
    ax.plot([phi, 180 - phi], [90, -90], "o", color=col, ms=7)
    print("  %-24s minimum direct inclination = %.2f deg; Earth's rotation contributes %.1f m/s eastward" % (name, phi, 465.1*np.cos(np.radians(phi))))
ax.axhline(0, color=INK, lw=0.5)
ax.set_xlim(0, 180); ax.set_ylim(-100, 100); ax.set_xlabel("target inclination $i$ (deg)"); ax.set_ylabel("northbound launch azimuth $A_z$ (deg east of north)")
ax.set_yticks([-90, -45, 0, 45, 90]); despine(ax)
ax.text(90, -8, "polar: due north (or, going south, due south)", fontsize=9, ha="center", va="top")
ax.text(2, 96, "no direct launch for $i<\\varphi_L$\n(dots: the limiting cases)", fontsize=9, va="top")
ax.legend(frameon=False, fontsize=9, loc="lower left")
ax.set_title("(a) $\\cos i=\\cos\\varphi_L\\,\\sin A_z$: the site latitude fixes the lowest inclination", fontsize=11)
# SHAR to SSO
bm = fig.add_subplot(gs[1]); draw_map(bm, lon0=68, lon1=92, lat0=0, lat1=22, grid=2)
lat_L, lon_L = 13.72, 80.23
i_s = sso_inclination(RE + 820)
u_L = np.pi - np.arcsin(np.sin(np.radians(lat_L))/np.sin(np.radians(i_s)))
Om = lon_L - np.degrees(np.arctan2(np.cos(np.radians(i_s))*np.sin(u_L), np.cos(u_L)))
tt = np.arange(-900, 900, 5.0)
aL = RE + 820
u0 = np.degrees(u_L)
tr = track(aL, 0.0, i_s, Om, 0.0, u0, tt)
plot_track(bm, tr["lon"], tr["lat_gc"], color=MAROON, lw=2.6)
bm.plot(lon_L, lat_L, "*", color=AMBER, mec=INK, ms=17, zorder=9); bm.text(lon_L + 0.5, lat_L + 0.6, "Sriharikota", fontsize=10, weight="bold")
k0 = np.argmin(abs(tt))
dl, dlon = tr["lat_gc"][k0 + 1] - tr["lat_gc"][k0], ((tr["lon"][k0 + 1] - tr["lon"][k0] + 180) % 360 - 180)*np.cos(np.radians(lat_L))
head = (np.degrees(np.arctan2(dlon, dl))) % 360
saz = np.degrees(np.arcsin(np.cos(np.radians(i_s))/np.cos(np.radians(lat_L)))) % 360
azin = 180 - np.degrees(np.arcsin(np.cos(np.radians(i_s))/np.cos(np.radians(lat_L))))*-1
print("\nSHAR polar/SSO launch: i = %.2f deg (820 km). sin Az = cos i / cos phi = %.4f -> Az = %.2f deg or %.2f deg (inertial)" % (
    i_s, np.cos(np.radians(i_s))/np.cos(np.radians(lat_L)), np.degrees(np.arcsin(np.cos(np.radians(i_s))/np.cos(np.radians(lat_L)))), 180 - np.degrees(np.arcsin(np.cos(np.radians(i_s))/np.cos(np.radians(lat_L))))))
print("  heading of the GROUND track over the launch site = %.2f deg (Earth's rotation included)" % head)
bm.text(80.6, 7.6, "Sri\nLanka", fontsize=9, style="italic", ha="center")
bm.set_title("(b) Ground track through the pad for a sun-synchronous orbit: heading %.0f° (south-south-west)" % head, fontsize=10.5)
bm.set_xticks(range(68, 93, 4)); bm.set_yticks(range(0, 23, 4))
fig.savefig("img/fig-gt-launch.png", dpi=160, bbox_inches="tight"); plt.close(fig)