"""The study's sample and its reference longitudes, from NASA JPL's DE440s.

Horospire ephemeris accuracy study, version 1.0. MIT licence (LICENSE-CODE.txt).

Usage:
    python reference.py DE440S_BSP OUT_DIR

Writes OUT_DIR/moments.json (the civil times, one a day) and
OUT_DIR/reference.json (ten longitudes per moment, degrees).

The sample: one moment a day from 1900-01-01 to 2050-12-31 (55,152 moments).
Day i is read at hour (7 i mod 24) and minute (13 i mod 60), so the sample
walks through every hour of the day instead of sitting on noon.

Time: a civil time before 1972 is read as UT1, because a clock time of that era
was kept in GMT; from 1972 it is read as UTC. Skyfield applies Delta T.

Reference: apparent geocentric ecliptic longitude on the true ecliptic and
equinox of date (Skyfield's ecliptic_frame): light time, annual aberration and
nutation applied, no topocentric parallax. Mars to Pluto are their system
barycentres, which is what DE440s carries.

Needs: skyfield 1.55, jplephem 2.24, numpy, and de440s.bsp from
https://ssd.jpl.nasa.gov/ftp/eph/planets/bsp/de440s.bsp
(sha256 c1c7feeab882263fc493a9d5a5b2ddd71b54826cdf65d8d17a76126b260a49f2).
"""
import json
import os
import sys
from datetime import datetime, timedelta, timezone

from skyfield.api import load, load_file
from skyfield.framelib import ecliptic_frame

BODIES = ["sun", "moon", "mercury", "venus", "mars", "jupiter", "saturn", "uranus", "neptune", "pluto"]
TARGETS = ["sun", "moon", "mercury", "venus", "mars barycenter", "jupiter barycenter",
           "saturn barycenter", "uranus barycenter", "neptune barycenter", "pluto barycenter"]


def moments():
    out = []
    day = datetime(1900, 1, 1, tzinfo=timezone.utc)
    end = datetime(2051, 1, 1, tzinfo=timezone.utc)
    i = 0
    while day < end:
        out.append(day.replace(hour=(i * 7) % 24, minute=(i * 13) % 60))
        day += timedelta(days=1)
        i += 1
    return out


def main():
    bsp, out_dir = sys.argv[1:3]
    os.makedirs(out_dir, exist_ok=True)
    eph = load_file(bsp)
    ts = load.timescale(builtin=True)
    earth = eph["earth"]
    when = moments()

    early = [w for w in when if w.year < 1972]
    late = [w for w in when if w.year >= 1972]
    parts = []
    if early:
        parts.append(ts.ut1([w.year for w in early], [w.month for w in early], [w.day for w in early],
                            [w.hour for w in early], [w.minute for w in early], [w.second for w in early]))
    if late:
        parts.append(ts.from_datetimes(late))

    columns = [[] for _ in BODIES]
    delta_t = []
    for t in parts:
        observer = earth.at(t)
        for k, name in enumerate(TARGETS):
            _, lon, _ = observer.observe(eph[name]).apparent().frame_latlon(ecliptic_frame)
            columns[k].extend(round(float(x) % 360.0, 7) for x in lon.degrees)
        delta_t.extend(round(float(x), 3) for x in (t.tt - t.ut1) * 86400.0)

    json.dump([w.strftime("%Y-%m-%dT%H:%M:%SZ") for w in when], open(os.path.join(out_dir, "moments.json"), "w"))
    json.dump({"bodies": BODIES, "lon": [list(row) for row in zip(*columns)], "delta_t_seconds": delta_t,
               "time_scale": ["UT1" if w.year < 1972 else "UTC" for w in when]},
              open(os.path.join(out_dir, "reference.json"), "w"))
    print(f"{len(when)} moments, {len(when) * len(BODIES)} placements")


if __name__ == "__main__":
    main()