#!/usr/bin/env python3 """How high was the mountain the Santa Fe camera was pointed at? The p. 169 frame is a long-lens shot of the Jemez skyline toward Los Alamos with the eclipsed Moon above it. That is enough to answer a question the video never asks: could that camera have recorded the selenelion at all? The Moon is its own ruler. Its angular diameter that morning is known to four decimals from the ephemeris, so fitting a circle to the lit outer arc calibrates the frame in degrees per pixel with nothing else assumed -- no focal length, no sensor size, no claim about where exactly the tripod stood. Measure down from the Moon's centre to the skyline directly beneath it and the ridge's apparent altitude falls out. The headline result is deliberately stated in a form that does not need the frame's timestamp at all. The Moon has to fall its own centre-to-skyline gap plus one semidiameter before the disc is wholly gone, and it falls at a rate the ephemeris gives. That converts the pixel measurement straight into *how much Moon was left in the shot* -- about seven minutes -- whenever the shutter actually opened. Only then are absolute clock times attached, and they are attached from the page's own independent dating of the frame (06:47, from its eclipse phase), not from anything measured here. A sensitivity sweep across every plausible frame time is printed alongside, because that is the number a reader should check. Two things this script deliberately does NOT claim: * It does not date the frame from the crescent. The lit thickness gives the umbral magnitude in principle, but the umbra's edge is diffuse and the answer swings from 0.82 to 0.88 -- five minutes of clock -- on the choice of threshold alone. It is printed as a loose consistency check and nothing rests on it. * It does not name the mountain. The implied summit height depends on a range to the ridge that the frame cannot supply, and there are nearer ridges in this view as well as the Jemez skyline behind them. The figure is printed for orientation only. The conclusion survives both of those. Across the whole plausible range the Moon is entirely behind the skyline before 07:00, and the Sun's upper limb does not clear even a flat horizon until 07:02:23 -- with the true sunrise later still, since the Sun came up behind the Sangre de Cristos at the cameraman's back. He stood inside the band and could not have seen it from where he stood. Input frame is a still from `moon observations` on the Shape Debate channel, carrying Steve Matthews' footage; it is not redistributed here. Point --frame at your own copy. python3 santafe_ridge.py --frame path/to/frame.png """ import argparse import json import math import os import numpy as np from PIL import Image, ImageDraw HERE = os.path.dirname(os.path.abspath(__file__)) ROOT = os.path.dirname(HERE) FRAMES = os.path.join(ROOT, 'docs', 'selenelion', 'selenelion-frames.json') OUT = os.path.join(ROOT, 'docs', 'selenelion', 'img', 'fig-santafe-ridge.jpg') # Where to look in the frame. These are generous boxes, not measurements: the # crescent and the skyline are both found inside them by threshold, not by hand. MOON_BOX = (950, 1150, 270, 360) # x0, x1, y0, y1 RIDGE_X = (940, 1140) SKY_PROBE_Y = 430 # clear sky above the ridge, below the Moon OBSERVER_M = 2130.0 # vantage north of Santa Fe, metres K_REFRACT = 7.0 / 6.0 # standard terrestrial refraction R_EARTH = 6371008.8 # ---------------------------------------------------------------- measurement def find_crescent(im): """Fit a circle to the Moon's lit outer arc, and measure the lit thickness. The crescent's outer boundary IS the lunar limb -- the umbra cuts the inner side, not the outer -- so a circle fit to the lower edge recovers the whole disc even though most of it is dark. """ a = np.asarray(im.convert('RGB')).astype(float) R, B = a[..., 0], a[..., 2] x0, x1, y0, y1 = MOON_BOX box = np.zeros(R.shape, bool) box[y0:y1, x0:x1] = True mask = ((R - B) > 25) & (R > 140) & box ys, xs = np.nonzero(mask) if len(xs) < 200: raise SystemExit('crescent not found in MOON_BOX -- wrong frame?') lo, hi = {}, {} for x, y in zip(xs, ys): lo[x] = min(lo.get(x, 1 << 30), y) hi[x] = max(hi.get(x, -1), y) cols = sorted(hi) # Trim the cusps, where the arc is cut by the umbra and the fit goes soft. inner = [x for x in cols if cols[0] + 8 <= x <= cols[-1] - 8] P = np.array([[x, hi[x]] for x in inner], float) x, y = P[:, 0], P[:, 1] A = np.c_[x, y, np.ones(len(x))] c = np.linalg.lstsq(A, x ** 2 + y ** 2, rcond=None)[0] cx, cy = c[0] / 2, c[1] / 2 r = math.sqrt(c[2] + cx ** 2 + cy ** 2) # Lit thickness on the centre line -> the uncovered fraction of the diameter. mid = int(round(cx)) near = [x for x in cols if abs(x - mid) <= 6] thick = float(np.mean([hi[x] - lo[x] for x in near])) return cx, cy, r, thick def find_skyline(im, x, sky_y=SKY_PROBE_Y): """First sustained drop below sky brightness, scanning down column x.""" lum = np.asarray(im.convert('L')).astype(float) sky = lum[sky_y:sky_y + 40, x].mean() col = lum[:, x] for y in range(sky_y, len(col) - 8): if all(col[y + k] < sky - 18 for k in range(6)): return y raise SystemExit(f'no skyline found in column {x}') # ---------------------------------------------------------------- ephemeris def load_track(): F = json.load(open(FRAMES))['frames'] def mins(t): h, m = t.split(':') return int(h) * 60 + int(m) F = [f for f in F if '05:30' <= f['t'] <= '07:10'] return (np.array([mins(f['t']) for f in F], float), np.array([f['malt'] for f in F]), np.array([f['sdm'] for f in F]), np.array([f['umag'] for f in F])) def hms(m): s = int(round(m * 60)) return f'{s // 3600:02d}:{s // 60 % 60:02d}:{s % 60:02d}' def build(frame_path, write_figure=True): im = Image.open(frame_path) cx, cy, r, thick = find_crescent(im) d_px = 2 * r ridge_y = find_skyline(im, int(round(cx))) drop_px = ridge_y - cy t, alt, sd, umag = load_track() # The semidiameter is flat to four decimals across the whole window, so # calibration does not depend on which minute the frame is. sd_frame = float(np.mean(sd)) deg_per_px = (2 * sd_frame) / d_px drop_deg = drop_px * deg_per_px def when(a): """Clock time at which the Moon's centre stands at altitude a.""" return float(np.interp(-a, -alt, t)) # altitude is descending print(f'moon edge fit centre ({cx:.1f}, {cy:.1f}) R {r:.2f} px ' f'D {d_px:.1f} px') print(f'scale {deg_per_px*1000:.4f} millidegrees per pixel') print(f'skyline below Moon {drop_px:.0f} px = {drop_deg:.3f}°') print() # --- the part that needs no timestamp ------------------------------------ # The descent rate is not constant, so integrate along the real altitude # curve rather than dividing by a single rate. span = [6 * 60 + 44 + 0.5 * k for k in range(20)] lead = [] for tf in span: a0 = float(np.interp(tf, t, alt)) lead.append((when(a0 - drop_deg + sd_frame) - tf, # lower limb touches when(a0 - drop_deg - sd_frame) - tf)) # disc fully gone gone = [g for _, g in lead] touch = [c for c, _ in lead] minutes_left = float(np.mean(gone)) print('WITHOUT USING THE FRAME TIME AT ALL') print(f' the Moon must fall {drop_deg:.3f}° (centre to skyline) plus its own ' f'radius {sd_frame:.3f}°') print(f' lower limb reaches the skyline {min(touch):.1f}-{max(touch):.1f} min ' f'after this frame') print(f' disc entirely gone {min(gone):.1f}-{max(gone):.1f} min ' f'after this frame') print(' (spread is over every plausible frame time; the answer barely moves)') print() # --- absolute times, from the page's independent dating ------------------ print('ABSOLUTE TIMES, taking the frame time as an input (not measured here)') print(f' {"frame":>8} {"Moon alt":>9} {"ridge":>7} {"touches":>9} {"hidden":>9}') for label in ('06:45', '06:47', '06:49', '06:52'): tf = float(np.interp(0, [0], [0])) if False else ( int(label[:2]) * 60 + int(label[3:])) ma = float(np.interp(tf, t, alt)) ridge = ma - drop_deg print(f' {label:>8} {ma:8.3f}° {ridge:6.3f}° ' f'{hms(when(ridge + sd_frame)):>9} {hms(when(ridge - sd_frame)):>9}') print(' (the page dates this frame 06:47 from its eclipse phase)') print() print(' Sunrise, flat horizon 07:02:23 Selenelion window 07:02:23 - 07:04:22') print(' In every row the Moon is gone before the Sun is up.') print() # --- checks that are printed but not leaned on --------------------------- covered = 1.0 - thick / d_px print(f'loose check (NOT a claim): lit thickness {thick:.0f} px -> umbral ' f'magnitude {covered:.2f}; threshold choice swings this 0.82-0.88.') t_ref = 6 * 60 + 47 ridge_ref = float(np.interp(t_ref, t, alt)) - drop_deg print(f'orientation only (NOT a claim): a {ridge_ref:.2f}° skyline implies a summit of') for km in (20, 30, 40): d = km * 1000.0 dh = d * math.tan(math.radians(ridge_ref)) + d * d / (2 * K_REFRACT * R_EARTH) print(f' {OBSERVER_M + dh:,.0f} m at {km} km', end=';') print('\n range to the ridge is not recoverable from the frame.') t_gone_ref = when(ridge_ref - sd_frame) # --- Earth's shadow, placed rather than fitted --------------------------- # Nothing here is measured off the image except the Moon's own diameter, # which sets the scale. The umbra's size and position come from the # ephemeris at 06:47 and are simply drawn where they fall. Whether the # result lands on the crescent is then a real test, not a fit. fr = json.load(open(FRAMES))['frames'] lo = [x for x in fr if x['t'] == '06:46'][0] hi = [x for x in fr if x['t'] == '06:48'][0] sep, pa = (lo['sep'] + hi['sep']) / 2, (lo['pa'] + hi['pa']) / 2 rum = (lo['rum'] + hi['rum']) / 2 R_sep, R_umb = sep / deg_per_px, rum / deg_per_px ang = math.radians(pa) ux, uy = cx + R_sep * math.sin(ang), cy - R_sep * math.cos(ang) edge_y = uy + R_umb # where its lower edge crosses print(f'\nEarth’s shadow, placed from the ephemeris (nothing fitted to the image):') print(f' umbra radius {R_umb:.1f} px, {R_sep:.1f} px away at PA {pa:+.1f}°') print(f' its lower edge crosses the Moon’s vertical at y = {edge_y:.1f}') print(f' the lit crescent begins at y = {cy - r + (d_px - thick):.1f} ' f'-> agreement {abs(edge_y - (cy - r + (d_px - thick))):.1f} px') if write_figure: draw_figure(im, cx, cy, r, ridge_y, drop_px, drop_deg, d_px, 2 * sd_frame, ridge_ref, minutes_left, t_gone_ref, (ux, uy, R_umb, pa)) return drop_deg, minutes_left # ---------------------------------------------------------------- the figure INK = (255, 255, 255) RED = (240, 84, 70) CYAN = (86, 224, 224) FONTDIR = '/usr/share/fonts/truetype/dejavu' def _font(size, bold=False): name = 'DejaVuSans-Bold.ttf' if bold else 'DejaVuSans.ttf' try: from PIL import ImageFont return ImageFont.truetype(os.path.join(FONTDIR, name), size) except Exception: from PIL import ImageFont return ImageFont.load_default() def _label(d, xy, text, fill, font, anchor='la'): """Text with a dark halo, so it survives whatever is behind it.""" x, y = xy for dx in (-2, -1, 0, 1, 2): for dy in (-2, -1, 0, 1, 2): if dx or dy: d.text((x + dx, y + dy), text, fill=(8, 10, 14), font=font, anchor=anchor) d.text((x, y), text, fill=fill, font=font, anchor=anchor) AMBER = (255, 198, 92) def draw_umbra(d, cx, cy, r, umbra, small, mid): """Earth's shadow where the ephemeris puts it, and its lean from vertical.""" ux, uy, R, pa = umbra d.ellipse([ux - R, uy - R, ux + R, uy + R], outline=AMBER, width=3) d.line([ux - 11, uy, ux + 11, uy], fill=AMBER, width=3) d.line([ux, uy - 11, ux, uy + 11], fill=AMBER, width=3) lx, ly = 1500, 262 d.line([lx - 8, ly + 10, ux + 14, uy + 8], fill=AMBER, width=1) _label(d, (lx, ly), 'centre of Earth’s shadow', AMBER, small) _label(d, (lx, ly + 22), 'placed from the ephemeris,', INK, small) _label(d, (lx, ly + 44), 'not fitted to the picture', INK, small) # the observer's vertical at the Moon, and the lean to that centre d.line([cx, cy, cx, cy - 210], fill=INK, width=2) d.line([cx, cy, ux, uy], fill=AMBER, width=2) _label(d, (150, 196), f'{abs(pa):.0f}° {"left" if pa < 0 else "right"} of straight up', AMBER, mid) _label(d, (150, 224), 'the shadow leans this far from the', INK, small) _label(d, (150, 246), 'observer’s vertical — which is why the', INK, small) _label(d, (150, 268), 'bite is on top, and near the top', INK, small) d.line([408, 208, cx - 14, cy - 96], fill=AMBER, width=1) def draw_figure(im, cx, cy, r, ridge_y, drop_px, drop_deg, d_px, d_deg, ridge_deg, minutes_left, t_gone, umbra=None): im = im.convert('RGB').copy() d = ImageDraw.Draw(im) big, mid, small = _font(30, True), _font(23), _font(20) if umbra: draw_umbra(d, cx, cy, r, umbra, small, mid) # the fitted lunar limb d.ellipse([cx - r, cy - r, cx + r, cy + r], outline=RED, width=3) d.line([cx - 10, cy, cx + 10, cy], fill=RED, width=2) d.line([cx, cy - 10, cx, cy + 10], fill=RED, width=2) _label(d, (cx + r + 14, cy - r + 2), f'the Moon’s own edge, fitted to the lit arc\n{d_px:.0f} px = {d_deg:.3f}°', RED, small) # the skyline x0, x1 = RIDGE_X d.line([x0 - 260, ridge_y, x1 + 240, ridge_y], fill=CYAN, width=3) _label(d, (x0 - 256, ridge_y + 10), f'the skyline toward Los Alamos, {ridge_deg:.2f}° up', CYAN, small) # the drop between them d.line([cx, cy, cx, ridge_y], fill=INK, width=2) for y in (cy, ridge_y): d.line([cx - 15, y, cx + 15, y], fill=INK, width=2) _label(d, (cx + 22, (cy + ridge_y) / 2 - 14), f'{drop_px:.0f} px = {drop_deg:.3f}°', INK, mid) # what the video's own annotation is, so nobody mistakes it for ours _label(d, (300, 648), 'the yellow ring and cursor are the video’s own annotation', (215, 215, 120), small) _label(d, (34, 26), f'{minutes_left:.0f} minutes of Moon left in this shot', INK, big) _label(d, (34, 66), 'It has to fall that gap plus its own radius to clear the skyline. Dating the frame ' f'at 06:47 from its eclipse phase, the disc is gone by {hms(t_gone)}.', INK, mid) _label(d, (34, 96), 'The Sun does not rise until 07:02:23, so the Moon had set before the ' 'selenelion window opened.', INK, mid) im.save(OUT, quality=86, optimize=True) print(f'\nwrote {OUT} ({os.path.getsize(OUT)/1024:.0f} KB)') if __name__ == '__main__': ap = argparse.ArgumentParser() ap.add_argument('--frame', required=True, help='the Santa Fe still to measure') ap.add_argument('--no-figure', action='store_true') a = ap.parse_args() build(a.frame, write_figure=not a.no_figure)