#!/usr/bin/env python3
"""Find EVERY selenelion visible from Santa Fe, without a hand-written list.

scripts/selenelion_dates.py measures the duration of each event in a list that
was assembled by hand. An adversarial review asked the obvious question about a
hand-assembled list -- is anything missing? -- and named 15 April 1995 as a
candidate. This script answers it mechanically: sweep every full Moon in the
window, keep the ones where Earth's umbra touches the disc at all, and for each
of those ask whether Santa Fe ever had the Sun's upper limb and the eclipsed
Moon's upper limb above opposite horizons at the same instant.

Nothing is filtered by hand. What comes out is the survey.

Conventions match scripts/selenelion_2011.py and selenelion_dates.py exactly:
upper limbs, each refracted at its own altitude (Bennett 1982, clamped where the
formula turns), apparent positions, sea-level horizon, Danjon-enlarged umbra,
and "eclipsed" means umbral magnitude > 0 -- the umbra touching the disc, not
the penumbra, which is invisible to the eye and to a camera.

    python3 selenelion_survey.py [start_year] [end_year]
"""

import math
import sys

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

R_EARTH, R_SUN, R_MOON = 6378.137, 696000.0, 1737.4
DANJON = 1.02
LAT, LON, ELEV = 35.6870, -105.9378, 2130.0

ts = load.timescale()
eph = load('de421.bsp')
earth, sun, moon = eph['earth'], eph['sun'], eph['moon']
topos = wgs84.latlon(LAT, LON, elevation_m=ELEV)
site = earth + topos


def bennett(h):
    return (1.0 / math.tan(math.radians(h + 7.31 / (h + 4.4)))) / 60.0


R_MAX = bennett(-1.2)
H_CLAMP = -1.2 - R_MAX


def refract(h_true):
    if h_true <= H_CLAMP:
        return h_true + R_MAX
    h = h_true
    for _ in range(60):
        h = h_true + bennett(h)
    return h


def umbral_magnitude(t):
    """Geocentric umbral magnitude: 0 = umbra tangent to the disc, 1 = totality."""
    e = earth.at(t)
    s = e.observe(sun).apparent().position.km
    m = e.observe(moon).apparent().position.km
    d_sun, d_moon = np.linalg.norm(s), np.linalg.norm(m)
    axis = -s / d_sun
    sep_km = float(np.linalg.norm(m - float(np.dot(m, axis)) * axis))
    apex = d_sun * R_EARTH / (R_SUN - R_EARTH)
    r_umbra_km = R_EARTH * (1.0 - d_moon / apex) * DANJON
    return (r_umbra_km + R_MOON - sep_km) / (2 * R_MOON)


def local_state(t):
    """Both upper limbs as seen from Santa Fe, plus the local umbral magnitude."""
    a = site.at(t)
    s, m = a.observe(sun).apparent(), a.observe(moon).apparent()
    s_alt, m_alt = s.altaz()[0].degrees, m.altaz()[0].degrees
    sd_s = math.degrees(math.asin(R_SUN / s.distance().km))
    sd_m = math.degrees(math.asin(R_MOON / m.distance().km))
    return refract(s_alt + sd_s), refract(m_alt + sd_m)


def umbral_eclipses(y0, y1):
    """Every full Moon in the window whose umbral magnitude peaks above zero."""
    t0, t1 = ts.utc(y0, 1, 1), ts.utc(y1 + 1, 1, 1)
    times, phases = almanac.find_discrete(t0, t1, almanac.moon_phases(eph))
    out = []
    for t, ph in zip(times, phases):
        if ph != 2:                      # 2 = full moon
            continue
        # peak umbral magnitude within +/- 4 h of syzygy, on a 4-minute grid
        offs = np.arange(-240, 241, 4) * 60.0
        best = max((umbral_magnitude(ts.tt_jd(t.tt + o / 86400.0)), o) for o in offs)
        if best[0] > 0:
            out.append((ts.tt_jd(t.tt + best[1] / 86400.0), best[0]))
    return out


def santa_fe_window(t_greatest):
    """Longest run with both upper limbs up AND the Moon in the umbra.

    Scanned +/- 4 h around greatest eclipse at 10 s, then each edge bisected.
    Returns (start, end, seconds, umbral magnitude at the midpoint) or None.
    """
    step = 10.0
    offs = np.arange(-4 * 3600, 4 * 3600 + step, step)

    def ok(off):
        t = ts.tt_jd(t_greatest.tt + off / 86400.0)
        s_limb, m_limb = local_state(t)
        return s_limb >= 0 and m_limb >= 0 and umbral_magnitude(t) > 0

    flags = [ok(o) for o in offs]
    if not any(flags):
        return None
    i0 = flags.index(True)
    i1 = len(flags) - 1 - flags[::-1].index(True)

    def edge(lo, hi, want_at_hi):
        for _ in range(30):
            mid = (lo + hi) / 2
            if ok(mid) == want_at_hi:
                hi = mid
            else:
                lo = mid
        return (lo + hi) / 2

    a = edge(offs[i0 - 1], offs[i0], True) if i0 > 0 else offs[i0]
    b = edge(offs[i1 + 1], offs[i1], True) if i1 < len(offs) - 1 else offs[i1]
    mid = ts.tt_jd(t_greatest.tt + ((a + b) / 2) / 86400.0)
    return (ts.tt_jd(t_greatest.tt + a / 86400.0),
            ts.tt_jd(t_greatest.tt + b / 86400.0),
            b - a, umbral_magnitude(mid))


def mst(t):
    d = t.utc_datetime()
    x = (d.hour * 3600 + d.minute * 60 + d.second - 7 * 3600) % 86400
    return f'{x // 3600:02d}:{x % 3600 // 60:02d}:{x % 60:02d}'


def side(t):
    """Which way the umbra's centre lies from the Moon's, in the observer's frame."""
    a = site.at(t)
    s, m = a.observe(sun).apparent(), a.observe(moon).apparent()
    e = earth.at(t)
    s_geo = e.observe(sun).apparent().position.km
    m_geo = e.observe(moon).apparent().position.km
    axis = -s_geo / np.linalg.norm(s_geo)
    u_pt = axis * np.linalg.norm(m_geo)
    from skyfield.positionlib import Geocentric
    obs = topos.at(t).position.km
    u_alt, u_az, _ = Geocentric((u_pt - obs) / 149597870.700, t=t).frame_latlon(topos)
    m_alt, m_az = m.altaz()[0].degrees, m.altaz()[1].degrees
    dalt = u_alt.degrees - m_alt
    daz = ((u_az.degrees - m_az + 540) % 360) - 180
    dx = daz * math.cos(math.radians(m_alt))
    pa = math.degrees(math.atan2(dx, dalt))
    # pa is measured from straight up in the observer's frame, positive toward
    # increasing azimuth. Eight sectors of 45 degrees, named as a reader would.
    names = ['the top', 'upper right', 'the right', 'lower right',
             'the bottom', 'lower left', 'the left', 'upper left']
    return pa, names[int(((pa % 360) + 22.5) // 45) % 8]


def main():
    y0 = int(sys.argv[1]) if len(sys.argv) > 1 else 1995
    y1 = int(sys.argv[2]) if len(sys.argv) > 2 else 2030
    print(f'Every selenelion visible from Santa Fe, {y0}-{y1}')
    print('swept mechanically: no hand-written event list\n')
    ecl = umbral_eclipses(y0, y1)
    print(f'{len(ecl)} umbral eclipses in the window\n')
    print(f'{"date":>13} {"MST window":>21} {"duration":>12} {"umag":>6} '
          f'{"tilt":>7}  where the shadow sat')
    total = 0.0
    hits = 0
    for t_g, _ in ecl:
        w = santa_fe_window(t_g)
        if w is None:
            continue
        t0, t1, sec, umag = w
        mid = ts.tt_jd((t0.tt + t1.tt) / 2)
        pa, name = side(mid)
        d = t0.utc_datetime()
        # morning (Moon setting) or evening (Moon rising)?
        a0 = site.at(t0).observe(moon).apparent().altaz()[0].degrees
        a1 = site.at(t1).observe(moon).apparent().altaz()[0].degrees
        kind = 'morning' if a1 < a0 else 'evening'
        hits += 1
        total += sec
        print(f'{d.strftime("%d %b %Y"):>13} {mst(t0)}-{mst(t1)} '
              f'{sec:8.0f} s {umag:6.2f} {pa:+6.1f}°  {name} ({kind})')
    print(f'\n{hits} events, {total / 60:.0f} min {total % 60:02.0f} s in total')


if __name__ == '__main__':
    main()