#!/usr/bin/env python3
"""Measure the Moon in the photograph on p.169, off a render of the book page.

The question is whether the dark part is the Moon's own night side (a phase) or
Earth's umbra (a deep partial eclipse). The tempting test -- fit a circle to the
outer edge of the bright sliver -- cannot answer it, because that edge is the
Moon's own limb in both cases. Two tests can:

  1. Horn span. A terminator bounds a lit hemisphere, so it meets the limb at
     two points exactly 180 degrees apart, whatever the phase. An umbra edge is
     a chord and its horns span less.
  2. Inner-edge curvature. A terminator is a half-ellipse centred on the Moon's
     own centre. The umbra edge is an arc of a much larger circle centred well
     outside the disc.

Usage:
    pdftoppm -f <pdf page> -l <pdf page> -r 300 -png book.pdf page
    python3 measure_p169_moon.py page-169.png
"""

import sys
import numpy as np
from PIL import Image


def fit_circle(pts):
    """Algebraic circle fit. Returns (cx, cy, r, rms residual)."""
    pts = np.asarray(pts, float)
    A = np.c_[2 * pts[:, 0], 2 * pts[:, 1], np.ones(len(pts))]
    b = (pts ** 2).sum(axis=1)
    sol, *_ = np.linalg.lstsq(A, b, rcond=None)
    cx, cy = sol[0], sol[1]
    r = np.sqrt(sol[2] + cx * cx + cy * cy)
    return cx, cy, r, (np.hypot(pts[:, 0] - cx, pts[:, 1] - cy) - r).std()


def measure(path, box=(0.30, 0.60, 0.0, 0.12)):
    a = np.asarray(Image.open(path).convert('RGB')).astype(float)
    h, w, _ = a.shape
    sub = a[int(box[2] * h):int(box[3] * h), int(box[0] * w):int(box[1] * w)]
    r, g, b = sub[..., 0], sub[..., 1], sub[..., 2]
    lum = 0.299 * r + 0.587 * g + 0.114 * b
    mask = (lum > lum.max() * 0.72) & ((r + g) / 2 - b > 25)
    ys, xs = np.nonzero(mask)
    if len(xs) < 200:
        sys.exit('found no lit region -- adjust the crop box')

    outer = [(x, ys[xs == x].max()) for x in range(xs.min() + 8, xs.max() - 7) if (xs == x).any()]
    inner = [(x, ys[xs == x].min()) for x in range(xs.min() + 10, xs.max() - 9) if (xs == x).any()]
    cx, cy, R, e1 = fit_circle(outer)
    ux, uy, RU, e2 = fit_circle(inner)

    print(f"Moon's limb   centre ({cx:7.1f},{cy:7.1f})  r {R:7.2f} px   rms {e1:.2f}")
    print(f"inner edge    centre ({ux:7.1f},{uy:7.1f})  r {RU:7.2f} px   rms {e2:.2f}")
    print(f"\n  inner-edge radius        {RU / R:6.3f} x the Moon's radius")
    print(f"  its centre lies          {np.hypot(ux - cx, uy - cy) / R:6.3f} Moon-radii from the Moon's centre")
    print(f"  at a position angle of   {np.degrees(np.arctan2(ux - cx, -(uy - cy))) % 360:6.1f} deg from straight up")

    onlimb = np.hypot(xs - cx, ys - cy) > R * 0.90
    th = np.degrees(np.arctan2(xs[onlimb] - cx, -(ys[onlimb] - cy))) % 360
    span = th.max() - th.min()
    print(f"\n  the lit arc runs         {span:6.1f} deg along the limb")
    print("                           a phase must run exactly 180 -- its horns are a diameter")

    col = ys[np.abs(xs - cx) < 3]
    lit = (col.max() - col.min()) / (2 * R)
    print(f"  lit across the middle    {lit * 100:6.1f}% of the diameter -> {100 - lit * 100:.0f}% covered")


if __name__ == '__main__':
    measure(sys.argv[1] if len(sys.argv) > 1 else 'page-169.png')