"""Ground-track note: repeat ground tracks.

    img/fig-gt-repeat.png   Landsat-like SSO (233 orbits in 16 days) and GPS (2 orbits per sidereal day)
    img/fig-gt-repeat-table.png is not made; the numbers of the validation table are printed.

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

# ---------------------------------------------------------------- validation table
print("REPEAT-ORBIT SOLVER vs published orbits")
rows = [("Landsat 8 / 9", 233, 16, "sso", None, "705 km (osculating, at the equator)"),
        ("Resourcesat / IRS-1C", 341, 24, "sso", None, "817 km"),
        ("Sentinel-2", 143, 10, "sso", None, "786 km"),
        ("Sentinel-1", 175, 12, "sso", None, "693 km"),
        ("Jason-3 (TOPEX class)", 127, 10, "nodal", 66.04, "1336 km, cycle 9.9156 d"),
        ("GPS", 2, 1, "nodal", 55.0, "26 560 km"),
        ("Galileo", 17, 10, "nodal", 56.0, "29 600 km"),
        ("GLONASS", 17, 8, "nodal", 64.8, "25 508 km")]
for name, k, D, kind, inc, pub in rows:
    if kind == "sso":
        a, i = repeat_a_sso(k, D)
    else:
        a, i = repeat_a(k, D, inc), inc
    Od, ud = nodal_rates(a, 0.0, i)
    cyc = D*2*np.pi/(OMEGA_E - Od)/DAY
    print("  %-22s k=%3d D=%2d  a=%9.2f km  h=%8.2f km  i=%7.3f  T=%8.3f min  cycle=%8.4f solar d  spacing=%6.1f km  [published: %s]" %
          (name, k, D, a, a - RE, i, period(a)/60, cyc, 2*np.pi*RE/k, pub))

# ---------------------------------------------------------------- the figure
kL, DL = 233, 16
aL, iL = repeat_a_sso(kL, DL)
dt = 20.0
tL = np.arange(0, DL*DAY + dt, dt)
trL = track(aL, 0.0, iL, 0.0, 0.0, 0.0, tL, j2=True)
print("\nLandsat-like: a = %.2f km, alt %.1f km, i = %.3f deg, T = %.3f min, nodal period T_Omega = %.3f min" %
      (aL, aL - RE, iL, period(aL)/60, 2*np.pi/nodal_rates(aL, 0, iL)[1]/60))
# ascending-node crossings: sign change of latitude, going north
lat, lon, tt = trL["lat_gc"], trL["lon"], trL["t"]
up = np.where((lat[:-1] < 0) & (lat[1:] >= 0))[0]
f = -lat[up]/(lat[up + 1] - lat[up])
lon_asc = (lon[up] + f*((lon[up + 1] - lon[up] + 180) % 360 - 180) + 180) % 360 - 180
t_asc = tt[up] + f*dt
day_of = (t_asc//DAY).astype(int)
order = np.argsort(lon_asc)
print("  %d ascending-node crossings in %.1f days; mean spacing of the sorted longitudes %.4f deg (360/233 = %.4f deg)" %
      (len(up), tt[-1]/DAY, np.mean(np.diff(np.sort(lon_asc))), 360/kL))
gaps = np.diff(np.sort(lon_asc)); print("  spacing min/max: %.4f / %.4f deg -> %.1f / %.1f km at the equator" % (gaps.min(), gaps.max(), gaps.min()*RE*np.pi/180, gaps.max()*RE*np.pi/180))
nb = np.diff(day_of[order])
print("  neighbouring tracks (in longitude) were flown %s days apart (mod 16)" % sorted(set(np.mod(nb, 16))))
print("  along one day: node shift per orbit %.3f deg, %.3f orbits per day" % (np.degrees((OMEGA_E - nodal_rates(aL, 0, iL)[0])*2*np.pi/nodal_rates(aL, 0, iL)[1]), DAY/(2*np.pi/nodal_rates(aL, 0, iL)[1])))
print("  sidelap for a 185 km swath: (185 - %.1f)/185 = %.1f %%" % (2*np.pi*RE/kL, 100*(185 - 2*np.pi*RE/kL)/185))

fig = plt.figure(figsize=(16.5, 9.4))
gs = fig.add_gridspec(2, 2, wspace=0.10, hspace=0.20)
a1 = fig.add_subplot(gs[0, 0]); draw_map(a1)
m1 = tL <= DAY
norb = tt[m1]/(2*np.pi/nodal_rates(aL, 0, iL)[1])
lc = plot_track(a1, lon[m1], trL["lat"][m1], c=norb, cmap="plasma", lw=1.2, vmin=0, vmax=norb.max())
a1.set_title("(a) day 1: 14.6 orbits, node shifts $24.7^\\circ$ west per orbit", fontsize=11)
a2 = fig.add_subplot(gs[0, 1]); draw_map(a2, labels=True)
cmap = plt.get_cmap("turbo")
for d in range(DL):
    m = (tt >= d*DAY) & (tt <= (d + 1)*DAY + dt)
    plot_track(a2, lon[m], trL["lat"][m], color=cmap(d/(DL - 1)), lw=0.7, alpha=0.9)
a2.set_title("(b) all 16 days: 233 orbits fill the globe with a spacing of 172 km", fontsize=11)
sm = cm.ScalarMappable(cmap=cmap, norm=plt.Normalize(1, 16)); cb = fig.colorbar(sm, ax=a2, fraction=0.022, pad=0.01); cb.set_label("day of the cycle", fontsize=10)
a3 = fig.add_subplot(gs[1, 0]); draw_map(a3, lon0=62, lon1=98, lat0=0, lat1=36, grid=5)
for d in range(DL):
    m = (tt >= d*DAY) & (tt <= (d + 1)*DAY + dt)
    plot_track(a3, lon[m], trL["lat"][m], color=cmap(d/(DL - 1)), lw=1.4, alpha=0.95)
a3.set_xticks(range(65, 100, 5)); a3.set_yticks(range(0, 36, 5))
a3.set_title("(c) zoom on India: the track next to yours is flown 7 days later", fontsize=11)
# GPS
aGP = repeat_a(2, 1, 55.0)
tG = np.arange(0, 3*DAY + 1, 30.0)
trG = track(aGP, 0.0, 55.0, 0.0, 0.0, 0.0, tG, j2=True)
a4 = fig.add_subplot(gs[1, 1]); draw_map(a4)
for d, (col, ls, lw) in enumerate(((BLUE, "-", 4.5), (AMBER, "--", 2.6), (MAROON, ":", 1.8))):
    m = (tG >= d*DAY) & (tG <= (d + 1)*DAY + 30)
    for seg in wrap_segments(trG["lon"][m], trG["lat"][m]):
        a4.plot(seg[:, 0], seg[:, 1], color=col, ls=ls, lw=lw, zorder=5)
    a4.plot([], [], color=col, ls=ls, lw=lw, label="day %d" % (d + 1))
a4.legend(loc="lower left", fontsize=9.5, frameon=True, facecolor="white", framealpha=0.9)
a4.set_title("(d) GPS, 26 560 km: 2 orbits per sidereal day, so every day retraces the same track", fontsize=11)
print("GPS: a = %.2f km, T = %.3f min; longitude of the first ascending node at t=0 vs after 1 day (nodal shift %.4f deg)" %
      (aGP, period(aGP)/60, np.degrees((OMEGA_E - nodal_rates(aGP, 0, 55.0)[0])*DAY) - 360))
fig.savefig("img/fig-gt-repeat.png", dpi=150, bbox_inches="tight"); plt.close(fig)