#!/usr/bin/env python3
"""
The 10 December 2011 lunar eclipse, seen from Santa Fe, New Mexico.

Globe Deconstruction? pp. 167-168 carries a photograph captioned "taken at
sunrise in the Santa Fe area, looking toward Los Alamos", with the claim that
"the heliocentric model dictates that the ball's shadow would be on the bottom
of the Moon. Here we see the opposite."

This script computes, from JPL DE421, where Earth's umbra actually sat on the
Moon's disc in that observer's own frame, and how much atmospheric refraction a
simultaneous view of the rising Sun and the setting eclipsed Moon required.

Method note. The umbra is a real cone in space, so we do not fudge it with a
position angle taken from a star chart. We build the shadow axis as the
anti-solar ray from Earth's centre, take the point on it at the Moon's
distance, and then observe BOTH that point and the Moon from the observer's
position, in the observer's horizontal frame. The offset that comes out is what
a camera on that spot records. Nothing depends on a parallactic-angle sign
convention.

Cross-checks printed at the end: greatest-eclipse time, gamma and umbral
magnitude against EclipseWise; peak visible magnitude against timeanddate.

    pip install skyfield numpy
    python3 selenelion_2011.py
"""

import numpy as np
from skyfield.api import load, wgs84
from skyfield.positionlib import Geocentric
from skyfield import almanac  # noqa: F401  (kept for readers extending this)

AU_KM = 149597870.700
R_EARTH, R_SUN, R_MOON = 6378.137, 696000.0, 1737.4

# Santa Fe, NM. Los Alamos is ~35 km WNW; at these altitudes the difference
# moves nothing here by more than a few hundredths of a degree.
LAT, LON, ELEV = 35.6870, -105.9378, 2130.0
MST = -7 * 60  # minutes from UT

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


# ---------------------------------------------------------------- refraction
def bennett(h_app_deg):
    """Bennett (1982): refraction in degrees for a body at apparent altitude h."""
    return (1.0 / np.tan(np.radians(h_app_deg + 7.31 / (h_app_deg + 4.4)))) / 60.0


# The deepest bend we are willing to model: the ray that arrives at an
# apparent -1.2 deg, which covers any horizon dip a ground observer can have.
# Below the true altitude that maps to it, hold the bending there -- further
# down, Bennett walks toward its pole at -4.4 deg and returns nonsense.
R_MAX = bennett(-1.2)                 # 0.8805 deg = 52.8'
H_CLAMP = -1.2 - R_MAX                # -2.0805 deg true altitude


def refract(h_true_deg):
    """True altitude -> apparent altitude, by iterating Bennett.

    Bennett's formula is written in terms of APPARENT altitude and is only
    meaningful at or above the horizon, staying usable a little way below it for
    an observer whose horizon is depressed. Fed a true altitude well below, the
    fixed-point iteration walks past the formula's pole at h = -4.4 and settles
    on a spurious root: it returns +55' of bending for a body 2.4 deg down, and
    NEGATIVE bending below about -4.9 deg. Neither is physical. So invert down
    to an apparent -1.2 deg and hold the bending constant past that, as the
    other scripts here do.

    Below the horizon the number is a modelling convenience with no observable
    meaning -- there is no ray from the body to the eye -- and it is never
    published as an apparent altitude. Where the page prints an altitude for a
    body below the horizon, it prints the geometric one.
    """
    if h_true_deg <= H_CLAMP:
        return h_true_deg + R_MAX
    h = h_true_deg
    for _ in range(60):
        h = h_true_deg + bennett(h)
    return h


# ------------------------------------------------------------------ geometry
def topocentric_altaz(vec_km, t):
    """alt/az of a geocentric POINT, seen from the observer (parallax included)."""
    d = np.asarray(vec_km) - topos.at(t).position.km
    alt, az, dist = Geocentric(d / AU_KM, t=t).frame_latlon(topos)
    return alt.degrees, az.degrees, dist.km


def direction_altaz(vec_km, t):
    """alt/az of a geocentric DIRECTION, in the observer's horizontal frame."""
    v = np.asarray(vec_km, dtype=float)
    v = v / np.linalg.norm(v) * 1.0e6
    alt, az, _ = Geocentric(v / AU_KM, t=t).frame_latlon(topos)
    return alt.degrees, az.degrees


def state(t):
    e = earth.at(t)
    m = e.observe(moon).apparent().position.km
    s = e.observe(sun).apparent().position.km
    r_moon, d_sun = np.linalg.norm(m), np.linalg.norm(s)

    # the umbra axis, at the Moon's distance
    axis_pt = (-s / d_sun) * r_moon

    m_alt, m_az, m_rng = topocentric_altaz(m, t)
    u_alt, u_az, _ = topocentric_altaz(axis_pt, t)
    s_alt, s_az, _ = topocentric_altaz(s, t)

    # offset of the umbra's centre from the Moon's centre, in the observer's
    # frame: PA measured from straight up, positive toward increasing azimuth
    # (that is, to the right of vertical when you are facing the Moon).
    d_up = u_alt - m_alt
    d_rt = (u_az - m_az) * np.cos(np.radians(m_alt))
    sep = np.hypot(d_up, d_rt)
    pa = np.degrees(np.arctan2(d_rt, d_up))

    # angular radii
    pi_moon = np.degrees(np.arcsin(R_EARTH / r_moon))
    pi_sun = np.degrees(np.arcsin(R_EARTH / d_sun))
    sd_sun = np.degrees(np.arcsin(R_SUN / d_sun))
    r_umbra = 1.02 * (pi_moon + pi_sun - sd_sun)      # Danjon 1/50 enlargement
    sd_moon = np.degrees(np.arcsin(R_MOON / m_rng))
    umag = (r_umbra + sd_moon - sep) / (2.0 * sd_moon)

    mg_alt, mg_az = direction_altaz(m, t)
    sg_alt, sg_az = direction_altaz(s, t)
    ug_alt, ug_az = direction_altaz(axis_pt, t)

    return dict(m_alt=m_alt, m_app=refract(m_alt), m_az=m_az,
                s_alt=s_alt, s_app=refract(s_alt), s_az=s_az,
                u_alt=u_alt, u_app=refract(u_alt), u_az=u_az,
                sep=sep, pa=pa, r_umbra=r_umbra, sd_moon=sd_moon, sd_sun=sd_sun, umag=umag,
                pi_moon=pi_moon,
                mg_alt=mg_alt, mg_az=mg_az, sg_alt=sg_alt, sg_az=sg_az,
                ug_alt=ug_alt, ug_az=ug_az)


def at_mst(hh, mm):
    return ts.utc(2011, 12, 10, 0, 0, ((hh * 60 + mm) - MST) * 60)


def hhmm(t):
    u = t.utc
    total = (u.hour * 60 + u.minute + MST) % 1440
    return f"{total // 60:02d}:{total % 60:02d}"


# ------------------------------------------------------- the survey on the page
def morning_selenelions(y0=1995, y1=2030):
    """Every morning selenelion from this site: an umbral eclipse with the Moon
    setting in the west and the Sun already risen in the east. Prints where the
    umbra sat on the disc in each case. Takes a minute or two."""
    obs = earth + topos

    # 1. coarse pass: when is the Moon within 1.6 deg of the antisolar point
    t = ts.utc(y0, 1, 1, np.arange(0, int((y1 - y0) * 365.25 * 8)) * 3)
    e = earth.at(t)
    m = e.observe(moon).apparent().position.km
    s_ = e.observe(sun).apparent().position.km
    sep = np.degrees(np.arccos(np.clip(
        ((m / np.linalg.norm(m, axis=0)) * (-s_ / np.linalg.norm(s_, axis=0))).sum(axis=0), -1, 1)))
    idx = np.where(sep < 1.6)[0]
    groups, cur = [], [idx[0]]
    for i in idx[1:]:
        if i - cur[-1] <= 2:
            cur.append(i)
        else:
            groups.append(cur)
            cur = [i]
    groups.append(cur)
    base = ts.utc(y0, 1, 1).tt

    def clock(pa):
        a = (pa + 360) % 360
        for lim, name in [(25, 'the top'), (70, 'upper right'), (110, 'the right side'),
                          (155, 'lower right'), (205, 'the bottom'), (250, 'lower left'),
                          (290, 'the left side'), (335, 'upper left')]:
            if a < lim:
                return name
        return 'the top'

    print(f"\nMorning selenelions from this site, {y0}-{y1}")
    print("  (umbral eclipse, Moon setting in the west, Sun already up in the east)")
    for g in groups:
        # 2. vectorised altaz over the group, to find candidate minutes cheaply
        secs = np.arange((g[0] - 1) * 10800, (g[-1] + 2) * 10800, 60.0)
        tt = ts.tt_jd(base + secs / 86400.0)
        ma, mz, _ = obs.at(tt).observe(moon).apparent().altaz()
        sa, sz, _ = obs.at(tt).observe(sun).apparent().altaz()
        ok = ((refract(ma.degrees) + 0.259 > 0) & (refract(sa.degrees) + 0.267 > 0)
              & (mz.degrees > 200) & (mz.degrees < 340))
        if not ok.any():
            continue
        # 3. the expensive call, only on the surviving minutes
        rows = [state(ts.tt_jd(base + x / 86400.0)) for x in secs[ok]]
        rows = [r for r in rows if r['umag'] > 0.0]
        if not rows:
            continue
        mid = rows[len(rows) // 2]
        u = ts.tt_jd(base + secs[ok][0] / 86400.0).utc
        print(f"  {u.year}-{u.month:02d}-{u.day:02d}  Moon az {mid['m_az']:5.1f}"
              f"  peak umbral magnitude {max(r['umag'] for r in rows):4.2f}"
              f"  umbra at {mid['pa']:+6.1f} deg from up -> {clock(mid['pa'])}")


# --------------------------------------------------------------------- run
if __name__ == '__main__':
    print("10 December 2011, from Santa Fe NM  %.4f N  %.4f W  %.0f m\n"
          % (LAT, -LON, ELEV))

    print("Where the umbra sat, in the observer's own frame")
    print(f"{'MST':>6} {'Moon alt':>9} {'az':>7} {'Sun alt':>9} {'az':>7}"
          f" {'offset':>7} {'PA from up':>11} {'covered':>8}")
    for hh, mm in [(4, 40), (5, 15), (5, 46), (6, 15), (6, 40), (6, 54), (7, 2)]:
        r = state(at_mst(hh, mm))
        side = 'left' if r['pa'] < 0 else 'right'
        print(f"{hh:02d}:{mm:02d}  {r['m_app']:8.2f} {r['m_az']:7.2f}"
              f" {r['s_app']:9.2f} {r['s_az']:7.2f} {r['sep']:7.3f}"
              f" {abs(r['pa']):6.1f} {side:<4} {r['umag'] * 100:7.0f}%")

    print("\nThe one-line rule, checked: the shadow's centre is the antisolar point")
    r = state(at_mst(7, 2))
    print(f"  Sun direction        alt {r['sg_alt']:+7.3f}  az {r['sg_az']:7.2f}")
    print(f"  antisolar (implied)  alt {-r['sg_alt']:+7.3f}  az {(r['sg_az'] + 180) % 360:7.2f}")
    print(f"  umbra axis (computed) alt {r['ug_alt']:+7.3f}  az {r['ug_az']:7.2f}")
    print(f"  Moon                 alt {r['mg_alt']:+7.3f}  az {r['mg_az']:7.2f}")
    print(f"  -> the Moon lay {r['ug_alt'] - r['mg_alt']:.3f} deg BELOW the shadow's centre")

    print("\nThe selenelion window, level horizon, standard refraction")
    print("  criterion: apparent altitude of the upper limb above zero, for both at once")

    def upper_limbs(t):
        # Each limb refracted AT ITS OWN altitude. Refracting the centre and then
        # adding the semidiameter treats the disc as rigid and puts the upper limb
        # about 2.2' too high near the horizon, because refraction is weaker up
        # there than at the centre -- the same rigid-disc error that used to be in
        # the shadow-edge figure. It widened this window by roughly fifty seconds.
        # Semidiameters are computed from the actual ranges, not assumed.
        r = state(t)
        return (refract(r['m_alt'] + r['sd_moon']),
                refract(r['s_alt'] + r['sd_sun']))

    def bisect(sign_fn, lo=14 * 3600.0, hi=14 * 3600.0 + 480.0):
        for _ in range(60):
            mid = 0.5 * (lo + hi)
            if sign_fn(ts.utc(2011, 12, 10, 0, 0, mid)) > 0:
                lo = mid
            else:
                hi = mid
        return 0.5 * (lo + hi)

    def as_mst(secs):
        total = secs / 60.0 - 420.0
        return f"{int(total // 60):02d}:{int(total % 60):02d}:{secs % 60:04.1f}"

    t_set = bisect(lambda t: upper_limbs(t)[0])
    t_rise = bisect(lambda t: -upper_limbs(t)[1])
    print(f"  Sun's upper limb clears the horizon  {as_mst(t_rise)} MST")
    print(f"  Moon's upper limb drops below it     {as_mst(t_set)} MST")
    print(f"  both visible at once for {t_set - t_rise:.0f} seconds")
    for s in (t_rise, t_set):
        r = state(ts.utc(2011, 12, 10, 0, 0, s))
        print(f"    {as_mst(s)}  umbra offset {r['sep']:.3f} deg,"
              f" {abs(r['pa']):.1f} deg {'left' if r['pa'] < 0 else 'right'} of up,"
              f" {r['umag'] * 100:.0f}% covered")

    print("\nHow much refraction a simultaneous view needed")
    prev = None
    for s in range(13 * 3600 + 30 * 60, 14 * 3600 + 30 * 60, 10):
        t = ts.utc(2011, 12, 10, 0, 0, s)
        r = state(t)
        d = r['m_alt'] - r['s_alt']
        if prev is not None and prev > 0 >= d:
            print(f"  their true altitudes cross at {hhmm(t)} MST, both at {r['m_alt']:+.3f} deg")
            print(f"  centres need   {-r['m_alt'] * 60:5.1f}' of lift each")
            # at the crossing the two true altitudes are equal by definition, so
            # the limb each body has to raise differs only by its own semidiameter
            h = -r['m_alt']
            print(f"  upper limbs need {(h - r['sd_moon']) * 60:5.1f}' (Moon) and "
                  f"{(h - r['sd_sun']) * 60:5.1f}' (Sun) -- semidiameters computed, "
                  f"not assumed")
            print(f"  standard refraction at the horizon supplies {bennett(0.0) * 60:5.1f}'")
            break
        prev = d

    print("\nThe same refraction, measurable in the frame")
    print(f"  Sun on the horizon: R(0.00 deg) = {bennett(0.0) * 60:.2f}',"
          f" R(0.53 deg) = {bennett(0.533) * 60:.2f}'")
    print(f"    -> a 32.0' Sun renders {32.0 - (bennett(0.0) - bennett(0.533)) * 60:.1f}' tall,"
          f" 32.0' wide")
    for h in (0.5, 1.0, 2.0, 5.0):
        sq = (bennett(h - 0.25) - bennett(h + 0.25)) * 60
        print(f"  Moon at {h:4.1f} deg: 31.0' disc squashed by {sq:4.2f}' ({sq / 31.0 * 100:4.1f}%)")

    print("\nCross-checks")
    best = None
    for s in range(13 * 3600, 16 * 3600, 30):
        t = ts.utc(2011, 12, 10, 0, 0, s)
        r = state(t)
        if best is None or r['sep'] < best[1]:
            best = (t, r['sep'], r)
    t, sep, r = best
    u = t.utc
    print(f"  greatest eclipse  {u.hour:02d}:{u.minute:02d} UT"
          f"   gamma {sep / r['pi_moon']:.4f}   U.Mag {r['umag']:.4f}")
    print(f"  EclipseWise       14:32 UT (14:33 TD)   gamma -0.3882   U.Mag 1.1061")
    r54 = state(at_mst(6, 54))
    print(f"  our magnitude at 06:54 MST {r54['umag']:.3f}"
          f"   timeanddate peak for this location 0.890 at 06:54")

    morning_selenelions()


# ------------------------------------------- data for the docs/moon-tilt widget
def export_widget_data(path=None, start_mst=(4, 0), end_mst=(7, 10), step_min=2):
    """Write the frame table the two-panel animation on the page is driven from.

    Everything the widget draws comes from here, so the picture cannot drift
    away from the computation. Vectors are ICRF, in Earth radii.
    """
    import json, os
    if path is None:
        path = os.path.join(os.path.dirname(os.path.dirname(os.path.abspath(__file__))),
                            'docs', 'selenelion', 'selenelion-frames.json')
    frames = []
    t0 = start_mst[0] * 60 + start_mst[1]
    t1 = end_mst[0] * 60 + end_mst[1]
    for mins in range(t0, t1 + 1, step_min):
        t = ts.utc(2011, 12, 10, 0, 0, (mins - MST) * 60)
        e = earth.at(t)
        m = e.observe(moon).apparent().position.km
        s = e.observe(sun).apparent().position.km
        o = topos.at(t).position.km
        r = state(t)
        frames.append(dict(
            t=f"{mins // 60:02d}:{mins % 60:02d}",
            moon=[round(v / R_EARTH, 4) for v in m],
            sun=[round(float(v), 6) for v in (s / np.linalg.norm(s))],
            obs=[round(v / R_EARTH, 4) for v in o],
            malt=round(r['m_app'], 3), maz=round(r['m_az'], 2),
            salt=round(r['s_app'], 3), saz=round(r['s_az'], 2),
            # the same altitudes with refraction switched off, so the widget
            # can draw the sky the geometry alone would have delivered
            maltg=round(r['m_alt'], 3), saltg=round(r['s_alt'], 3),
            sdm=round(r['sd_moon'], 4), sds=round(r['sd_sun'], 4),
            # the umbra's own axis: where Earth's shadow centre sits in the sky,
            # refracted and not, so the shadow can be drawn against the Moon's track
            ualt=round(r['u_app'], 3), ualtg=round(r['u_alt'], 3),
            uaz=round(r['u_az'], 2),
            # The shadow's lower edge, refracted AT ITS OWN ALTITUDE. Refracting
            # the axis and then subtracting the radius treats the umbra as a rigid
            # disc and understates the edge's lift by up to 6' near the horizon,
            # because the edge hangs below the axis and so sits in thicker air.
            uedgeg=round(r['u_alt'] - r['r_umbra'], 4),
            uedge=round(refract(r['u_alt'] - r['r_umbra']), 4),
            sep=round(r['sep'], 4), pa=round(r['pa'], 2),
            umag=round(r['umag'], 4),
            rum=round(r['r_umbra'], 4), rmo=round(r['sd_moon'], 4)))
    meta = dict(
        site=dict(name='Santa Fe, NM', lat=LAT, lon=LON, elev_m=ELEV),
        date='2011-12-10', tz='MST (UT-7)',
        note='ICRF vectors in Earth radii; sep/pa are the umbra centre offset from '
             'the Moon centre in the observer horizontal frame, pa measured from '
             'straight up, positive toward increasing azimuth. Angles in degrees.',
        source='scripts/selenelion_2011.py, JPL DE421',
        frames=frames)
    blob = json.dumps(meta, separators=(',', ':'))
    with open(path, 'w') as f:
        f.write(blob)
    print(f"\nwrote {len(frames)} frames to {path}")

    # the page carries its own copy so it works as a single file; keep them in step
    page = os.path.join(os.path.dirname(path), 'index.html')
    if os.path.exists(page):
        import re
        html = open(page, encoding='utf-8').read()
        tag = '<script id="selenelion-frames" type="application/json">'
        new = re.sub(re.escape(tag) + r'.*?</script>',
                     lambda m: tag + blob + '</script>', html, count=1, flags=re.S)
        if new != html:
            open(page, 'w', encoding='utf-8').write(new)
            print('refreshed the inline copy in', page)
    return path


# --------------------------------------- the two stations in the p.169 source video
def two_stations():
    """The video the book's QR link reaches shows two observers on the same
    morning: the New Mexico footage the book crops p.169 from, and a second
    clip from Monks Mound at Cahokia, just east of St Louis. Same Moon, same shadow, same minute -- and
    the umbra is tilted about thirty degrees differently, because the two
    observers' verticals point in different directions. Nobody set that test
    up; it is simply in the video."""
    global topos
    keep = topos
    print("\nTwo stations, 10 December 2011")
    for name, lat, lon, elev, tz in [('Santa Fe NM', 35.687, -105.938, 2130, -7 * 60),
                                     ('Monks Mound, Cahokia IL', 38.6606, -90.0619, 160, -6 * 60)]:
        topos = wgs84.latlon(lat, lon, elevation_m=elev)
        print(f"\n  {name}")
        print(f"    {'local':>7} {'moon alt':>9} {'sun alt':>8} {'covered':>9} {'umbra tilt':>12}")
        for mins in range(6 * 60 + 20, 7 * 60 + 13, 6):
            t = ts.utc(2011, 12, 10, 0, 0, (mins - tz) * 60)
            r = state(t)
            if refract(r['m_alt'] + r['sd_moon']) < 0:
                break
            side = 'left of up' if r['pa'] < 0 else 'right of up'
            cov = max(r['umag'], 0.0) * 100
            print(f"    {mins // 60:02d}:{mins % 60:02d}   {refract(r['m_alt']):9.2f}"
                  f" {refract(r['s_alt']):8.2f} {cov:8.1f}% {abs(r['pa']):6.1f} {side}")
    topos = keep