#!/usr/bin/env python3
"""Dataset for the true-scale mechanics scene: full Moon to solar eclipse.

The widget this feeds draws NOTHING to a false scale. Every body is its real
size at its real distance, and the camera does the work of getting the reader
from a view of the Earth's orbit down to Erie's own sky. So this file's one job
is to hand over positions -- in kilometres, in one frame, for one fortnight --
plus the handful of derived lines the picture needs: the two orbits, the plane
of the ecliptic, the line of nodes, the shadow, and the four contacts.

THE FORTNIGHT. It runs from the full Moon of 25 March 2024 (the night Jeremy
McGarry filmed his "height match") to fourth contact at Erie on 8 April. Half
an orbit of the Moon: it starts on the far side of the Earth from the Sun and
swings round to sit between them. That half-turn is the whole answer to "how
did it come from behind" -- the Moon is on a narrower, faster orbit than the
one the Earth is on, and is carried round the Earth's night side onto its day
side in a fortnight while the Sun, seen from the Earth, moves fourteen degrees.

FRAME. Geocentric inertial (GCRS), kilometres. The Earth SPINS inside it: each
frame carries the exact ITRS-to-GCRS rotation as nine row-major numbers, so
anything on the surface hangs off that and the Moon, Sun, shadow and orbits do
not turn. The Sun-centred opening is the same frame seen from a camera locked
to the Sun's position -- there is no second coordinate system.

TWO TIME AXES, ONE CLOCK. Everything is keyed to hours since 25 March 00:00 UT.
The coarse track samples the fortnight every thirty minutes (seven degrees of
Earth-turn between samples, so straight-line interpolation of the rotation
matrix is within a quarter of a per cent of orthonormal, and the Earth is
eleven pixels across wherever those frames are drawn); the fine track samples
eclipse afternoon every minute and carries the shadow and the view from Erie.

THE MOON'S ORBIT is not drawn as an ellipse. It is the Moon's own geocentric
track for the 27.3 days centred on the eclipse, so the loop the reader sees is
the loop the Moon actually made. Its plane is fitted to that track and the
nodes are where the track crosses the ecliptic -- found, not assumed. That the
line of nodes points at the Sun on 8 April is then a result, and the test suite
treats it as one.

    python3 scripts/build_eclipse_mechanics.py

Writes docs/eclipse-trajectory/data/eclipse-mechanics.json.
"""

import csv
import json
import os
import sys
from datetime import datetime

import numpy as np
from skyfield import almanac
from skyfield.api import load, wgs84
from skyfield.framelib import ecliptic_frame, itrs

ROOT = os.path.dirname(os.path.dirname(os.path.abspath(__file__)))
sys.path.insert(0, os.path.join(ROOT, 'scripts'))
import mcgarry_measured_vs_model as M      # noqa: E402  the plate solution

OUT = os.path.join(ROOT, 'docs', 'eclipse-trajectory', 'data', 'eclipse-mechanics.json')
ECLIPSE_CSV = os.path.join(ROOT, 'docs', 'eclipse-trajectory', 'data',
                           'erie-eclipse-day-pixels.csv')
# His frame at totality: camera clock 14:36:11, which the eclipse's own 41.4
# minute offset puts at 15:17:35 EDT. A raw 1280x720 screen-recording frame;
# the camera image sits inside CROP, the rest is the player's chrome.
HIS_FRAME = os.path.join(ROOT, 'docs', 'eclipse-trajectory', 'img', 'his-frame-totality.jpg')
HIS_FRAME_CAM = (14, 36, 11)
HIS_CROP = (30, 30, 1250, 645)

LAT, LON = 42 + 8 / 60.0, -(80 + 5 / 60.0)      # Erie, off his own slide
R_EARTH, R_MOON, R_SUN = 6371.0, 1737.4, 696000.0
TZ = -4                                          # EDT
CLOCK_APR = 41.4      # his April clock offset, the one the eclipse itself measures

T0 = (2024, 3, 25, 0, 0)          # the fortnight starts here, 00:00 UT
COARSE_STEP_H = 0.5                # thirty minutes
COARSE_END_H = 14 * 24 + 22        # 8 April 22:00 UT
FINE_FROM_UT, FINE_TO_UT = 16.0, 21.0     # 8 April, 12:00 to 17:00 EDT
FINE_STEP_MIN = 1.0
ORBIT_DAYS = 27.32                 # one sidereal month, centred on the eclipse
EARTH_ORBIT_DAYS = 366

SITES = [('Mazatlan', 23.22, -106.42), ('Dallas', 32.78, -96.80),
         ('Indianapolis', 39.77, -86.16), ('Erie', 42.13, -80.08),
         ('Montreal', 45.50, -73.57)]

ts = load.timescale()
eph = load(os.path.join(ROOT, 'de421.bsp'))
EARTH, SUN, MOON = eph['earth'], eph['sun'], eph['moon']
SITE = EARTH + wgs84.latlon(LAT, LON)
T0_TT = ts.utc(*T0)


def t_h(h):
    """hours since T0 -> Skyfield time"""
    return ts.utc(T0[0], T0[1], T0[2], T0[3] + h, T0[4])


def hours_of(t):
    return float((t.tt - T0_TT.tt) * 24.0)


def vecs(t):
    """Geocentric Sun and Moon in GCRS km, light-time corrected."""
    return (np.asarray(EARTH.at(t).observe(SUN).position.km, float),
            np.asarray(EARTH.at(t).observe(MOON).position.km, float))


def elongation(s, m):
    """Sun-Earth-Moon angle, degrees: 180 at full Moon, 0 at conjunction."""
    return float(np.degrees(np.arccos(np.clip(float(s @ m) / np.linalg.norm(s) / np.linalg.norm(m), -1, 1))))


def rot_itrs_to_gcrs(t):
    return np.asarray(itrs.rotation_at(t), float).T


def shadow(s, m):
    """Umbra axis, apex distance from the Moon, and where the cone meets the Earth."""
    axis = m - s
    axis = axis / np.linalg.norm(axis)
    apex = R_MOON * np.linalg.norm(m - s) / (R_SUN - R_MOON)
    b = 2.0 * float(m @ axis)
    c = float(m @ m) - R_EARTH ** 2
    disc = b * b - 4 * c
    hit = None
    if disc > 0:
        k = (-b - np.sqrt(disc)) / 2.0
        if k > 0:
            p = m + k * axis
            r = float(R_MOON * (1 - k / apex))
            n = p / np.linalg.norm(p)
            cosi = float(-(axis @ n))
            hit = {'k': float(k), 'p': [round(float(x), 1) for x in p],
                   'radius': round(r, 2),
                   'incidence': round(float(np.degrees(np.arccos(np.clip(cosi, -1, 1)))), 2),
                   'footprint': round(r / max(cosi, 1e-3), 1)}
    return axis, apex, hit


def latlon_of_itrs(p):
    x, y, z = p
    return (float(np.degrees(np.arcsin(z / np.linalg.norm(p)))),
            float(np.degrees(np.arctan2(y, x))))


def erie_view(t):
    """What Erie sees: both bodies' alt/az, apparent radii, separation and the
    position angle of the Moon from the Sun (north through east)."""
    sv = SITE.at(t).observe(SUN).apparent()
    mv = SITE.at(t).observe(MOON).apparent()
    sa, sz = sv.altaz()[0].degrees, sv.altaz()[1].degrees
    ma, mz = mv.altaz()[0].degrees, mv.altaz()[1].degrees
    rs = float(np.degrees(np.arcsin(R_SUN / sv.distance().km)))
    rm = float(np.degrees(np.arcsin(R_MOON / mv.distance().km)))
    # position angle on the sky, from the Sun to the Moon: 0 = north, 90 = east
    dra = np.radians(mv.radec()[0]._degrees - sv.radec()[0]._degrees)
    dd = np.radians(mv.radec()[1].degrees - sv.radec()[1].degrees)
    pa = float(np.degrees(np.arctan2(np.sin(dra) * np.cos(np.radians(mv.radec()[1].degrees)),
                                     np.sin(dd)))) % 360
    sep = float(sv.separation_from(mv).degrees)
    # the fraction of the Sun's disc the Moon covers: two-circle lens area
    obsc = obscuration(sep, rs, rm)
    return {'sun_alt': round(float(sa), 3), 'sun_az': round(float(sz), 3),
            'obsc': round(obsc, 4),
            'moon_alt': round(float(ma), 3), 'moon_az': round(float(mz), 3),
            'r_sun': round(rs, 4), 'r_moon': round(rm, 4),
            'sep': round(sep, 4), 'pa': round(pa, 2)}


def obscuration(d, rs, rm):
    """Fraction of the Sun's disc (radius rs) covered by the Moon's (rm) at
    centre separation d. Plain plane geometry; the angles are small."""
    if d >= rs + rm:
        return 0.0
    if d <= abs(rm - rs):
        return 1.0 if rm >= rs else (rm / rs) ** 2
    a = (d * d + rs * rs - rm * rm) / (2 * d)
    b = d - a
    area = (rs * rs * np.arccos(np.clip(a / rs, -1, 1)) - a * np.sqrt(max(rs * rs - a * a, 0))
            + rm * rm * np.arccos(np.clip(b / rm, -1, 1)) - b * np.sqrt(max(rm * rm - b * b, 0)))
    return float(min(1.0, area / (np.pi * rs * rs)))


def label_ut(t):
    d = t.utc_datetime()
    return d.strftime('%-d %b %H:%M UT')


def label_edt(t):
    d = t.utc_datetime()
    m = (d.hour + TZ) * 60 + d.minute
    return f'{m // 60:02d}:{m % 60:02d} EDT'


def contacts():
    """The four contacts and totality at Erie, from apparent radii, to the second."""
    def state(sec):
        t = ts.utc(2024, 4, 8, 0, 0, sec)
        v = erie_view(t)
        return v['sep'], v['r_sun'], v['r_moon']

    def cross(lo, hi, fn):
        for _ in range(40):
            mid = (lo + hi) / 2
            if fn(mid) > 0:
                lo = mid
            else:
                hi = mid
        return (lo + hi) / 2

    ext = lambda s: state(s)[0] - (state(s)[1] + state(s)[2])        # noqa: E731
    inn = lambda s: state(s)[0] - abs(state(s)[2] - state(s)[1])     # noqa: E731
    g_lo, g_hi = 17 * 3600, 21.5 * 3600
    secs = np.arange(g_lo, g_hi, 30.0)
    seps = np.array([state(s)[0] for s in secs])
    i = int(np.argmin(seps))
    # greatest eclipse: minimum separation, by golden section on a smooth bowl
    a, b = secs[max(i - 2, 0)], secs[min(i + 2, len(secs) - 1)]
    for _ in range(60):
        m1, m2 = a + (b - a) * 0.382, a + (b - a) * 0.618
        if state(m1)[0] < state(m2)[0]:
            b = m2
        else:
            a = m1
    greatest = (a + b) / 2
    c1 = cross(greatest - 2 * 3600, greatest, ext)
    c4 = cross(greatest, greatest + 2 * 3600, lambda s: -ext(s))
    c2 = cross(greatest - 600, greatest, inn)
    c3 = cross(greatest, greatest + 600, lambda s: -inn(s))
    out = {}
    for k, s in (('c1', c1), ('c2', c2), ('greatest', greatest), ('c3', c3), ('c4', c4)):
        t = ts.utc(2024, 4, 8, 0, 0, s)
        out[k] = {'h': round(hours_of(t), 5), 'edt': label_edt(t),
                  'ut_sec': round(float(s), 1)}
    out['totality_seconds'] = round(float(c3 - c2), 1)
    return out


def his_camera_and_frame(p):
    """The plate solution the last act looks through, and where his frame at
    totality actually has the Sun, so the page can say how far apart they are.
    The solution is the March one, fitted on the Sun with the Moon held out;
    nothing is refitted for eclipse day."""
    from PIL import Image
    im = np.asarray(Image.open(HIS_FRAME).convert('RGB'), float)
    lum = im.sum(axis=2)
    x0, y0, x1, y1 = HIS_CROP
    sub = lum[y0:y1, x0:x1]
    # the eclipsed Sun is the brightest thing in the sky half of the frame; take
    # the centroid of everything within ten per cent of the peak up there,
    # rather than the single brightest pixel of a JPEG
    sky = sub[: (y1 - y0) // 2]
    ys, xs = np.where(sky >= 0.9 * sky.max())
    ix, iy = float(xs.mean()), float(ys.mean())
    hh, mm, ss = HIS_FRAME_CAM
    edt_min = hh * 60 + mm + ss / 60 + CLOCK_APR
    t = ts.utc(2024, 4, 8, 4, edt_min)
    v = erie_view(t)
    qx, qy = M.project(p, np.array([v['sun_alt']]), np.array([v['sun_az']]))
    off_px = float(np.hypot(qx[0] - (ix + x0), qy[0] - (iy + y0)))
    return ({'yaw': round(float(p[0]), 4), 'pitch': round(float(p[1]), 4), 'roll': round(float(p[2]), 4),
             'f': round(float(p[3]), 4), 'k1': round(float(p[4]), 5),
             'cx': round(float(p[5]), 3), 'cy': round(float(p[6]), 3),
             'width': 1280, 'height': 720, 'crop': list(HIS_CROP),
             'model': 'equidistant fisheye, r = f (theta + k1 theta^3), one radial term, '
                      'solved on the Sun in March with the Moon held out; nothing refitted'},
            {'file': 'img/his-frame-totality.jpg', 'camera_clock': '14:36:11',
             'h': round(hours_of(t), 5),
             'edt': f'{int(edt_min) // 60:02d}:{int(edt_min) % 60:02d}:{int(round((edt_min % 1) * 60)):02d} EDT',
             'sun_px_his': [round(ix + x0, 1), round(iy + y0, 1)],
             'sun_px_model': [round(float(qx[0]), 1), round(float(qy[0]), 1)],
             'off_px': round(off_px, 1), 'off_deg': round(float(np.degrees(off_px / p[3])), 2),
             'sun_alt': v['sun_alt'], 'sun_az': v['sun_az'], 'obsc': v['obsc']})


def measured_eclipse_day():
    """His eclipse-afternoon Sun, through the March plate solution unchanged.

    Bright blobs only, and only while the Sun is predicted to be inside the
    frame -- after that the linker is holding on to cloud. Alt/az of where his
    pixels land on the sky, so the widget can put them next to the model's Sun."""
    rs = M.rows()
    px = np.array([r['px'] for r in rs]); py = np.array([r['py'] for r in rs])
    alt, az, is_sun = M.sky(rs, M.CLOCK_MIN)
    p, kept, err = M.solve(alt, az, px, py, is_sun)
    global HIS_P
    HIS_P = p
    rows = []
    for r in csv.DictReader(open(ECLIPSE_CSV)):
        d = datetime.fromisoformat(r['raw_stamp_edt'])
        rows.append((d.day, d.hour * 60 + d.minute + d.second / 60.0,
                     float(r['pixel_x']), float(r['pixel_y']), int(r['blob_px'])))
    mins = np.array([r[1] + CLOCK_APR for r in rows])
    ex = np.array([r[2] for r in rows]); ey = np.array([r[3] for r in rows])
    npx = np.array([r[4] for r in rows])
    t = M.ts.utc(2024, 4, [r[0] for r in rows], 4, mins)
    sa, sz = M.altaz(M.SUN, t)
    qx, qy = M.project(p, sa, sz)
    ealt, eaz = M.unproject(p, ex, ey)
    ok = (npx >= 2000) & (qx < 1270)
    out = []
    for i in np.where(ok)[0]:
        tt = ts.utc(2024, 4, rows[i][0], 4, float(mins[i]))
        sep = float(np.degrees(np.arccos(np.clip(
            np.sin(np.radians(sa[i])) * np.sin(np.radians(ealt[i]))
            + np.cos(np.radians(sa[i])) * np.cos(np.radians(ealt[i]))
            * np.cos(np.radians(sz[i] - eaz[i])), -1, 1))))
        out.append({'h': round(hours_of(tt), 4), 'alt': round(float(ealt[i]), 3),
                    'az': round(float(eaz[i]), 3), 'off': round(sep, 3)})
    return out, {'n_track': len(rows), 'n_usable': int(ok.sum()),
                 'median_deg': round(float(np.median([o['off'] for o in out])), 2)}


def main():
    erie = np.asarray(wgs84.latlon(LAT, LON).itrs_xyz.km, float)
    la, lo = np.radians(LAT), np.radians(LON)
    erie_up = [float(np.cos(la) * np.cos(lo)), float(np.cos(la) * np.sin(lo)), float(np.sin(la))]

    # ---- the coarse fortnight
    coarse = []
    for h in np.arange(0.0, COARSE_END_H + 1e-9, COARSE_STEP_H):
        t = t_h(float(h))
        s, m = vecs(t)
        # the Sun-to-Earth vector IS the Earth's heliocentric position
        e_helio = np.asarray(SUN.at(t).observe(EARTH).position.km, float)
        coarse.append({'h': round(float(h), 4), 'label': label_ut(t),
                       'elong': round(elongation(s, m), 2),
                       'sun': [round(float(x), 1) for x in s],
                       'moon': [round(float(x), 1) for x in m],
                       'earth_helio': [round(float(x), 1) for x in e_helio],
                       'rot': [round(float(x), 9) for x in rot_itrs_to_gcrs(t).reshape(9)]})

    # ---- the fine afternoon
    fine = []
    day_h = 14 * 24
    for ut in np.arange(FINE_FROM_UT, FINE_TO_UT + 1e-9, FINE_STEP_MIN / 60.0):
        h = day_h + float(ut)
        t = t_h(h)
        s, m = vecs(t)
        axis, apex, hit = shadow(s, m)
        if hit:
            Rg2i = rot_itrs_to_gcrs(t).T
            pla, plo = latlon_of_itrs(Rg2i @ np.asarray(hit['p'], float))
            hit['lat'], hit['lon'] = round(pla, 3), round(plo, 3)
        fine.append({'h': round(h, 5), 'label': label_edt(t),
                     'elong': round(elongation(s, m), 3),
                     'sun': [round(float(x), 1) for x in s],
                     'moon': [round(float(x), 1) for x in m],
                     'rot': [round(float(x), 9) for x in rot_itrs_to_gcrs(t).reshape(9)],
                     'axis': [round(float(x), 7) for x in axis],
                     'apex': round(float(apex), 1),
                     'hit': hit,
                     'view': erie_view(t)})

    # ---- the Moon's own loop, one sidereal month centred on the eclipse
    con = contacts()
    hg = con['greatest']['h']
    loop = []
    for h in np.arange(hg - ORBIT_DAYS * 12, hg + ORBIT_DAYS * 12 + 1e-9, 1.0):
        t = t_h(float(h))
        _, m = vecs(t)
        loop.append((float(h), m))
    moon_orbit = [[round(float(x), 1) for x in m] for _, m in loop]
    # its plane, fitted; its tilt to the ecliptic; its nodes, found
    P = np.array([m for _, m in loop])
    cen = P.mean(axis=0)
    _, _, vt = np.linalg.svd(P - cen)
    normal = vt[2]
    ecl = np.asarray(ecliptic_frame.rotation_at(t_h(hg)), float)   # rows: ecliptic axes in ICRS
    ecl_z = ecl[2]
    if float(normal @ ecl_z) < 0:
        normal = -normal
    incl = float(np.degrees(np.arccos(np.clip(float(normal @ ecl_z), -1, 1))))
    lat_ecl = [float(np.degrees(np.arcsin((m / np.linalg.norm(m)) @ ecl_z))) for _, m in loop]
    nodes = []
    for i in range(1, len(loop)):
        if lat_ecl[i - 1] == 0 or (lat_ecl[i - 1] < 0) != (lat_ecl[i] < 0):
            f = lat_ecl[i - 1] / (lat_ecl[i - 1] - lat_ecl[i])
            h = loop[i - 1][0] + f * (loop[i][0] - loop[i - 1][0])
            _, m = vecs(t_h(h))
            d = m / np.linalg.norm(m)
            nodes.append({'h': round(h, 3), 'label': label_ut(t_h(h)),
                          'kind': 'ascending' if lat_ecl[i] > lat_ecl[i - 1] else 'descending',
                          'dir': [round(float(x), 7) for x in d]})
    # the node nearest the eclipse, and how far its line is from the Sun then
    sg, mg = vecs(t_h(hg))
    near_node = min(nodes, key=lambda n: abs(n['h'] - hg))
    nd = np.asarray(near_node['dir'], float)
    node_sun = float(np.degrees(np.arccos(abs(float(nd @ (sg / np.linalg.norm(sg)))))))

    # ---- the Earth's orbit, a year of it, heliocentric
    earth_orbit = []
    for d in range(-EARTH_ORBIT_DAYS // 2, EARTH_ORBIT_DAYS // 2 + 1):
        t = t_h(hg + 24.0 * d)
        e = np.asarray(SUN.at(t).observe(EARTH).position.km, float)
        earth_orbit.append([round(float(x), 0) for x in e])

    # ---- phases
    tf, te = t_h(-24), t_h(COARSE_END_H + 24)
    ph_t, ph_k = almanac.find_discrete(tf, te, almanac.moon_phases(eph))
    phases = [{'h': round(hours_of(t), 3), 'label': label_ut(t),
               'phase': almanac.MOON_PHASES[int(k)]} for t, k in zip(ph_t, ph_k)]

    # ---- the ground track and the five sites, as before
    track = []
    for ut in np.arange(15.0, 22.0 + 1e-9, 1 / 60.0):
        t = t_h(day_h + float(ut))
        s, m = vecs(t)
        R = rot_itrs_to_gcrs(t)
        _, _, hit = shadow(R.T @ s, R.T @ m)
        if hit and hit['radius'] > 0:
            pla, plo = latlon_of_itrs(np.asarray(hit['p'], float))
            track.append({'h': round(day_h + float(ut), 4), 'lat': round(pla, 3),
                          'lon': round(plo, 3), 'r': hit['radius']})
    # Who the umbra lands on. A site is in totality when its perpendicular
    # distance from the shadow axis is less than the cone's radius at that
    # point -- the cone's own definition, applied to the site's own position,
    # second by second. An earlier version measured each site against the
    # nearest ONE-MINUTE SAMPLE of the ground track (the umbra moves sixty
    # kilometres a minute, so up to thirty kilometres of error) after turning
    # the sphere hit into latitude/longitude geocentrically and reading it
    # back geodetically (up to twenty more at this latitude). The two errors
    # put Montreal 101 km outside an 85 km half-width. It was 77 km inside,
    # by eight kilometres, with a minute and a half of totality -- which is
    # what Montreal saw on the day.
    def cone_margin(h, v):
        t = t_h(h)
        Rg2i = rot_itrs_to_gcrs(t).T
        s, m = vecs(t)
        s, m = Rg2i @ s, Rg2i @ m
        a = (m - s) / np.linalg.norm(m - s)
        apex = R_MOON * np.linalg.norm(m - s) / (R_SUN - R_MOON)
        k = float((v - m) @ a)
        perp = float(np.linalg.norm((v - m) - k * a))
        return perp, float(R_MOON * (1 - k / apex))

    def site_pass(v):
        hs = np.arange(day_h + 17.5, day_h + 21.0, 1 / 3600)
        ms = [cone_margin(h, v) for h in hs]
        i = int(np.argmin([p - r for p, r in ms]))
        inside = [h for h, (p, r) in zip(hs, ms) if p < r]
        return {'closest_h': round(float(hs[i]), 5), 'miss_km': round(ms[i][0], 1),
                'umbra_radius_km': round(ms[i][1], 2), 'inside': bool(ms[i][0] < ms[i][1]),
                'totality_s': round((inside[-1] - inside[0]) * 3600, 0) if inside else 0}

    sites = []
    for name, sla, slo in SITES:
        v = np.asarray(wgs84.latlon(sla, slo).itrs_xyz.km, float)
        sites.append({'name': name, 'lat': sla, 'lon': slo,
                      'xyz': [round(float(x), 1) for x in v], **site_pass(v)})
    erie_pass = site_pass(erie)
    best = min(track, key=lambda q: abs(q['h'] - erie_pass['closest_h']))
    d_best = erie_pass['miss_km']

    meas, meas_stats = measured_eclipse_day()
    his_cam, his_frame = his_camera_and_frame(HIS_P)

    out = {
        't0': '2024-03-25T00:00:00Z',
        'frame': ('Geocentric inertial (GCRS), kilometres, hours since t0. The Earth '
                  'spins inside it; each frame carries the exact ITRS-to-GCRS rotation. '
                  'Nothing is scaled: every radius and distance is the true one.'),
        'site': {'label': 'Erie, Pennsylvania', 'lat': LAT, 'lon': LON,
                 'xyz': [round(float(x), 1) for x in erie], 'up': erie_up},
        'radii_km': {'earth': R_EARTH, 'moon': R_MOON, 'sun': R_SUN},
        'ecliptic': {'x': [round(float(x), 9) for x in ecl[0]],
                     'y': [round(float(x), 9) for x in ecl[1]],
                     'z': [round(float(x), 9) for x in ecl[2]]},
        'moon_orbit': {'points': moon_orbit, 'from_h': round(loop[0][0], 6),
                       'to_h': round(loop[-1][0], 6), 'step_h': 1.0,
                       'normal': [round(float(x), 7) for x in normal],
                       'inclination_deg': round(incl, 3),
                       'nodes': nodes,
                       'node_line_to_sun_deg_at_greatest': round(node_sun, 3)},
        'earth_orbit': {'points': earth_orbit, 'step_days': 1,
                        'centre_h': round(hg, 6)},
        'phases': phases,
        'contacts': con,
        'coarse': {'step_h': COARSE_STEP_H, 'frames': coarse},
        'fine': {'step_min': FINE_STEP_MIN, 'frames': fine},
        'track': track,
        'sites': sites,
        'closest': {'h': erie_pass['closest_h'], 'lat': best['lat'], 'lon': best['lon'],
                    'umbra_radius_km': erie_pass['umbra_radius_km'], 'miss_km': d_best,
                    'totality_s': erie_pass['totality_s']},
        'measured': {'eclipse_day': meas, **meas_stats},
        'his_camera': his_cam,
        'his_frame': his_frame,
    }
    with open(OUT, 'w') as fh:
        json.dump(out, fh, separators=(',', ':'))

    # ---- self-checks, printed and asserted
    print(f"wrote {os.path.relpath(OUT, ROOT)}  ({os.path.getsize(OUT) / 1024:.0f} KB)")
    print(f"  coarse {len(coarse)} frames every {COARSE_STEP_H * 60:.0f} min; "
          f"fine {len(fine)} frames every {FINE_STEP_MIN:.0f} min")
    full = [p for p in phases if p['phase'] == 'Full Moon']
    new = [p for p in phases if p['phase'] == 'New Moon']
    print(f"  phases: {', '.join(p['phase'] + ' ' + p['label'] for p in phases)}")
    s0, m0 = vecs(t_h(full[0]['h']))
    elong = float(np.degrees(np.arccos(float(s0 @ m0) / np.linalg.norm(s0) / np.linalg.norm(m0))))
    print(f"  at full Moon the Moon is {elong:.1f} deg from the Sun (opposite)")
    assert elong > 175, elong
    print(f"  Moon's orbit: tilt {incl:.2f} deg to the ecliptic; nodes "
          + ', '.join(f"{n['kind']} {n['label']}" for n in nodes))
    assert 4.9 < incl < 5.4, incl
    print(f"  at greatest eclipse the line of nodes is {node_sun:.2f} deg from the Sun's direction")
    assert node_sun < 6.0, node_sun
    # New Moon is the GEOCENTRIC conjunction; Erie's greatest eclipse is an hour
    # later because Erie is off the axis and sees the Moon shifted by parallax.
    # The moment the axis passes nearest the Earth's centre is what New Moon
    # should agree with, to within minutes.
    def _axis_miss(f):
        m, a = np.asarray(f['moon'], float), np.asarray(f['axis'], float)
        return float(np.linalg.norm(m - (m @ a) * a))
    fglob = min(fine, key=_axis_miss)
    print(f"  new Moon {new[0]['label']}; axis nearest the Earth's centre {label_ut(t_h(fglob['h']))} "
          f"({_axis_miss(fglob):.0f} km off it); greatest at Erie {con['greatest']['edt']} "
          f"({label_ut(t_h(hg))}), an hour later by parallax; totality {con['totality_seconds']:.0f} s")
    assert abs(new[0]['h'] - fglob['h']) < 0.2, (new[0]['h'], fglob['h'])
    assert 0.5 < hg - new[0]['h'] < 1.5, (hg, new[0]['h'])
    # the page counts hours from the node in two ways; both are printed so the
    # figure is traceable. Nodes here are instantaneous zero crossings of
    # ecliptic latitude on the true ecliptic of date, sampled hourly.
    asc = [n for n in nodes if n['kind'] == 'ascending' and abs(n['h'] - hg) < 48][0]
    print(f"  ascending node {asc['label']}: {fglob['h'] - asc['h']:.2f} h before the axis passes "
          f"closest to the Earth's centre ({label_ut(t_h(fglob['h']))}), {hg - asc['h']:.2f} h before greatest "
          f"eclipse at Erie ({con['greatest']['edt']})")
    assert 5.5 < fglob['h'] - asc['h'] < 6.5 and 6.5 < hg - asc['h'] < 7.5, (fglob['h'] - asc['h'], hg - asc['h'])
    out['greatest_global'] = {'h': fglob['h'], 'label': label_ut(t_h(fglob['h'])),
                              'axis_miss_km': round(_axis_miss(fglob), 1)}
    with open(OUT, 'w') as fh:
        json.dump(out, fh, separators=(',', ':'))
    ds = [np.linalg.norm(np.asarray(f['moon'], float)) for f in coarse]
    print(f"  Moon's distance over the fortnight {min(ds):.0f}-{max(ds):.0f} km")
    for f in coarse[::200]:
        R = np.asarray(f['rot'], float).reshape(3, 3)
        assert np.allclose(R @ R.T, np.eye(3), atol=1e-6)
    print("  who the umbra lands on (perpendicular distance from the axis against the cone's radius there):")
    for st in sites:
        print(f"     {st['name']:13s} {st['miss_km']:6.1f} km off the axis, cone {st['umbra_radius_km']:5.1f} km -> "
              f"{'totality ' + str(int(st['totality_s'])) + ' s' if st['inside'] else 'partial only'}")
    assert erie_pass['inside'] and abs(erie_pass['totality_s'] - con['totality_seconds']) < 5, \
        (erie_pass, con['totality_seconds'])
    mont = [st for st in sites if st['name'] == 'Montreal'][0]
    assert mont['inside'] and 60 < mont['totality_s'] < 120, mont
    ff = [f for f in fine if f['hit']]
    print(f"  umbra on the Earth in {len(ff)} of {len(fine)} fine frames; footprint "
          f"{min(f['hit']['footprint'] for f in ff):.0f}-{max(f['hit']['footprint'] for f in ff):.0f} km")
    vg = min(fine, key=lambda f: abs(f['h'] - hg))['view']
    print(f"  from Erie at greatest: Sun alt {vg['sun_alt']:.1f} az {vg['sun_az']:.1f}, "
          f"separation {vg['sep']:.3f} deg, Moon/Sun radii {vg['r_moon']:.3f}/{vg['r_sun']:.3f}")
    assert vg['r_moon'] > vg['r_sun'] and vg['sep'] < 0.02
    v1 = min(fine, key=lambda f: abs(f['h'] - con['c1']['h']))['view']
    print(f"  at first contact the Moon stands at position angle {v1['pa']:.0f} deg from the Sun "
          f"(0 north, 90 east, 180 south, 270 west)")
    print(f"  his eclipse-afternoon Sun: {meas_stats['n_usable']} of {meas_stats['n_track']} usable, "
          f"median {meas_stats['median_deg']} deg off the model")
    print(f"  his frame at totality ({his_frame['edt']}): Sun at pixel {his_frame['sun_px_his']}, "
          f"the model through his lens puts it at {his_frame['sun_px_model']} -- "
          f"{his_frame['off_px']} px, {his_frame['off_deg']} deg")
    assert his_frame['off_px'] < 20 and his_frame['obsc'] == 1.0, his_frame


if __name__ == '__main__':
    main()