#!/usr/bin/env python3 """Dataset for the animated Erie sky-track: the model's month, and his three days on it. The chart this feeds shows the same thing the Reykjavik graphic shows, for his site and his weeks: one arc per day for the Moon, the Sun's arc barely moving, and a dot per day at the eclipse's own clock time so the Moon's march across the Sun is visible as a row of dots rather than asserted in prose. What it adds is his own footage. Every measured pixel position is carried back onto the sky through the plate solution in mcgarry_measured_vs_model.py -- solved on the SUN ALONE, so the Moon points are held out and land on a curve the fit never saw. Those become dots on the model's own arcs for 23, 24 and 25 March. Two honesty notes travel in the JSON so the widget can print them rather than the page having to remember: clock_min his stated overlay correction. This footage cannot determine it (a clock error and a yaw error trade against each other at 0.30 deg/min with a nearly flat residual), so it is set to his value and declared, not fitted. solved_on 'sun'. The Moon residual is a held-out number and is the one worth quoting. python3 scripts/build_erie_skytrack.py Writes docs/eclipse-trajectory/data/erie-skytrack.json. """ import csv import json import os import sys from datetime import datetime import numpy as np sys.path.insert(0, os.path.dirname(os.path.abspath(__file__))) import mcgarry_measured_vs_model as M # noqa: E402 (path set above) from skyfield.framelib import ecliptic_frame as ECL # noqa: E402 OUT = os.path.join(M.ROOT, 'docs', 'eclipse-trajectory', 'data', 'erie-skytrack.json') # 18 March to 8 April: his reference day, his three filmed days, the Moon's dip # to its lowest arc of the month, and the climb back onto the Sun. SPAN = [(3, d) for d in range(18, 32)] + [(4, d) for d in range(1, 9)] APRIL_CSV = os.path.join(M.ROOT, 'docs', 'eclipse-trajectory', 'data', 'erie-april-pixels.csv') # The April segment's clock offset is the one the eclipse itself measures, not # the 43 minutes he states for March. Both are declared, neither is fitted here. CLOCK_APR = 41.4 ECLIPSE_CSV = os.path.join(M.ROOT, 'docs', 'eclipse-trajectory', 'data', 'erie-eclipse-day-pixels.csv') MARK_MIN = 15 * 60 + 18 # 15:18 EDT -- greatest eclipse STEP = 6 # minutes between samples along a daily arc def altaz_at(body, mo, day, minutes): """Refracted alt/az at Erie for a list of EDT minutes on one date.""" t = M.ts.utc(2024, mo, day, 4, list(minutes)) return M.altaz(body, t) def arc(body, mo, day): """One day's track, sampled while the body is up. Returns EDT minutes + sky.""" mins = list(range(0, 1440, STEP)) alt, az = altaz_at(body, mo, day, mins) # 2 dp = 36 arcsec, far finer than the chart can draw; 3 dp doubled the file return [{'m': m, 'az': round(float(z), 2), 'alt': round(float(a), 2)} for m, a, z in zip(mins, alt, az) if a > -1.0] def april_rows(): """The 7 April Sun track, same shape as the March rows. Eclipse afternoon is deliberately absent. Through the partial phases Erie is under broken cloud and no smooth track survives; the one clean post-totality track sits 1.6 deg from the Sun with 8% of its points inside a degree, out in the frame corner where the one-term lens model is extrapolating and the March fit already rejects points. It is not good enough to publish, and leaving it out is the honest choice rather than a gap. """ out = [] for r in csv.DictReader(open(APRIL_CSV)): d = datetime.fromisoformat(r['raw_stamp_edt']) out.append(dict(day=d.day, min=d.hour * 60 + d.minute + d.second / 60.0, body=r['body'], px=float(r['pixel_x']), py=float(r['pixel_y']))) return out def eclipse_day_stats(p): """What his eclipse-afternoon footage gives, and why it stays off the chart. The Sun IS recovered after totality -- this is not a null result. But the residual is systematic rather than scattered: the track runs out to the right edge of a 1280-wide frame, and the one-term radial term the lens model uses under-bends there, so the measured Sun sits a few pixels inside the predicted one and the offset grows with radius. Quoting the whole track would also fold in points taken after the Sun has left the frame entirely, where the linker is holding on to cloud, so the summary is over bright blobs while the Sun is still predicted to be inside the frame. """ 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']))) day = [r[0] for r in rows] mins = [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, day, 4, mins) sa, sz = M.altaz(M.SUN, t) qx, qy = M.project(p, sa, sz) # where the model puts the Sun ealt, eaz = M.unproject(p, ex, ey) sep = angsep(sa, sz, ealt, eaz) ok = (npx >= 2000) & (qx < 1270) # bright, and Sun still in frame rad = np.hypot(ex - p[5], ey - p[6]) return { 'n_track': len(rows), 'n_usable': int(ok.sum()), 'from_edt': round(min(mins), 1), 'to_edt': round(max(mins), 1), 'median_deg': round(float(np.median(sep[ok])), 2), 'within_1deg_pct': round(float((sep[ok] < 1.0).mean() * 100)), 'within_2deg_pct': round(float((sep[ok] < 2.0).mean() * 100)), 'leaves_frame_edt': round(float(np.array(mins)[qx >= 1270].min()), 1) if (qx >= 1270).any() else None, 'dx_inner_px': round(float(np.median((ex - qx)[ok & (rad < 500)])), 1), 'dx_outer_px': round(float(np.median((ex - qx)[ok & (rad >= 520)])), 1), 'note': ('the Sun is recovered, but the residual is systematic and grows ' 'with radius -- the one-term lens model under-bending in the ' 'frame corner, not scatter and not cloud'), } def contacts(): """The four contacts and the totality window, from apparent radii. These were constants in an earlier version of this file. They are computed now, because the chart shades a band between them and a shaded band is a claim about when the eclipse was, which ought to come from the ephemeris rather than from something typed in once and never checked again. """ def state(minute): t = M.ts.utc(2024, 4, 8, 4, minute) sv = M.SITE.at(t).observe(M.SUN).apparent() mv = M.SITE.at(t).observe(M.MOON).apparent() sa, sz = sv.altaz()[0].degrees, sv.altaz()[1].degrees ma, mz = mv.altaz()[0].degrees, mv.altaz()[1].degrees sep = np.hypot((mz - sz) * np.cos(np.radians(sa)), ma - sa) rs = np.degrees(np.arcsin(696000.0 / sv.distance().km)) rm = np.degrees(np.arcsin(1737.4 / mv.distance().km)) return sep, rs, rm, sz grid = np.arange(13 * 60, 17 * 60, 0.05) v = [state(x) for x in grid] partial = [g for g, q in zip(grid, v) if q[0] < q[1] + q[2]] total = [g for g, q in zip(grid, v) if q[0] < q[2] - q[1]] great = grid[int(np.argmin([q[0] for q in v]))] return { 'date': '2024-04-08', 'begin_min': round(float(partial[0]), 2), 'end_min': round(float(partial[-1]), 2), 'greatest_min': round(float(great), 2), 'total_begin_min': round(float(total[0]), 2) if total else None, 'total_end_min': round(float(total[-1]), 2) if total else None, 'total_seconds': round(float(total[-1] - total[0]) * 60) if total else 0, } def culminate(pts): """The day's meridian crossing, or None when the day has not got one. A lunar day is about 24h50m, so roughly once a month a calendar day contains no crossing at all -- 22 March 2024 at Erie is that day, sitting between crossings at 23:24 on the 21st and 00:04 on the 23rd. Taking a plain maximum over the day returns the midnight boundary sample instead, which is the Moon already an hour past the meridian and heading down. Plotted, it is a dot stranded off the column of real culminations. Better to say there is not one. """ if len(pts) < 3: return None i = max(range(len(pts)), key=lambda k: pts[k]['alt']) if i == 0 or i == len(pts) - 1: return None return pts[i] def solve_camera(): """The published plate solution: fit on the Sun, Moon held out.""" 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) malt, maz = M.unproject(p, px, py) # Residuals in the SAME convention the published figure uses: pixel error # divided by the focal length. For an equidistant fisheye r = f*theta that is # near enough exact, and using it here means the chart and the table on the # page cannot drift apart. The true great-circle separation is carried # alongside as a cross-check; the two agree to about half an arcminute. deg = err / p[3] * 57.29577951 sep = angsep(alt, az, malt, maz) return rs, p, kept, is_sun, malt, maz, deg, sep def elongation(mo, day, minutes): """Moon's elongation east of the Sun in ecliptic longitude, 0-360. The plain angular separation is NOT a progress measure: it climbs to 180 at full Moon and falls back, so it says 164 degrees on 23 March and 174 on the 25th while the Moon is in fact marching steadily the whole time. Elongation runs once round the circle per synodic month and does not fold, so 360 - elongation is a genuine countdown to new Moon -- which is the eclipse. """ t = M.ts.utc(2024, mo, day, 4, list(minutes)) lm = M.SITE.at(t).observe(M.MOON).apparent().frame_latlon(ECL)[1].degrees ls = M.SITE.at(t).observe(M.SUN).apparent().frame_latlon(ECL)[1].degrees return (lm - ls) % 360.0 def angsep(a1, z1, a2, z2): """Great-circle separation on the sky, degrees.""" a1, z1, a2, z2 = map(np.radians, (a1, z1, a2, z2)) c = np.sin(a1) * np.sin(a2) + np.cos(a1) * np.cos(a2) * np.cos(z1 - z2) return np.degrees(np.arccos(np.clip(c, -1, 1))) def main(): global ECL_CONTACTS ECL_CONTACTS = contacts() rs, p, kept, is_sun, malt, maz, deg, sep = solve_camera() # measured positions, bucketed by the date they were actually taken meas = {} for i, r in enumerate(rs): key = f"2024-03-{r['day']:02d}" b = meas.setdefault(key, {'sun': [], 'moon': []}) b[r['body']].append({'m': round(r['min'] + M.CLOCK_MIN, 2), 'az': round(float(maz[i]), 3), 'alt': round(float(malt[i]), 3), 'off': round(float(deg[i]), 3)}) for b in meas.values(): for k in b: b[k].sort(key=lambda d: d['m']) # 7 April, through the SAME solution with nothing refitted ar = april_rows() apx = np.array([r['px'] for r in ar]); apy = np.array([r['py'] for r in ar]) at = M.ts.utc(2024, 4, [r['day'] for r in ar], 4, [r['min'] + CLOCK_APR for r in ar]) asa, asz = M.altaz(M.SUN, at) aalt, aaz = M.unproject(p, apx, apy) asep = angsep(asa, asz, aalt, aaz) for i, r in enumerate(ar): b = meas.setdefault(f"2024-04-{r['day']:02d}", {'sun': [], 'moon': []}) b[r['body']].append({'m': round(r['min'] + CLOCK_APR, 2), 'az': round(float(aaz[i]), 3), 'alt': round(float(aalt[i]), 3), 'off': round(float(asep[i]), 3)}) for b in meas.values(): for k in b: b[k].sort(key=lambda q: q['m']) days = [] for mo, d in SPAN: key = f'2024-{mo:02d}-{d:02d}' malt_mark, maz_mark = altaz_at(M.MOON, mo, d, [MARK_MIN]) salt_mark, saz_mark = altaz_at(M.SUN, mo, d, [MARK_MIN]) marc = arc(M.MOON, mo, d) sarc = arc(M.SUN, mo, d) stop = max(sarc, key=lambda q: q['alt']) top = culminate(marc) elong = float(elongation(mo, d, [MARK_MIN])[0]) days.append({ 'date': key, 'label': f'{d} {"Mar" if mo == 3 else "Apr"}', 'moon': marc, 'sun': sarc, 'culmination': round(top['alt'], 2) if top else None, 'culm': {'moon': ({'az': top['az'], 'alt': top['alt'], 'm': top['m']} if top else None), 'sun': {'az': stop['az'], 'alt': stop['alt'], 'm': stop['m']}}, 'elong': round(elong, 2), 'to_go': round(360.0 - elong if elong > 0 else 0.0, 2), 'mark': { 'moon': {'az': round(float(maz_mark[0]), 3), 'alt': round(float(malt_mark[0]), 3)}, 'sun': {'az': round(float(saz_mark[0]), 3), 'alt': round(float(salt_mark[0]), 3)}, 'gap': round(float(angsep(malt_mark, maz_mark, salt_mark, saz_mark)[0]), 3), }, 'measured': meas.get(key), }) ns, nm = int(is_sun.sum()), int((~is_sun).sum()) out = { 'site': {'lat': M.LAT, 'lon': M.LON, 'label': 'Erie, Pennsylvania'}, 'mark_min': MARK_MIN, 'mark_label': '15:18 EDT', 'clock_min': M.CLOCK_MIN, 'clock_note': ('his own stated overlay correction; this footage cannot ' 'determine it, so it is declared rather than fitted'), 'camera': {n: round(float(v), 4) for n, v in zip(M.PNAME, p)}, 'fit': { 'solved_on': 'sun', 'n_sun': ns, 'n_moon': nm, 'n_sun_kept': int(kept.sum()), 'sun_median_deg': round(float(np.median(deg[kept])), 3), 'moon_median_deg': round(float(np.median(deg[~is_sun])), 3), 'moon_within_1deg_pct': round(float((deg[~is_sun] < 1.0).mean() * 100), 1), 'residual_convention': 'pixel error / focal length, as on the published figure', 'moon_median_greatcircle_deg': round(float(np.median(sep[~is_sun])), 3), }, 'april': {'clock_min': CLOCK_APR, 'n_sun': len(ar), 'median_deg': round(float(np.median(asep)), 3), 'within_1deg_pct': round(float((asep < 1.0).mean() * 100), 1), 'note': ('the March solution applied unchanged two weeks later; ' 'eclipse afternoon itself is not measurable from this ' 'footage -- cloud through the partial phases, and the ' 'one clean post-totality track sits in the frame corner ' 'where the lens model extrapolates')}, 'eclipse_day': eclipse_day_stats(p), 'eclipse': ECL_CONTACTS, 'days': days, } with open(OUT, 'w') as fh: json.dump(out, fh, separators=(',', ':')) print(f'wrote {OUT} {os.path.getsize(OUT) / 1024:.0f} KB') print(f" days {len(days)} measured days {sorted(meas)}") e = out['eclipse'] def _hm(x): return f'{int(x) // 60:02d}:{int(x) % 60:02d}' print(f" contacts: C1 {_hm(e['begin_min'])} C2 {_hm(e['total_begin_min'])} " f"greatest {_hm(e['greatest_min'])} C3 {_hm(e['total_end_min'])} " f"C4 {_hm(e['end_min'])} totality {e['total_seconds']} s") ed = out['eclipse_day'] print(f" 8 Apr post-totality: {ed['n_track']} tracked, {ed['n_usable']} usable, " f"median {ed['median_deg']} deg, {ed['within_2deg_pct']}% inside 2 deg") print(f" systematic dx {ed['dx_inner_px']} px inner -> {ed['dx_outer_px']} px outer; " f"Sun leaves frame at {int(ed['leaves_frame_edt']) // 60}:" f"{int(ed['leaves_frame_edt']) % 60:02d} EDT") print(f" 7 Apr, March solution unchanged: {len(ar)} Sun, " f"median {out['april']['median_deg']} deg, " f"{out['april']['within_1deg_pct']}% inside 1 deg") print(f" fit on sun: {ns} pts ({out['fit']['n_sun_kept']} kept), " f"median {out['fit']['sun_median_deg']} deg") print(f" held-out moon: {nm} pts, median {out['fit']['moon_median_deg']} deg, " f"{out['fit']['moon_within_1deg_pct']}% inside 1 deg " f"(great-circle cross-check {out['fit']['moon_median_greatcircle_deg']})") print(f" gap at {out['mark_label']}: " f"{days[0]['mark']['gap']:.1f} deg on {days[0]['label']} -> " f"{days[-1]['mark']['gap']:.4f} deg on {days[-1]['label']}") tg = [x['to_go'] for x in days] st = [round(tg[i] - tg[i + 1], 2) for i in range(len(tg) - 1)] print(f" still to go: {tg[0]:.1f} deg -> {tg[-1]:.2f} deg, " f"steps {min(st):.1f}-{max(st):.1f} deg/day, monotone={all(x > 0 for x in st)}") cs = [x for x in days if x['culmination'] is not None] print(f" culmination: {days[0]['culmination']:.1f} -> " f"{min(x['culmination'] for x in cs):.1f} -> {days[-1]['culmination']:.1f} " f"({len(cs)} of {len(days)} days have a meridian crossing; " f"none on {[x['label'] for x in days if x['culmination'] is None]})") azs = [abs(x['culm']['moon']['az'] - 180) for x in cs] print(f" every crossing is on the meridian: max |az-180| = {max(azs):.1f} deg") if __name__ == '__main__': main()