"""Ground-track note: the orbit stays put, the Earth turns, and the ground track appears.

Left panel   an inertial view.  The orbit plane (blue ring) is fixed in space; the globe turns beneath it
             (about 15 deg/h).  The red trace is painted on the rotating surface.
Right panel  the same trace on an Earth-fixed map: latitude and longitude of the sub-satellite point.

Run from the deck directory:   python codes/anim-gt-build.py [i_deg altitude_km n_orbits]
Writes img/anim-gt-build.gif
"""
import sys
import numpy as np
import matplotlib.pyplot as plt
from matplotlib.patches import Circle
from _gt import *
from _gif import save_gif
from _timeutil import RE, OMEGA_E, D2R
from _manstyle import CREAM, INK, MAROON, BLUE, GREEN, AMBER, PURPLE

inc = float(sys.argv[1]) if len(sys.argv) > 1 else 60.0
alt = float(sys.argv[2]) if len(sys.argv) > 2 else 1200.0
norb = float(sys.argv[3]) if len(sys.argv) > 3 else 4.0
a = RE + alt
T = period(a)
fps, seconds, hold = 15, 15, 20
n = fps*seconds
times = np.concatenate([np.linspace(0, norb*T, n), np.full(hold, norb*T)])
tdense = np.linspace(0, norb*T, 2500)
TR = track(a, 0.0, inc, 20.0, 0.0, 0.0, tdense)
tr_all = track(a, 0.0, inc, 20.0, 0.0, 0.0, times)
Om = np.radians(20.0)

# camera
el, az = np.radians(35), np.radians(-20)
cam = np.array([np.cos(el)*np.cos(az), np.cos(el)*np.sin(az), np.sin(el)])
right = np.cross([0, 0, 1], cam); right /= np.linalg.norm(right)
upv = np.cross(cam, right)

def R3(th):
    c, s = np.cos(th), np.sin(th)
    return np.array([[c, -s, 0], [s, c, 0], [0, 0, 1]])

def project(P, rt, up, cm):
    return P@rt, P@up, P@cm

land_ecef = [np.column_stack([np.cos(np.radians(r[:, 1]))*np.cos(np.radians(r[:, 0])),
                              np.cos(np.radians(r[:, 1]))*np.sin(np.radians(r[:, 0])),
                              np.sin(np.radians(r[:, 1]))]) for r in land()]

fig = plt.figure(figsize=(13.4, 5.5))
gs = fig.add_gridspec(1, 2, width_ratios=[1, 1.75], wspace=0.06)
g = fig.add_subplot(gs[0]); g.set_aspect("equal"); g.axis("off"); g.set_xlim(-1.75, 1.75); g.set_ylim(-1.75, 1.75)
m = fig.add_subplot(gs[1]); draw_map(m)
g.add_patch(Circle((0, 0), 1, color="#cfe3f5", ec=INK, lw=1.3, zorder=1))
land_lines = [g.plot([], [], color="#7d7660", lw=0.9, zorder=2)[0] for _ in land_ecef]
trace_g, = g.plot([], [], color=MAROON, lw=2.4, zorder=6)
# orbit ring (inertial)
ang = np.linspace(0, 2*np.pi, 400)
Ro = a/RE
ring = Ro*np.array([np.cos(ang)*np.cos(Om) - np.sin(ang)*np.cos(np.radians(inc))*np.sin(Om),
                    np.cos(ang)*np.sin(Om) + np.sin(ang)*np.cos(np.radians(inc))*np.cos(Om),
                    np.sin(ang)*np.sin(np.radians(inc))]).T
rx, ry, rz = project(ring, right, upv, cam)
hidden = (rz < 0) & (np.hypot(rx, ry) < 1)
rxf, ryf = rx.copy(), ry.copy(); rxf[hidden] = np.nan; ryf[hidden] = np.nan
g.plot(rxf, ryf, color=BLUE, lw=1.4, zorder=5)
sat_g, = g.plot([], [], "o", color=BLUE, ms=9, mec=INK, zorder=9)
nadir_g, = g.plot([], [], color=BLUE, lw=1.0, ls=":", zorder=8)
sub_g, = g.plot([], [], "o", color=MAROON, ms=6, mec="white", zorder=9)
g.text(-1.7, 1.7, "Inertial view: the orbit plane (blue) is fixed,\nthe Earth turns beneath it", fontsize=10, va="top")
trace_m, = m.plot([], [], color=MAROON, lw=2.0, zorder=6)
sub_m, = m.plot([], [], "o", color=MAROON, ms=8, mec="white", zorder=9)
info = m.text(0.01, 0.02, "", transform=m.transAxes, fontsize=10.5, family="monospace", va="bottom", zorder=12,
              bbox=dict(fc="white", ec="none", alpha=0.85, pad=3))
m.set_title("Earth-fixed map: latitude and longitude of the sub-satellite point", fontsize=11)

def frame(k):
    t = times[k]
    th = OMEGA_E*t
    rt, up, cm = R3(th).T@right, R3(th).T@upv, R3(th).T@cam        # camera axes expressed in the Earth-fixed frame
    for line, P in zip(land_lines, land_ecef):
        x, y, z = P@rt, P@up, P@cm
        x = np.where(z > 0, x, np.nan); y = np.where(z > 0, y, np.nan)
        line.set_data(x, y)
    # trace so far (ECEF unit vectors)
    mk = tdense <= t
    lat, lon = np.radians(TR["lat_gc"][mk]), np.radians(TR["lon"][mk])
    P = np.column_stack([np.cos(lat)*np.cos(lon), np.cos(lat)*np.sin(lon), np.sin(lat)])
    x, y, z = P@rt, P@up, P@cm
    x = np.where(z > 0.02, x, np.nan); y = np.where(z > 0.02, y, np.nan)
    trace_g.set_data(x, y)
    # satellite (inertial)
    sat = tr_all["eci"][:, k]/RE
    sx, sy, sz = sat@right, sat@upv, sat@cam
    hid = (sz < 0) and (np.hypot(sx, sy) < 1)
    sat_g.set_data([] if hid else [sx], [] if hid else [sy])
    e_ = np.array([np.cos(np.radians(tr_all["lat_gc"][k]))*np.cos(np.radians(tr_all["lon"][k])),
                   np.cos(np.radians(tr_all["lat_gc"][k]))*np.sin(np.radians(tr_all["lon"][k])),
                   np.sin(np.radians(tr_all["lat_gc"][k]))])
    ex, ey, ez = e_@rt, e_@up, e_@cm
    sub_g.set_data([ex] if ez > 0 else [], [ey] if ez > 0 else [])
    nadir_g.set_data([sx, ex] if (not hid and ez > 0) else [], [sy, ey] if (not hid and ez > 0) else [])
    # map
    lonm = (TR["lon"][mk] + 180) % 360 - 180
    latm = TR["lat"][mk]
    cut = np.where(np.abs(np.diff(lonm)) > 180)[0] + 1
    lonp, latp = np.insert(lonm.astype(float), cut, np.nan), np.insert(latm.astype(float), cut, np.nan)
    trace_m.set_data(lonp, latp)
    sub_m.set_data([(tr_all["lon"][k] + 180) % 360 - 180], [tr_all["lat"][k]])
    info.set_text("t = %5.1f min   orbit %.2f\nlat %+6.1f°  lon %+7.1f°" % (t/60, t/T, tr_all["lat"][k], (tr_all["lon"][k] + 180) % 360 - 180))
    return []

save_gif(fig, frame, len(times), fps, "img/anim-gt-build.gif", dpi=90)
print("i = %.1f deg, alt = %.0f km, T = %.2f min, %g orbits = %.1f h; Earth turns %.1f deg meanwhile; node shift per orbit %.2f deg" % (
    inc, alt, T/60, norb, norb*T/3600, np.degrees(OMEGA_E*norb*T), np.degrees(OMEGA_E*T)))