#!/usr/bin/env python3
"""The 8 April 2024 eclipse from Erie, Pennsylvania: where the bite was.

The claim this answers is that the body which crossed the Sun did not travel
along the trajectory measured for the Sun and the Moon, and that this is
impossible. The measurement behind it is real -- a fixed time-lapse camera, one
frame a minute, for weeks -- and the trajectory it establishes is right. What
has gone wrong is which motion the approach is compared against.

  The daily arc is Earth's rotation. EVERY object in the sky shares it: Sun,
  Moon, Jupiter, a star. It is not a property of the Moon and it cannot be
  what carries the Moon onto the Sun.

  The approach is the Moon's orbital motion RELATIVE to the Sun. That belongs
  to the Moon alone, and it is what closes the gap.

They are different motions and there is no reason for them to agree. This
script prints both, for his own site and his own day, so the size of the
expected disagreement can be read rather than argued about -- and then prints
where the model says the bite must sit at each of the moments he photographed,
which is the part his own footage can check.

    python3 eclipse_trajectory.py

Position angles are measured at the Sun, from straight up, positive toward
increasing azimuth -- so 0 is the top of the disc, 90 the side toward the
setting Sun, 180 the bottom. Facing the afternoon Sun in the northern
hemisphere, increasing azimuth is to the observer's right, so 148 deg reads as
"lower right" in a photograph taken with the camera level.
"""

import math

from skyfield.api import load, wgs84

LAT, LON, ELEV = 42.129, -80.085, 180.0        # Erie, Pennsylvania
DATE = (2024, 4, 8)
TZ = -4                                        # EDT

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


def t_local(minute_of_day):
    """Skyfield time from a local clock minute."""
    return ts.utc(*DATE, (minute_of_day // 60) - TZ, minute_of_day % 60)


def altaz(body, t):
    a = SITE.at(t).observe(body).apparent().altaz()
    return a[0].degrees, a[1].degrees


def tangent(alt1, az1, alt0, az0):
    """Offset of (1) from (0) in the sky's tangent plane at (0)."""
    return (az1 - az0) * math.cos(math.radians(alt0)), alt1 - alt0


def great_circle(alt1, az1, alt0, az0):
    """True angular distance on the sky. The tangent-plane offset above is fine
    at the sub-degree separations of the bite work and four degrees short at
    the thirty-seven-degree sweep of an afternoon, so the sweep uses this."""
    def vec(alt, az):
        a, z = math.radians(alt), math.radians(az)
        return (math.cos(a) * math.cos(z), math.cos(a) * math.sin(z), math.sin(a))
    p, q = vec(alt1, az1), vec(alt0, az0)
    return math.degrees(math.acos(max(-1.0, min(1.0, sum(x * y for x, y in zip(p, q))))))


def pa(dx, dy):
    """Position angle from straight up, positive toward increasing azimuth."""
    return math.degrees(math.atan2(dx, dy))


def separation_and_bite(t):
    """How far the Moon is from the Sun, and in which direction it sits.

    The direction is where the bite is: the Moon covers the part of the Sun it
    overlaps, so the notch appears on the side the Moon's centre lies toward.
    """
    sa, sz = altaz(SUN, t)
    ma, mz = altaz(MOON, t)
    dx, dy = tangent(ma, mz, sa, sz)
    return math.hypot(dx, dy), pa(dx, dy)


def directions(t, dt_min=10):
    """The two motions, as position angles.

    (a) where the Sun goes -- the daily arc, the line drawn on his slide
    (b) where the Moon goes RELATIVE to the Sun -- the approach
    """
    t2 = ts.tt_jd(t.tt + dt_min / 1440.0)
    sa0, sz0 = altaz(SUN, t);  sa1, sz1 = altaz(SUN, t2)
    ma0, mz0 = altaz(MOON, t); ma1, mz1 = altaz(MOON, t2)

    arc = pa(*tangent(sa1, sz1, sa0, sz0))

    d0 = tangent(ma0, mz0, sa0, sz0)
    d1 = tangent(ma1, mz1, sa1, sz1)
    approach = pa(d1[0] - d0[0], d1[1] - d0[1])

    gap = abs((approach - arc + 180) % 360 - 180)
    return arc, approach, gap


def greatest_eclipse():
    """Local clock minute of least separation, to the minute."""
    return min(range(13 * 60, 17 * 60), key=lambda m: separation_and_bite(t_local(m))[0])


def clock(m):
    return f'{m // 60:02d}:{m % 60:02d}'


def side(angle):
    named = [(0, 'top'), (45, 'top right'), (90, 'right'), (135, 'lower right'),
             (180, 'bottom'), (-135, 'lower left'), (-90, 'left'), (-45, 'top left')]
    return min(named, key=lambda n: abs((angle - n[0] + 180) % 360 - 180))[1]


def main():
    gm = greatest_eclipse()
    print(f'Erie, Pennsylvania  {LAT}N {abs(LON)}W   8 April 2024')
    print(f'Greatest eclipse at {clock(gm)} EDT '
          f'(his own slide says the Moon must match at 15:18)\n')

    print('THE TWO MOTIONS -- the arc he drew, and the one that closes the gap')
    print(f"  {'EDT':>6} {'Sun alt':>8} {'daily arc':>11} {'approach':>10} {'gap':>7}")
    for m in (gm - 78, gm - 50, gm - 16, gm, gm + 44, gm + 72):
        t = t_local(m)
        arc, approach, gap = directions(t)
        print(f'  {clock(m):>6} {altaz(SUN, t)[0]:7.1f}° {arc:10.0f}° {approach:9.0f}° {gap:6.0f}°')
    print('\n  The two are near enough opposite. Relative to the Sun the Moon creeps')
    print('  eastward while the pair together sweep westward, so the approach runs')
    print('  BACK along the arc rather than forward down it.\n')

    # the same point without any position angles: how far each motion gets
    span = 150                                        # minutes of partial phase
    t0, t1 = t_local(gm - span // 2), t_local(gm + span // 2)
    sa0, sz0 = altaz(SUN, t0); sa1, sz1 = altaz(SUN, t1)
    swept = great_circle(sa1, sz1, sa0, sz0)
    closed = separation_and_bite(t0)[0] + separation_and_bite(t1)[0]
    print(f'ACROSS THE {span}-MINUTE PARTIAL PHASE (both on the sky, neither in his frame)')
    print(f'  the pair sweeps {swept:5.1f}° across the sky')
    print(f'  the Moon moves  {closed:5.1f}° relative to the Sun')
    print(f'  ratio {swept / closed:.0f}:1 -- a line the length of the first cannot')
    print('  describe a motion the length of the second.\n')

    print('WHERE THE BITE IS, at the moments he photographed')
    print(f"  {'EDT':>6} {'vs totality':>12} {'separation':>11} {'direction':>10}   in the frame")
    for d in (-55, -35, -17, -6, 0, 3, 15, 35, 55):
        t = t_local(gm + d)
        sep, bite = separation_and_bite(t)
        label = 'totality' if d == 0 else f'{d:+d} min'
        print(f'  {clock(gm + d):>6} {label:>12} {sep * 60:10.1f}′ {bite:9.0f}°   {side(bite)}')
    print('\n  His 14:28 photograph has the bite at the lower right, which is what')
    print('  the line above it says. The direction then rotates through about 160')
    print('  degrees across the event and finishes at the top -- and he took four')
    print('  frames after totality, so his own set can check that.')


if __name__ == '__main__':
    main()