"""Post-hoc, and a correction to my own sealed prose.

The sealed rule takes the phase from `ephem.Moon(t).elong` and explains the
departure of the null fraction from 0.5 by the moon's elliptical orbit. The
shape table in results.txt shows the 0-5 degree bin short of its expected
count by a third, which is far too big for an orbital-rate effect, so both
halves of that sentence were checked. Two separate things are wrong with it.

1. `elong` is the ANGULAR SEPARATION OF MOON AND SUN ON THE SKY, not the
   phase angle. The moon's orbit is inclined about 5.1 degrees to the
   ecliptic, so at the instant of new moon the two bodies are usually several
   degrees apart on the sky and the separation almost never reaches 0; the
   same geometry keeps it from reaching 180 at full moon. Time is pushed out
   of the bins nearest syzygy and into the rest, which is the shape table.

2. THE ELLIPTICAL ORBIT IS THE WRONG MECHANISM for the remaining departure
   from 0.5. The equation of the centre -- the big term that comes from the
   ellipse -- has the mean anomaly as its argument, and the anomalistic month
   is not the synodic month, so over fifty-three years it averages out of a
   statistic binned by phase and contributes nothing. The term that cannot
   average out is the one whose argument IS the phase: the Variation, which
   Tycho Brahe found, 0.658 degrees of sin(2D). It makes the moon run ahead
   of its mean place near new and full moon and behind at the quarters, so
   the phase angle sweeps FASTER through syzygy and less time is spent there.
   The function `rate_by_phase` below measures that directly.

Neither touches the verdict, and the reason is the one part of the design that
was right: f_exp is computed by running the SAME function over a time grid, so
whatever elong does, the null expectation does it too. The half-circle split at
d <= 45 is 45 degrees away from the distorted region. What they do mean is that
the sealed rule measured a slightly different thing from what it said, and
explained a real number by a mechanism that is not the one at work.

    python geometry.py > geometry.txt
"""

import math
import pathlib

import ephem
import numpy as np

import controls
import rule

DATA = pathlib.Path.home() / "data" / "quake-moon"


def dlon(times):
    """Difference in ecliptic longitude, moon minus sun: the phase angle.
    0 at new moon, 180 at full. This is the variable the tidal claim is
    actually about."""
    out = np.empty(len(times))
    for i, t in enumerate(times):
        d = ephem.Date(rule._to_pyephem(t))
        m = ephem.Ecliptic(ephem.Moon(d))
        s = ephem.Ecliptic(ephem.Sun(d))
        out[i] = math.degrees(m.lon - s.lon) % 360.0
    return out


def grid_dlon():
    cache = DATA / "grid_hour_dlon.npy"
    if cache.exists():
        return np.load(cache)
    t0 = np.datetime64(rule.START[:-1], "s")
    t1 = np.datetime64(rule.END[:-1], "s")
    n = int((t1 - t0) / np.timedelta64(1, "s")) // 3600 + 1
    g = t0 + (np.arange(n) * 3600).astype("timedelta64[s]")
    th = dlon(g)
    np.save(cache, th)
    return th


def rate_by_phase(g_dl):
    """How fast the phase angle advances, as a function of where it is.

    If the Variation is the mechanism, the sweep is fastest at syzygy and
    slowest at the quarters, and the deficit of syzygy time follows from it.
    """
    step = np.diff(g_dl) % 360.0          # degrees per hour, wrapped
    mid = (g_dl[:-1] + step / 2.0) % 360.0
    d = rule.syzygy_distance(mid)
    edges = np.arange(0, 100, 10)
    idx = np.digitize(d, edges) - 1
    mean_all = step.mean()
    print("\nmean advance of the phase angle, by distance from syzygy")
    print(f"  the mean over the whole window is {mean_all:.5f} deg/hour "
          f"(a synodic month of {360/mean_all/24:.3f} days)")
    for b in range(len(edges) - 1):
        m = step[idx == b].mean()
        print(f"  {edges[b]:>2}-{edges[b+1]:>2} deg from syzygy   "
              f"{m:.5f} deg/h   {100*(m/mean_all - 1):+.3f} % of mean")
    fast = step[d <= 45].mean()
    slow = step[d > 45].mean()
    print(f"  syzygy half {fast:.5f} deg/h against quadrature half "
          f"{slow:.5f}: {100*(fast/slow - 1):+.3f} % faster")
    print(f"  time in the syzygy half predicted from those two rates: "
          f"{(1/fast)/((1/fast)+(1/slow)):.6f}")


def main():
    times, lat, lon, dep, mag = controls.load()
    mains, _ = rule.decluster(times, lat, lon, mag)

    g_el = controls.grid_phases(3600, "hour")
    g_dl = grid_dlon()
    d_el = rule.syzygy_distance(g_el)
    d_dl = rule.syzygy_distance(g_dl)

    print("fraction of 1973-2025 spent within 45 degrees of syzygy")
    print(f"  sky separation (elong, what the seal used)  {(d_el <= 45).mean():.6f}")
    print(f"  ecliptic longitude difference (phase angle) {(d_dl <= 45).mean():.6f}")
    print(f"  exactly half the circle would be            0.500000")

    print("\nwhere the two differ: time per 5-degree bin of syzygy distance")
    print("  bin      elong     phase angle   ratio")
    edges = np.arange(0, 95, 5)
    he, _ = np.histogram(d_el, bins=edges)
    hd, _ = np.histogram(d_dl, bins=edges)
    for i in range(len(he)):
        print(f"  {edges[i]:>2}-{edges[i+1]:>2}  {he[i]/he.sum():.5f}   "
              f"{hd[i]/hd.sum():.5f}     {he[i]/hd[i]:.3f}")

    rate_by_phase(g_dl)

    print("\nthe primary statistic on the phase angle instead of the sky "
          "separation\n  (post-hoc; the sealed verdict stands on the sealed "
          "variable)")
    f_exp = float((d_dl <= 45).mean())
    th = dlon(times[mains])
    rng = np.random.default_rng(rule.SEED + 2)
    e = rule.excess(th, f_exp)
    lo, hi = rule.bootstrap_excess(th, f_exp, rng)
    inside = int((rule.syzygy_distance(th) <= 45).sum())
    print(f"  n = {len(th)}   inside = {inside} ({100*inside/len(th):.2f} %)"
          f"   expected {100*f_exp:.2f} %")
    print(f"  E = {100*e:+.2f} %   95 % [{100*lo:+.2f}, {100*hi:+.2f}]"
          f"   -> {rule.verdict(lo, hi)}")


if __name__ == "__main__":
    main()