"""Time note, part 1: the cost of a second, solar vs sidereal day, the equation of time,
and the geostationary altitude (why 24 h is the wrong period).

Run from the deck directory:   python codes/fig-time-days.py
Writes img/fig-time-cost.png, img/fig-time-solar-sidereal.png, img/fig-time-eot.png,
img/fig-time-geo.png and prints every number quoted in the note.
"""
import numpy as np
import matplotlib.pyplot as plt
from matplotlib.patches import Circle, FancyArrowPatch, Wedge
from _manstyle import CREAM, INK, MAROON, BLUE, GREEN, AMBER, PURPLE, despine
from _timeutil import (MU_E, RE, C_KMS, AU, MU_SUN, T_STELLAR, T_SIDEREAL, DAY, OMEGA_E,
                       Y_TROPICAL, jd_gregorian, calendar_from_jd)

GREY = "#8b8378"

# ================================================================== 1. what one second costs
speeds = {   # km/s
    "Earth's surface at the equator": OMEGA_E*RE,
    "geostationary satellite": np.sqrt(MU_E/42164.17),
    "International Space Station": np.sqrt(MU_E/(RE + 420)),
    "Earth or Mars around the Sun": 0.5*(np.sqrt(MU_SUN/AU) + np.sqrt(MU_SUN/(1.5237*AU))),
    "Parker Solar Probe (perihelion, 2024)": 192.0,
    "a radio signal": C_KMS,
}
cols = {"Earth's surface at the equator": GREY, "geostationary satellite": GREEN, "International Space Station": BLUE,
        "Earth or Mars around the Sun": MAROON, "Parker Solar Probe (perihelion, 2024)": PURPLE,
        "a radio signal": INK}
print("SPEEDS km/s: earth %.3f mars %.3f" % (np.sqrt(MU_SUN/AU), np.sqrt(MU_SUN/(1.5237*AU))),
      {k: round(v, 3) for k, v in speeds.items()})

SHORT = {"Earth's surface at the equator": "equator (rotation)", "geostationary satellite": "GEO satellite", "International Space Station": "ISS",
         "Earth or Mars around the Sun": "Earth or Mars", "Parker Solar Probe (perihelion, 2024)": "Parker Solar Probe"}
fig, (a1, a2) = plt.subplots(1, 2, figsize=(14.2, 5.3), gridspec_kw=dict(width_ratios=[1.5, 1], wspace=0.95))
dt = np.logspace(-9, 4, 300)
for k, v in speeds.items():
    a1.loglog(dt, v*1e3*dt, color=cols[k], lw=2.2, ls="-" if k != "a radio signal" else "--")
    if k != "a radio signal":
        a1.text(1.25e4, v*1e3*1e4, SHORT[k], color=cols[k], fontsize=9.5, va="center", ha="left", clip_on=False)
a1.text(2e-8, 3.2e1, "a radio signal", color=INK, fontsize=9.5, rotation=36, rotation_mode="anchor")
for y, s_ in ((1, "1 m"), (1e3, "1 km"), (1e4, "10 km")):
    a1.axhline(y, color=INK, lw=0.7, ls=":")
    a1.text(1.3e-9, y*1.25, s_, fontsize=9, color=INK)
a1.axhline(3.3895e6, color=INK, lw=0.7, ls=":"); a1.text(1.3e-9, 3.3895e6*1.25, "Mars radius (3390 km)", fontsize=9)
for x, s_ in ((1e-9, "1 ns"), (1e-6, "1 µs"), (1e-3, "1 ms"), (1, "1 s"), (60, "1 min"), (3600, "1 h")):
    a1.axvline(x, color=INK, lw=0.5, ls=":", alpha=0.6)
    a1.text(x*1.15, 2.2e9, s_, fontsize=9, rotation=90, va="top")
a1.set_xlim(1e-9, 1e4); a1.set_ylim(1e-3, 3e9)
a1.set_xlabel("timing error $\\delta t$ (s)"); a1.set_ylabel("position error  $v\\,\\delta t$  (m)")
a1.set_title("Position error = speed × timing error", fontsize=12); despine(a1)

names = list(speeds)[::-1]
vals = [speeds[k] for k in names]
a2.barh(range(len(names)), vals, color=[cols[k] for k in names], height=0.62)
a2.set_xscale("log"); a2.set_xlim(0.1, 2e6)
for j, (k, v) in enumerate(zip(names, vals)):
    a2.text(v*1.25, j, f"{v:,.3g} km" if v < 1000 else f"{v:,.0f} km", va="center", fontsize=9.5)
a2.set_yticks(range(len(names))); a2.set_yticklabels(names, fontsize=9.3)
a2.set_xlabel("distance covered in one second (km)")
a2.set_title("What one second is worth", fontsize=12); despine(a2)
fig.tight_layout(); fig.savefig("img/fig-time-cost.png", dpi=160, bbox_inches="tight"); plt.close(fig)

# ================================================================== 2. solar vs sidereal day, geometry
DELTA = np.radians(25)          # exaggerated: the true value is 0.9856 deg
fig = plt.figure(figsize=(15.4, 4.9))
gs = fig.add_gridspec(1, 4, width_ratios=[1.2, 1, 1, 1], wspace=0.04)

# --- panel 0: the orbit
ax = fig.add_subplot(gs[0]); ax.set_aspect("equal"); ax.axis("off")
ax.add_patch(Circle((0, 0), 0.20, color="#f5b301", ec=AMBER, lw=1.2, zorder=3))
ax.text(0, 0, "Sun", ha="center", va="center", fontsize=9.5, zorder=4)
th = np.linspace(0, 2*np.pi, 400); R0 = 1.0
ax.plot(R0*np.cos(th), R0*np.sin(th), color=GREY, lw=0.8, ls="--")
A0 = np.radians(200)
P0 = R0*np.array([np.cos(A0), np.sin(A0)])
ang1 = A0 + DELTA
P1 = R0*np.array([np.cos(ang1), np.sin(ang1)])
STAR = np.degrees(A0) - 180           # a star that lines up with the Sun on day 0 (direction 20 deg)
u_star = np.array([np.cos(np.radians(STAR)), np.sin(np.radians(STAR))])
for P, col in ((P0, BLUE), (P1, GREEN)):
    ax.add_patch(Circle(P, 0.085, color="#7fb3e6", ec=col, lw=1.5, zorder=4))
    ax.plot([P[0], 0], [P[1], 0], color=col, lw=1.0, zorder=2)
    ax.annotate("", xy=P + 0.55*u_star, xytext=P, arrowprops=dict(arrowstyle="-|>", color=PURPLE, lw=1.3, ls="--"))
ax.text(P0[0] - 0.14, P0[1] - 0.02, "Earth, day 0", color=BLUE, ha="right", va="center", fontsize=9.5)
ax.text(P1[0] - 0.12, P1[1] - 0.13, "Earth, one\nsidereal day later", color=GREEN, ha="right", va="top", fontsize=9.5)
ax.text(-1.02, 0.20, "to a distant star: same direction", color=PURPLE, fontsize=8.6, ha="left", va="bottom")
ax.add_patch(Wedge((0, 0), 0.55, np.degrees(A0), np.degrees(ang1), color=AMBER, alpha=0.35, zorder=1))
mid = A0 + DELTA/2
ax.text(0.72*np.cos(mid), 0.72*np.sin(mid), r"$\Delta$", color=AMBER, fontsize=15, ha="center", va="center")
ax.text(0.0, -1.62, "Orbital motion in one day:\n$\\Delta=0.9856^\\circ$  (drawn 25 times too big)", ha="center", fontsize=9.6)
ax.set_xlim(-1.45, 1.7); ax.set_ylim(-1.85, 1.3)
ax.set_title("Top view from the north ecliptic pole", fontsize=10.5)

# --- panels 1-3: the Earth disc with an observer's flag
def earth_panel(ax, flag_deg, sun_deg, star_deg, title, sub, sun_lab="Sun", star_lab="star"):
    ax.set_aspect("equal"); ax.axis("off")
    ax.add_patch(Circle((0, 0), 1, color="#cfe3f5", ec=INK, lw=1.4))
    ax.add_patch(Wedge((0, 0), 1, sun_deg + 90, sun_deg + 270, color="#31435f", alpha=0.55))   # night side
    f = np.radians(flag_deg)
    ax.plot([np.cos(f), 1.30*np.cos(f)], [np.sin(f), 1.30*np.sin(f)], color=MAROON, lw=2.8, zorder=6)
    ax.plot(np.cos(f), np.sin(f), "o", color=MAROON, ms=8, zorder=7)
    ax.text(1.55*np.cos(f) - 0.05, 1.55*np.sin(f) - 0.3, "you", color=MAROON, ha="center", va="center", fontsize=10.5)
    for d, c, lab, ls, off in ((sun_deg, AMBER, sun_lab, "-", 0.0), (star_deg, PURPLE, star_lab, "--", 0.0)):
        r = np.radians(d)
        ax.annotate("", xy=(2.3*np.cos(r), 2.3*np.sin(r)), xytext=(1.03*np.cos(r), 1.03*np.sin(r)),
                    arrowprops=dict(arrowstyle="-|>", color=c, lw=2.0, ls=ls))
        if lab:
            ax.text(2.55*np.cos(r), 2.55*np.sin(r), lab, color=c, ha="center", va="center", fontsize=10.5)
    ax.text(0, -1.5, sub, ha="center", va="top", fontsize=9.6)
    ax.set_xlim(-1.6, 3.0); ax.set_ylim(-2.6, 3.0)
    ax.set_title(title, fontsize=10.5)

sun0, sun1 = STAR, np.degrees(ang1) - 180
earth_panel(fig.add_subplot(gs[1]), sun0, sun0, sun0, "Day 0: Sun and star\ncross your meridian together",
            "noon, and the star transits", sun_lab="Sun and star", star_lab=None)
ax_ = fig.axes[-1]
earth_panel(fig.add_subplot(gs[2]), STAR + 360, sun1, STAR, "One sidereal day later: Earth has\nturned exactly $360^\\circ$",
            "the star transits again, but the Sun\nis still $\\Delta$ short of your meridian")
earth_panel(fig.add_subplot(gs[3]), sun1 + 360, sun1, STAR, "One solar day later: Earth has\nturned $360^\\circ+\\Delta$",
            "the Sun is on your meridian again:\nnoon to noon")
fig.savefig("img/fig-time-solar-sidereal.png", dpi=160, bbox_inches="tight"); plt.close(fig)

# extra time needed
extra = (DAY - T_SIDEREAL)
print("SOLAR - SIDEREAL = %.4f s = %d m %.3f s" % (extra, extra//60, extra % 60))
print("orbital motion per day (deg): tropical %.5f  sidereal %.5f" % (360/Y_TROPICAL, 360/365.256363))
print("rotation per solar day wrt equinox: %.6f deg; extra rotation to bring the Sun back: %.4f deg" %
      (360*DAY/T_SIDEREAL, 360*DAY/T_SIDEREAL - 360))
print("1/T_sid = 1/T_solar + 1/Y_tropical ->  T_sid = %.5f s" % (DAY/(1 + 1/Y_TROPICAL)))
print("stellar day %.5f s ;  sidereal day %.5f s ; difference %.2f ms" % (T_STELLAR, T_SIDEREAL, 1e3*(T_SIDEREAL - T_STELLAR)))
print("omega_E = %.10e rad/s = %.7f deg/day(solar) = %.7f deg/h" % (OMEGA_E, np.degrees(OMEGA_E)*DAY, np.degrees(OMEGA_E)*3600))
print("equatorial surface speed %.4f km/s" % (OMEGA_E*RE))

# ================================================================== 3. equation of time
def sun_low(jd):
    n = np.asarray(jd, float) - 2451545.0
    L = np.radians(280.460 + 0.9856474*n)
    g = np.radians(357.528 + 0.9856003*n)
    C = np.radians(1.915)*np.sin(g) + np.radians(0.020)*np.sin(2*g)
    lam = L + C
    eps = np.radians(23.439 - 4.0e-7*n)
    alpha = np.arctan2(np.cos(eps)*np.sin(lam), np.cos(lam))
    delta = np.arcsin(np.sin(eps)*np.sin(lam))
    wrap = lambda x: (x + np.pi) % (2*np.pi) - np.pi
    return dict(eot=np.degrees(wrap(L - alpha)), ecc=-np.degrees(C), obl=np.degrees(wrap(lam - alpha)), dec=np.degrees(delta))

jd0 = jd_gregorian(2026, 1, 1)
days = np.arange(0, 366, 0.25)
S = sun_low(jd0 + days)
eot_min = S["eot"]*4                                   # 1 deg of hour angle = 4 min
dlen = -np.gradient(S["eot"]*240, days)                # apparent day - 86400 s: a sundial that gains (EoT rising) has SHORTER days
fig, (b1, b2, b3) = plt.subplots(1, 3, figsize=(15.4, 4.5), gridspec_kw=dict(width_ratios=[1.35, 1.1, 0.9]))
mstart = np.cumsum([0, 31, 28, 31, 30, 31, 30, 31, 31, 30, 31, 30])
mlab = list("JFMAMJJASOND")
MON = ["Jan", "Feb", "Mar", "Apr", "May", "Jun", "Jul", "Aug", "Sep", "Oct", "Nov", "Dec"]
b1.plot(days, eot_min, color=INK, lw=2.6, label="equation of time")
b1.plot(days, S["ecc"]*4, color=BLUE, lw=1.4, ls="--", label="Earth's elliptic orbit (Kepler II)")
b1.plot(days, S["obl"]*4, color=GREEN, lw=1.4, ls="--", label="tilt of the ecliptic")
b1.axhline(0, color=INK, lw=0.6)
b1.set_xticks(mstart + 15); b1.set_xticklabels(mlab); b1.set_xlim(0, 365)
for x in mstart: b1.axvline(x, color=INK, lw=0.3, alpha=0.4)
b1.set_ylabel("sundial time − clock time (min)"); b1.legend(frameon=False, fontsize=8.6, loc="upper left"); despine(b1)
b1.set_title("Equation of time, 2026", fontsize=11.5)
imax, imin = np.argmax(eot_min), np.argmin(eot_min)
b1.annotate("%+.1f min" % eot_min[imax], (days[imax], eot_min[imax]), xytext=(days[imax] - 45, eot_min[imax] + 1.2), fontsize=9, arrowprops=dict(arrowstyle="-", lw=0.6))
b1.annotate("%+.1f min" % eot_min[imin], (days[imin], eot_min[imin]), xytext=(days[imin] + 22, eot_min[imin] - 0.4), fontsize=9, arrowprops=dict(arrowstyle="-", lw=0.6))
print("EoT extrema (min): max %+.2f on day %.1f (%s), min %+.2f on day %.1f (%s)" %
      (eot_min[imax], days[imax], calendar_from_jd(jd0 + days[imax]), eot_min[imin], days[imin], calendar_from_jd(jd0 + days[imin])))

b2.plot(days, 86400 + dlen, color=INK, lw=2.4)
b2.axhline(86400, color=MAROON, lw=1.0, ls="--"); b2.text(4, 86400.8, "mean solar day = 86 400 s", color=MAROON, fontsize=9)
b2.set_xticks(mstart + 15); b2.set_xticklabels(mlab); b2.set_xlim(0, 365)
b2.set_ylabel("length of the apparent solar day (s)"); despine(b2)
b2.set_title("Sundial days are not all equal", fontsize=11.5)
jl, js = np.argmax(dlen), np.argmin(dlen)
b2.annotate("longest: 86 400 %+.1f s" % dlen[jl], (days[jl], 86400 + dlen[jl]), xytext=(days[jl] - 175, 86400 + dlen[jl] - 12), fontsize=9, arrowprops=dict(arrowstyle="-", lw=0.6))
b2.annotate("shortest: 86 400 %+.1f s" % dlen[js], (days[js], 86400 + dlen[js]), xytext=(days[js] - 200, 86400 + dlen[js] - 1.5), fontsize=9, arrowprops=dict(arrowstyle="-", lw=0.6))
print("apparent solar day: longest %.2f s (day %.0f), shortest %.2f s (day %.0f); mean of the year %.4f s" %
      (86400 + dlen[jl], days[jl], 86400 + dlen[js], days[js], 86400 + dlen[:-1].mean()))

b3.plot(eot_min, S["dec"], color=INK, lw=2.4)
for k, m0 in enumerate(mstart):
    j = int(m0*4)
    b3.plot(eot_min[j], S["dec"][j], "o", color=AMBER, ms=6, zorder=5)
    b3.text(eot_min[j] + (0.8 if eot_min[j] > -8 else -0.8), S["dec"][j] + 0.6, MON[k], fontsize=8.3, ha="left" if eot_min[j] > -8 else "right")
b3.axvline(0, color=INK, lw=0.5); b3.axhline(0, color=INK, lw=0.5)
b3.set_xlabel("equation of time (min)"); b3.set_ylabel("solar declination (deg)")
b3.set_title("The analemma", fontsize=11.5); despine(b3)
fig.tight_layout(); fig.savefig("img/fig-time-eot.png", dpi=160, bbox_inches="tight"); plt.close(fig)

# ================================================================== 4. geostationary altitude
def a_of_T(T):
    return (MU_E*(T/(2*np.pi))**2)**(1/3)
aG, a24 = a_of_T(T_STELLAR), a_of_T(DAY)
print("GEO radius for T = stellar day:  a = %.3f km, altitude %.3f km" % (aG, aG - RE))
print("GEO radius for T = 24 h:         a = %.3f km, altitude %.3f km   (higher by %.2f km)" % (a24, a24 - RE, a24 - aG))
print("period at a=42164.17: %.4f s" % (2*np.pi*np.sqrt(aG**3/MU_E)))
drift = -360*(DAY/T_STELLAR - 1)       # deg per solar day, west negative
print("a 24 h-period satellite drifts %.4f deg/day (%.1f deg in 30 d, 365.25 d -> %.1f deg)" % (drift, 30*drift, 365.25*drift))
print("orbit speed GEO %.4f km/s ; 1 deg of longitude on GEO = %.1f km" % (np.sqrt(MU_E/aG), aG*np.pi/180))

fig, (g1, g2) = plt.subplots(1, 2, figsize=(12.6, 4.5))
h = np.linspace(35300, 36400, 400)
T = 2*np.pi*np.sqrt((RE + h)**3/MU_E)
g1.plot(h, T, color=INK, lw=2.4)
g1.axhline(DAY, color=MAROON, ls="--", lw=1.2); g1.text(36430, DAY + 10, "24 h = 86 400 s", color=MAROON, fontsize=10, ha="right", va="bottom")
g1.axhline(T_STELLAR, color=GREEN, ls="--", lw=1.2); g1.text(36430, T_STELLAR + 10, "one rotation of the Earth = 86 164.1 s\n(23 h 56 m 4.1 s)", color=GREEN, fontsize=10, ha="right", va="bottom")
g1.plot([aG - RE], [T_STELLAR], "o", color=GREEN, ms=9); g1.plot([a24 - RE], [DAY], "o", color=MAROON, ms=9)
g1.annotate("", xy=(a24 - RE, 86140), xytext=(aG - RE, 86140), arrowprops=dict(arrowstyle="<->", color=INK, lw=1.2))
g1.text(0.5*(a24 + aG) - RE, 86105, "%.0f km" % (a24 - aG), ha="center", fontsize=10)
g1.text(aG - RE - 12, T_STELLAR + 45, "%.0f km" % (aG - RE), ha="right", color=GREEN, fontsize=10)
g1.text(a24 - RE - 12, DAY + 45, "%.0f km" % (a24 - RE), ha="right", color=MAROON, fontsize=10)
g1.set_ylim(86050, 86560); g1.set_xlabel("altitude (km)"); g1.set_ylabel("orbital period (s)")
g1.set_title("A 24-hour orbit is not a geostationary orbit", fontsize=11.5); despine(g1)

tt = np.linspace(0, 60, 400)
lon0 = 83.0
g2.plot(tt, np.full_like(tt, lon0), color=GREEN, lw=2.6, label="$T=86\\,164.1$ s: stays at $83^\\circ$E")
g2.plot(tt, lon0 + drift*tt, color=MAROON, lw=2.6, label="$T=86\\,400$ s: slides west $%.3f^\\circ$/day" % abs(drift))
g2.set_xlabel("days after launch"); g2.set_ylabel("longitude of the sub-satellite point (deg E)")
g2.text(30, lon0 + drift*30 + 4, "$-59^\\circ$ after 60 days;\na full circle in one year", color=MAROON, fontsize=10)
g2.legend(frameon=False, fontsize=9.5, loc="lower left"); despine(g2)
g2.set_title("What the 77 km error does", fontsize=11.5)
fig.tight_layout(); fig.savefig("img/fig-time-geo.png", dpi=160, bbox_inches="tight"); plt.close(fig)