#!/usr/bin/env python3
"""The moon tilt illusion, in closed form, applied to Example 1 on p.168.

That example is the panorama from Earth Science Stack Exchange question 14809,
"What is this sun and moon photographic anomaly?", captioned in the book as it
is on the site: "Photo taken 16 May 2016". The questioner's own diagram carries
the estimate that "the sun would have to be about 45deg above and slightly in
front of the moon", offered as an impossibility.

The angles are published. Myers-Beaghton, A. K., & Myers, A. L. (2014),
"The Moon Tilt Illusion", KoG 18, 53-59, gives:

    Eq. (3)   tan a = (cos nm tan ns - sin nm cos dphi) / sin dphi
    Eq. (12)  tan b = -(sin nm - sin ns) / (cos ns sin dphi)
    Eq. (13)  delta = a - b

where nm and ns are the altitudes of moon and sun, dphi is the absolute
azimuth difference, a is the observed direction of the incoming light measured
from the horizontal, and b is the direction a viewer naively expects. Positive
is above the horizontal.

This script checks that implementation against the paper's own worked example
before using it, then runs it over the evening of 16 May 2016 in Scotland.

    pip install skyfield numpy
    python3 moon_tilt_scotland.py
"""

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

D, G = np.radians, np.degrees


def alpha(eta_m, eta_s, dphi):
    """Observed slope of the incoming light, degrees above the horizontal."""
    m, s, d = D(eta_m), D(eta_s), D(abs(dphi))
    return G(np.arctan2(np.cos(m) * np.tan(s) - np.sin(m) * np.cos(d), np.sin(d)))


def beta(eta_m, eta_s, dphi):
    """Slope the observer naively expects, degrees above the horizontal."""
    m, s, d = D(eta_m), D(eta_s), D(abs(dphi))
    return G(np.arctan2(-(np.sin(m) - np.sin(s)), np.cos(s) * np.sin(d)))


def check():
    """The paper's own Figure 1: eta_m 45, eta_s -15, dphi 128."""
    a, b = alpha(45, -15, 128), beta(45, -15, 128)
    print("Check against the paper's worked example (their Figure 1)")
    print(f"  alpha {a:6.1f}   paper says 17")
    print(f"  beta  {b:6.1f}   paper says -52")
    print(f"  delta {a - b:6.1f}   paper says 69")
    ok = abs(a - 17) < 0.6 and abs(b + 52) < 0.6
    print("  ->", "implementation agrees" if ok else "MISMATCH -- do not use")
    return ok


def scotland(lat=56.40, lon=-3.43, name='central Scotland'):
    ts = load.timescale()
    eph = load('de421.bsp')
    earth, sun, moon = eph['earth'], eph['sun'], eph['moon']
    obs = earth + wgs84.latlon(lat, lon)

    e = earth.at(ts.utc(2016, 5, 16, 17))
    elong = e.observe(moon).apparent().separation_from(e.observe(sun).apparent()).degrees
    print(f"\n16 May 2016, {name} ({lat:.2f} N, {abs(lon):.2f} W)")
    print(f"  Moon-Sun elongation {elong:.0f} deg, illuminated fraction "
          f"{(1 - np.cos(D(elong))) / 2 * 100:.0f}%")
    print(f"\n{'BST':>6} {'Moon alt':>9} {'az':>7} {'Sun alt':>9} {'az':>7}"
          f" {'|dphi|':>7} {'sep':>6} {'alpha':>7} {'beta':>7} {'delta':>7}")
    out = []
    for minute in range(14 * 60, 23 * 60 + 1, 30):
        t = ts.utc(2016, 5, 16, 0, minute)
        ma, mz, _ = obs.at(t).observe(moon).apparent().altaz()
        sa, sz, _ = obs.at(t).observe(sun).apparent().altaz()
        nm, ns = ma.degrees, sa.degrees
        if nm < 0 or ns < -0.8:
            continue
        dphi = abs((mz.degrees - sz.degrees + 540) % 360 - 180)
        sep = G(np.arccos(np.sin(D(nm)) * np.sin(D(ns))
                          + np.cos(D(nm)) * np.cos(D(ns)) * np.cos(D(dphi))))
        a, b = alpha(nm, ns, dphi), beta(nm, ns, dphi)
        out.append((minute, a, b))
        print(f"{(minute // 60 + 1) % 24:02d}:{minute % 60:02d}"
              f" {nm:9.1f} {mz.degrees:7.1f} {ns:9.1f} {sz.degrees:7.1f}"
              f" {dphi:7.1f} {sep:6.1f} {a:7.1f} {b:7.1f} {a - b:7.1f}")
    if out:
        A = [r[1] for r in out]
        print(f"\n  alpha runs {max(A):.0f} deg down to {min(A):.0f} deg across the evening.")
        # when does alpha pass through the questioner's own estimate of 45 deg?
        for (m0, a0, _), (m1, a1, _) in zip(out, out[1:]):
            if (a0 - 45) * (a1 - 45) <= 0:
                f = (45 - a0) / (a1 - a0) if a1 != a0 else 0
                mm = m0 + f * (m1 - m0)
                print(f"  it passes 45 deg -- the figure written on the questioner's own "
                      f"diagram -- at about {(int(mm) // 60 + 1) % 24:02d}:{int(mm) % 60:02d} BST.")
                break


if __name__ == '__main__':
    if check():
        for lat, lon, nm in [(56.40, -3.43, 'central Scotland'),
                             (55.95, -3.19, 'Edinburgh'),
                             (57.48, -4.22, 'Inverness')]:
            scotland(lat, lon, nm)