"""The reference corpus of doc 40 §5, from Skyfield over JPL DE440 — the truth `:almanac` is tested against.

    python3 tools/corpus/generate.py --kernel ~/Developer/temp/dawngloaming-corpus/de440s.bsp [--only london …]

Run it with exactly the packages in requirements.txt beside it; it refuses anything else, and it refuses a
kernel whose SHA-256 is not the one recorded in Plan/almanac/SOURCES.md. It writes Plan/almanac/corpus/,
one JSON line per place and local day, plus seasons.jsonl and meta.json.

What it computes, and how the app must agree (doc 40 §2, §5.3):
  * Positions are apparent and topocentric, on the WGS84 ellipsoid, with no refraction: every refraction
    and semidiameter lives in the thresholds. The moon's line is −0.5666° − 0.2725 × HP at every instant.
  * An instant is read as UT1 — civil time, which UTC keeps within 0.9 s of — and TT follows from
    Skyfield's built-in ΔT.
  * A day is the half-open [local midnight, next local midnight) of IANA tz 2025b, kept beside the corpus as
    Plan/almanac/tzdb-2025b.bin; the app may bundle a newer release so long as every corpus zone reads the
    same in it (T28). Where a zone springs forward over midnight the day starts where the gap ends, as in
    java.time.
  * Crossings are found on a five-minute grid and bisected to a millisecond. Every extremum of the margin
    that comes within 0.1° of the line is refined as well, so a pair of crossings closer together than the
    grid is still found. A day records its nearest approach as "nearest", and one within 0.001° of the
    line is marked "tangent": its state is not something two ephemerides can be asked to agree on
    (doc 40 §5.5).

Dedicated to the public domain under CC0 1.0 (https://creativecommons.org/publicdomain/zero/1.0/): to the extent
possible under law, Lustig Fairchild Ltd. has waived all copyright and related rights to this script. Published
from Bulgaria, with the almanac it writes (https://lustigfairchild.com/dawn-and-gloaming/almanac/).
"""

import argparse
import datetime as dt
import hashlib
import importlib.metadata
import json
import math
import pathlib
import zoneinfo

import numpy as np

KERNEL_SHA256 = "c1c7feeab882263fc493a9d5a5b2ddd71b54826cdf65d8d17a76126b260a49f2"
PINNED = {"skyfield": "1.55", "numpy": "2.2.6", "jplephem": "2.24", "sgp4": "2.27", "tzdata": "2025.2"}
GENERATOR = "1"

# Doc 40 §5.1. Elevation is carried for La Paz alone; the app computes every place at sea level, and the
# corpus test asserts the difference stays inside the bands.
PLACES = {
    "singapore": (1.3521, 103.8198, "Asia/Singapore", 0),
    "makkah": (21.4225, 39.8262, "Asia/Riyadh", 0),
    "london": (51.5074, -0.1278, "Europe/London", 0),
    "wellington": (-41.2866, 174.7756, "Pacific/Auckland", 0),
    "reykjavik": (64.1466, -21.9426, "Atlantic/Reykjavik", 0),
    "tromso": (69.6492, 18.9553, "Europe/Oslo", 0),
    "mcmurdo": (-77.8419, 166.6863, "Antarctica/McMurdo", 0),
    "la-paz": (-16.4897, -68.1193, "America/La_Paz", 3640),
}
FIXTURES = {
    "oulu-band": ((65.0, 25.0, "Europe/Helsinki", 0), ["2026-07-05", "2026-07-06", "2026-07-07"]),
    "tromso-band": ((69.6846, 25.0, "Europe/Helsinki", 0), ["2026-07-25", "2026-07-26"]),
    "tromso-grazing": ((69.6846, 18.9553, "Europe/Oslo", 0), ["2026-07-25", "2026-07-26"]),
}
SWEPT = {"london", "tromso"}
YEARS = (2026, 2028)
RANGE_DATES = ["1900-01-01", "2000-01-01", "2050-01-01", "2100-12-31"]
SEASON_YEARS = (1900, 2000, 2026, 2027, 2028, 2050, 2100)

SUN = {"horizon": -0.8333, "civil": -6.0, "nautical": -12.0, "astronomical": -18.0, "golden": 6.0, "blue": -4.0}
DEPRESSIONS = [13.0, 15.0, 16.0, 18.2, 18.5, 19.5, 20.0]  # doc 40 §4, less the 12° and 18° lines already here
DEPRESSION_PLACES = {"makkah", "london"}
QUARTERS = ["NEW", "FIRST_QUARTER", "FULL", "LAST_QUARTER"]
SEASONS = ["MARCH_EQUINOX", "JUNE_SOLSTICE", "SEPTEMBER_EQUINOX", "DECEMBER_SOLSTICE"]

GRID = 300.0  # seconds
NEAR = 0.1  # degrees: an extremum of the grid this close to the line is refined
ILL_CONDITIONED_BELOW = 0.02  # degrees a minute (A12)
TANGENT_WITHIN = 0.001  # degrees (doc 40 §5.5)
EARTH_RADIUS_KM = 6378.137
UTC = dt.timezone.utc


def check_environment(kernel):
    for name, version in PINNED.items():
        found = importlib.metadata.version(name)
        if found != version:
            raise SystemExit(f"{name} {found} is installed; the corpus is generated with {version} only")
    digest = hashlib.sha256(kernel.read_bytes()).hexdigest()
    if digest != KERNEL_SHA256:
        raise SystemExit(f"{kernel} has SHA-256 {digest}, not the recorded {KERNEL_SHA256}; refusing it")
    zoneinfo.reset_tzpath(to=[])  # the tzdata package only, never the machine's own tz
    import tzdata

    if tzdata.IANA_VERSION != "2025b":
        raise SystemExit(f"tzdata carries {tzdata.IANA_VERSION}; the corpus is cut with 2025b")


def start_of(date, zone):
    """The first instant of a local day, as java.time's LocalDate.atStartOfDay(ZoneId) gives it."""
    midnight = dt.datetime.combine(date, dt.time(0))
    first = midnight.replace(tzinfo=zone, fold=0).astimezone(UTC)
    second = midnight.replace(tzinfo=zone, fold=1).astimezone(UTC)
    if first <= second or first.astimezone(zone).replace(tzinfo=None) == midnight:
        return min(first, second).timestamp()
    lo, hi = second.timestamp(), first.timestamp()  # midnight falls in a gap: the day starts where it ends

    def offset(s):
        return dt.datetime.fromtimestamp(s, UTC).astimezone(zone).utcoffset()

    before = offset(lo)
    while hi - lo > 0.5:
        mid = (lo + hi) / 2
        lo, hi = (mid, hi) if offset(mid) == before else (lo, mid)
    return float(round(hi))


def days_of(zone, first, last):
    dates = [first + dt.timedelta(days=k) for k in range((last - first).days + 1)]
    starts = [start_of(d, zone) for d in dates + [last + dt.timedelta(days=1)]]
    return dates, starts[:-1], starts[1:]


def iso(seconds):
    return dt.datetime.fromtimestamp(math.floor(seconds + 0.5), UTC).strftime("%Y-%m-%dT%H:%M:%SZ")


class Sky:
    def __init__(self, kernel):
        from skyfield.api import load, wgs84

        self.ts = load.timescale(builtin=True)
        self.eph = load(str(kernel))
        self.earth, self.sun, self.moon = self.eph["earth"], self.eph["sun"], self.eph["moon"]
        self.wgs84 = wgs84

    def times(self, seconds):
        return self.ts.ut1_jd(np.asarray(seconds, dtype=float) / 86400.0 + 2440587.5)

    def seconds(self, t):
        return (t.ut1 - 2440587.5) * 86400.0

    def observer(self, place):
        lat, lon, _, elevation = place
        return self.earth + self.wgs84.latlon(lat, lon, elevation_m=elevation)

    def altaz(self, observer, body, seconds):
        seconds = np.atleast_1d(np.asarray(seconds, dtype=float))
        alt, az, distance = observer.at(self.times(seconds)).observe(body).apparent().altaz()
        return alt.degrees, az.degrees, distance.km

    def ecliptic(self, body, seconds):
        from skyfield.framelib import ecliptic_frame

        position = self.earth.at(self.times(seconds)).observe(body).apparent()
        lat, lon, distance = position.frame_latlon(ecliptic_frame)
        return lon.degrees, lat.degrees, distance.km


def margin(altitude, distance, line):
    """Degrees above the line; a line of None is the moon's, which moves with its distance."""
    if line is None:
        return altitude - (-0.5666 - 0.2725 * np.degrees(np.arcsin(EARTH_RADIUS_KM / distance)))
    return altitude - line


def bisect(f, lo, hi, rising, resolution=0.001):
    lo, hi = np.array(lo, dtype=float), np.array(hi, dtype=float)
    while len(lo) and np.max(hi - lo) > resolution:
        mid = (lo + hi) / 2
        after = (f(mid) > 0) == rising
        lo, hi = np.where(after, lo, mid), np.where(after, mid, hi)
    return (lo + hi) / 2


def golden(f, lo, hi, resolution=0.1):
    """Where f is greatest on each [lo, hi]; f takes and returns arrays."""
    ratio = (math.sqrt(5) - 1) / 2
    lo, hi = np.array(lo, dtype=float), np.array(hi, dtype=float)
    while len(lo) and np.max(hi - lo) > resolution:
        x1, x2 = hi - ratio * (hi - lo), lo + ratio * (hi - lo)
        right = f(x1) < f(x2)
        lo, hi = np.where(right, x1, lo), np.where(right, hi, x2)
    return (lo + hi) / 2


class Crossings:
    """Every crossing of every line of one body over [start, end], found on one shared grid."""

    def __init__(self, sky, observer, body, lines, start, end):
        self.lines = lines
        grid = np.arange(start, end + GRID, GRID)
        grid[-1] = min(grid[-1], end)
        alt, _, distance = sky.altaz(observer, body, grid)
        self.grid = grid
        self.values = {key: margin(alt, distance, line) for key, line in lines.items()}

        keys = list(lines)
        brackets, turns = [], []
        for j, key in enumerate(keys):
            v = self.values[key]
            above = v > 0
            brackets += [(grid[i], grid[i + 1], above[i + 1], j) for i in np.nonzero(above[:-1] != above[1:])[0]]
            dv = np.diff(v)
            turns += [(i, j) for i in np.nonzero((dv[:-1] * dv[1:] <= 0) & (np.abs(v[1:-1]) < NEAR))[0] + 1]

        def f_of(items):
            line = np.array([np.nan if lines[keys[j]] is None else lines[keys[j]] for j in items])
            moon = lines[keys[items[0]]] is None if len(items) else False

            def f(seconds):
                a, _, d = sky.altaz(observer, body, seconds)
                return margin(a, d, None) if moon else a - line

            return f

        self.extrema = {key: [] for key in keys}
        if turns:
            index, which = map(np.array, zip(*turns))
            sign = np.array([1.0 if self.values[keys[j]][i + 1] < self.values[keys[j]][i] else -1.0 for i, j in turns])
            f = f_of(which)
            t = golden(lambda s: sign * f(s), grid[index - 1], grid[index + 1])
            depth = f(t)
            for i, j, te, ve in zip(index, which, t, depth):
                v = self.values[keys[j]]
                self.extrema[keys[j]].append((float(te), float(ve)))
                if (ve > 0) != (v[i] > 0) and (v[i - 1] > 0) == (v[i] > 0) == (v[i + 1] > 0):
                    brackets += [(grid[i - 1], te, ve > 0, j), (te, grid[i + 1], v[i + 1] > 0, j)]

        self.found = {key: [] for key in keys}
        if brackets:
            lo, hi, rising, which = map(np.array, zip(*brackets))
            f = f_of(which)
            at = bisect(f, lo, hi, rising.astype(bool))
            rate = f(at + 30) - f(at - 30)  # degrees over the minute around the crossing
            _, azimuth, _ = sky.altaz(observer, body, at)
            for a, r, q, z, j in zip(at, rising, rate, azimuth, which):
                self.found[keys[j]].append({"at": float(a), "rising": bool(r), "rate": float(q), "az": float(z)})
            for events in self.found.values():
                events.sort(key=lambda e: e["at"])

    def state(self, key, start, end):
        events = [e for e in self.found[key] if start <= e["at"] < end]
        if events:
            out = {
                "state": "CROSSINGS",
                "events": [
                    {
                        "dir": "RISING" if e["rising"] else "SETTING",
                        "at": iso(e["at"]),
                        "certainty": "ILL_CONDITIONED" if abs(e["rate"]) < ILL_CONDITIONED_BELOW else "NORMAL",
                        "rate": round(e["rate"], 4),
                        "az": round(e["az"], 3),
                    }
                    for e in events
                ],
            }
        else:
            inside = np.nonzero((self.grid >= start) & (self.grid < end))[0]
            out = {"state": "ALWAYS_ABOVE" if self.values[key][inside[len(inside) // 2]] > 0 else "ALWAYS_BELOW"}
        near = [v for t, v in self.extrema[key] if start <= t < end]
        if near:
            out["nearest"] = round(min(near, key=abs), 5)
            if abs(out["nearest"]) < TANGENT_WITHIN:
                out["tangent"] = True
        return out, events


def daylight(state, events, start, end):
    """Time the sun's centre spends above the horizon line in [start, end): a sum, so a day with two
    sunsets is well defined, a polar day is the whole day and a polar night none (doc 40 §5.3)."""
    if not events:
        return end - start if state["state"] == "ALWAYS_ABOVE" else 0.0
    edges = ([start] if not events[0]["rising"] else []) + [e["at"] for e in events]
    if events[-1]["rising"]:
        edges.append(end)
    return sum(b - a for a, b in zip(edges[::2], edges[1::2]))


def extremes(sky, observer, starts, ends, sign):
    """The sun's highest (sign 1) or lowest (sign −1) point in each day, found as the app's `highest` and
    `lowest` find it: a ten-minute grid over [start, end − 1 ms], then a golden-section search around the
    best sample. Returns the instants and the altitudes."""
    starts, lasts = np.array(starts), np.array(ends) - 0.001
    samples = np.maximum(2, np.ceil((lasts - starts) / 600).astype(int))
    k = np.minimum(np.arange(samples.max() + 1)[None, :], samples[:, None])
    grid = starts[:, None] + (lasts - starts)[:, None] * k / samples[:, None]
    alt = sky.altaz(observer, sky.sun, grid.ravel())[0].reshape(grid.shape)
    best = np.argmax(sign * alt, axis=1)
    rows = np.arange(len(starts))
    lo, hi = grid[rows, np.maximum(0, best - 1)], grid[rows, np.minimum(samples, best + 1)]
    peak = golden(lambda s: sign * sky.altaz(observer, sky.sun, s)[0], lo, hi)
    at = np.where(sign * sky.altaz(observer, sky.sun, peak)[0] >= sign * alt[rows, best], peak, grid[rows, best])
    return at, sky.altaz(observer, sky.sun, at)[0]


def moon_phase(sky, seconds):
    sun_lon, sun_lat, sun_r = sky.ecliptic(sky.sun, seconds)
    moon_lon, moon_lat, moon_r = sky.ecliptic(sky.moon, seconds)

    def vector(lon, lat, r):
        lon, lat = np.radians(lon), np.radians(lat)
        return np.array([r * np.cos(lat) * np.cos(lon), r * np.cos(lat) * np.sin(lon), r * np.sin(lat)])

    s, m = vector(sun_lon, sun_lat, sun_r), vector(moon_lon, moon_lat, moon_r)
    to_sun = s - m
    cos_i = -np.sum(to_sun * m, axis=0) / (np.linalg.norm(to_sun, axis=0) * np.linalg.norm(m, axis=0))
    return (moon_lon - sun_lon) % 360, (1 + cos_i) / 2


def events_between(sky, function, start, end):
    from skyfield import almanac

    t, kind = almanac.find_discrete(sky.times([start])[0], sky.times([end])[0], function(sky.eph))
    return list(zip(sky.seconds(t), kind))


def seasons(sky, year):
    from skyfield import almanac

    first = dt.datetime(year, 1, 1, tzinfo=UTC).timestamp()
    last = dt.datetime(year + 1, 1, 1, tzinfo=UTC).timestamp()
    return [(SEASONS[k], at) for at, k in events_between(sky, almanac.seasons, first, last)]


def span(sky, name, place, first, last, keep):
    """The rows of one place for the dates in keep, computed over the whole of [first, last]; also the
    state kind of every sun line on every day of it, which is what finds the polar transitions."""
    from skyfield import almanac

    zone = zoneinfo.ZoneInfo(place[2])
    observer = sky.observer(place)
    dates, starts, ends = days_of(zone, first, last)
    sun_lines = dict(SUN)
    if name in DEPRESSION_PLACES:
        sun_lines.update({f"{angle:.1f}": -angle for angle in DEPRESSIONS})
    sun = Crossings(sky, observer, sky.sun, sun_lines, starts[0], ends[-1])
    moon = Crossings(sky, observer, sky.moon, {"moon": None}, starts[0], ends[-1])

    kept = [i for i, d in enumerate(dates) if d in keep]
    kept_starts, kept_ends = [starts[i] for i in kept], [ends[i] for i in kept]
    noon, noon_alt = extremes(sky, observer, kept_starts, kept_ends, 1.0)
    midnight, midnight_alt = extremes(sky, observer, kept_starts, kept_ends, -1.0)
    phase, lit = moon_phase(sky, np.array(kept_starts))
    quarters = events_between(sky, almanac.moon_phases, starts[0], ends[-1] + 40 * 86400)

    rows, kinds, lengths = [], {key: [] for key in SUN}, []
    for i, date in enumerate(dates):
        states = {key: sun.state(key, starts[i], ends[i]) for key in sun_lines}
        lengths.append(daylight(*states["horizon"], starts[i], ends[i]))
        for key in SUN:
            kinds[key].append(states[key][0]["state"])
        if date not in keep:
            continue
        k = kept.index(i)
        first_day = date == dt.date(1900, 1, 1)  # its yesterday is outside the range (A19)
        row = {
            "place": name,
            "date": date.isoformat(),
            "zone": place[2],
            "lat": place[0],
            "lon": place[1],
            "day": {"start": iso(starts[i]), "end": iso(ends[i])},
            "sun": {key: states[key][0] for key in SUN},
            "moon": {
                "horizon": moon.state("moon", starts[i], ends[i])[0],
                "phase_deg": round(float(phase[k]), 3),
                "illuminated": round(float(lit[k]), 5),
                "next_phases": [{"phase": QUARTERS[q], "at": iso(t)} for t, q in quarters if t >= starts[i]][:4],
            },
            "solar": {
                "noon": iso(noon[k]),
                "noon_alt": round(float(noon_alt[k]), 4),
                "midnight": iso(midnight[k]),
                "midnight_alt": round(float(midnight_alt[k]), 4),
                "day_length_s": round(lengths[-1]),
                "delta_yesterday_s": None if first_day or i == 0 else round(lengths[-1] - lengths[-2]),
            },
        }
        if name in DEPRESSION_PLACES:
            row["depression"] = {f"{a:.1f}": states[f"{a:.1f}"][0] for a in DEPRESSIONS}
        rows.append(row)
    return rows, dates, kinds


def selection(sky, name, place):
    """Doc 40 §5.2's rows for one place: each year is computed whole, and the dates it asks for kept."""
    zone = zoneinfo.ZoneInfo(place[2])
    rows = []
    for year in YEARS:
        first, last = dt.date(year - 1, 12, 31), dt.date(year, 12, 31)
        dates, starts, ends = days_of(zone, first, last)
        every, _, kinds = span(sky, name, place, first, last, keep=set(dates[1:]))
        named = {d for d, start, end in zip(dates, starts, ends) if end - start != 86400}  # every DST change
        for series in kinds.values():
            for i in range(1, len(series)):
                if series[i] != series[i - 1]:  # every change of state, ±3 days
                    named.update(dates[i] + dt.timedelta(days=k) for k in range(-3, 4))
        named.update(dt.datetime.fromtimestamp(at, UTC).astimezone(zone).date() for _, at in seasons(sky, year))
        for row in every:
            date = dt.date.fromisoformat(row["date"])
            if date in named:
                rows.append(row)
            elif name in SWEPT and date.year == year:
                row.pop("depression", None)  # the sweep is for the horizon and twilights; the angles have enough
                rows.append(row)
    for text in RANGE_DATES:
        d = dt.date.fromisoformat(text)
        rows += span(sky, name, place, d - dt.timedelta(days=1), d, {d})[0]
    return rows


def write(path, rows):
    rows.sort(key=lambda r: r["date"])
    with open(path, "w", encoding="utf-8", newline="\n") as f:
        for row in rows:
            f.write(json.dumps(row, ensure_ascii=False, separators=(",", ":")) + "\n")
    return len(rows)


def main():
    parser = argparse.ArgumentParser(description=__doc__.splitlines()[0])
    parser.add_argument("--kernel", required=True, type=pathlib.Path)
    parser.add_argument("--out", type=pathlib.Path, default=pathlib.Path(__file__).parents[2] / "Plan" / "almanac")
    parser.add_argument("--only", nargs="+", choices=[*PLACES, *FIXTURES], help="regenerate these; keep the rest")
    args = parser.parse_args()
    check_environment(args.kernel)
    sky = Sky(args.kernel)

    corpus = args.out / "corpus"
    corpus.mkdir(parents=True, exist_ok=True)
    for name in args.only or [*PLACES, *FIXTURES]:
        if name in PLACES:
            rows = selection(sky, name, PLACES[name])
        else:
            place, texts = FIXTURES[name]
            days = [dt.date.fromisoformat(d) for d in texts]
            rows = span(sky, name, place, days[0] - dt.timedelta(days=1), days[-1], set(days))[0]
        print(name, write(corpus / f"{name}.jsonl", rows), flush=True)

    with open(args.out / "seasons.jsonl", "w", encoding="utf-8", newline="\n") as f:
        for year in SEASON_YEARS:
            for kind, at in seasons(sky, year):
                f.write(json.dumps({"year": year, "kind": kind, "at": iso(at)}) + "\n")

    import skyfield

    iers = pathlib.Path(skyfield.__file__).parent / "data" / "iers.npz"
    meta = {
        "generator": GENERATOR,
        "packages": PINNED,
        "kernel": "de440s.bsp",
        "kernel_sha256": KERNEL_SHA256,
        "delta_t": f"Skyfield built-in; iers.npz SHA-256 {hashlib.sha256(iers.read_bytes()).hexdigest()}",
        "tz": "IANA 2025b",
        "rows": {p.stem: len(p.read_text(encoding="utf-8").splitlines()) for p in sorted(corpus.glob("*.jsonl"))},
    }
    (args.out / "meta.json").write_text(json.dumps(meta, indent=2) + "\n", encoding="utf-8")


if __name__ == "__main__":
    main()