"""The leap-second essay: every number it cites from the IERS EOP 20 C04 series, computed here.
Input: sources/eopc04.1962-now (downloaded 2026-10-05 from https://hpiers.obspm.fr/iers/eop/eopc04/eopc04.1962-now).
Columns used: YR MM DD, UT1-UTC (s), LOD (s) = excess length of day over 86,400 SI seconds.
Run: python lod_c7655.py > lod_c7655.out.txt"""
import collections
import os

HERE = os.path.dirname(os.path.abspath(__file__))
rows = []
for line in open(os.path.join(HERE, "sources", "eopc04.1962-now"), encoding="utf-8"):
    if line.startswith("#") or not line.strip():
        continue
    p = line.split()
    rows.append((int(p[0]), int(p[1]), int(p[2]), float(p[7]), float(p[12])))   # y, m, d, ut1_utc, lod
print("rows", len(rows), "first", rows[0][:3], "last", rows[-1][:3])

# 1. annual mean excess LOD (ms) and days shorter than 86,400 s
by_year = collections.defaultdict(list)
for y, m, d, u, lod in rows:
    by_year[y].append(lod)
print("\n# year  mean_excess_LOD_ms  days_LOD<0  n_days  accumulated_s(sum LOD)")
for y in sorted(by_year):
    v = by_year[y]
    if y >= 1972 and (y % 4 == 0 or y >= 2010):
        print(y, round(1000 * sum(v) / len(v), 3), sum(1 for x in v if x < 0), len(v), round(sum(v), 3))

# 2. the ten shortest days in the record
print("\n# ten shortest days (excess LOD, ms)")
for y, m, d, u, lod in sorted(rows, key=lambda r: r[4])[:10]:
    print("%04d-%02d-%02d %.4f" % (y, m, d, 1000 * lod))

# 3. first year ever with a negative annual-mean excess LOD
neg = [y for y in sorted(by_year) if sum(by_year[y]) < 0]
print("\n# years with negative annual mean excess LOD:", neg)

# 4. UT1-UTC path since the last leap second (2017-01-01)
print("\n# UT1-UTC on 1 Jan / 1 Jul since 2016-12-31")
for y, m, d, u, lod in rows:
    if (y, m, d) == (2016, 12, 31) or (y >= 2017 and d == 1 and m in (1, 7)):
        print("%04d-%02d-%02d %+.4f" % (y, m, d, u))
print("last row: %04d-%02d-%02d UT1-UTC %+.7f" % (rows[-1][0], rows[-1][1], rows[-1][2], rows[-1][3]))
post = [r for r in rows if (r[0], r[1], r[2]) >= (2017, 1, 1)]
lo = min(post, key=lambda r: r[3])
print("minimum UT1-UTC since 2017-01-01: %04d-%02d-%02d %+.4f" % (lo[0], lo[1], lo[2], lo[3]))

# 5. leap-second cadence context: the longest gap between leap seconds (from sources/Leap_Second.dat)
mjd = []
for line in open(os.path.join(HERE, "sources", "Leap_Second.dat"), encoding="utf-8"):
    if line.strip() and not line.lstrip().startswith("#"):
        p = line.split()
        mjd.append((float(p[0]), int(p[3]), int(p[2]), int(p[1]), int(p[4])))
print("\n# leap-second table rows", len(mjd), "first", mjd[0][1:], "last", mjd[-1][1:])
gaps = [(mjd[i + 1][0] - mjd[i][0], mjd[i][1:4], mjd[i + 1][1:4]) for i in range(len(mjd) - 1)]
g = max(gaps)
print("longest gap in table: %.0f days, %s -> %s" % g)
print("gap since last entry to 2026-10-05 (MJD 61318): %.0f days" % (61318 - mjd[-1][0]))

# 6. like-for-like: Jan 1 - Sep 5 of every year (the seasonal cycle makes a partial year incomparable to a full one)
print("\n# Jan 1 - Sep 5 mean excess LOD (ms), days LOD<0, and UT1-UTC change over the window (s)")
for y in range(2016, 2027):
    w = [r for r in rows if r[0] == y and (r[1], r[2]) <= (9, 5)]
    print(y, round(1000 * sum(r[4] for r in w) / len(w), 3), sum(1 for r in w if r[4] < 0), len(w),
          "%+.4f" % (w[-1][3] - w[0][3]))

# 7. the shortest day of each year since 2020
print("\n# shortest day per year (ms)")
for y in range(2020, 2027):
    r = min((r for r in rows if r[0] == y), key=lambda r: r[4])
    print("%04d-%02d-%02d %.4f" % (r[0], r[1], r[2], 1000 * r[4]))

# 8. UT1-UTC peak since the leap second's recovery and its date
pk = max((r for r in rows if r[0] >= 2023), key=lambda r: r[3])
print("\nmax UT1-UTC since 2023: %04d-%02d-%02d %+.4f" % (pk[0], pk[1], pk[2], pk[3]))
print("2025-09-05 UT1-UTC: %+.4f" % next(r[3] for r in rows if r[:3] == (2025, 9, 5)))