"""How strong are the planets' tides on Earth, next to the Moon's and the Sun's?

Standard library only. The peak tidal (differential) acceleration a body of mass M at distance d
exerts across Earth's radius R, along the line to the body, is 2 G M R / d^3 (the leading term
of the tidal expansion). GM values: JPL, DE440 (https://ssd.jpl.nasa.gov/astro_par.html);
R = 6371.0084 km (JPL planetary physical parameters). Distances: our Horizons tables
(heliocentric states; geocentric vector = body - Earth), hourly, for every year we fetched.

Usage: python3 -I tides.py <horizons dir> [out.json]
"""
import json, math, pathlib, sys

AU = 149597870.700                       # km
GM = {"Sun": 1.32712440041279419e11, "Moon": 4902.800118, "Mercury": 22031.868551,
      "Venus": 324858.592000, "Mars": 42828.375816, "Jupiter": 126712764.1}   # km^3/s^2
CODE = {"Mercury": "199", "Venus": "299", "Moon": "301", "Mars": "499", "Jupiter": "599"}
R = 6371.0084                            # km


def positions(path):
    text = pathlib.Path(path).read_text(encoding="utf-8")
    body = text.split("$$SOE", 1)[1].split("$$EOE", 1)[0]
    out = {}
    for line in body.strip().splitlines():
        f = [x.strip() for x in line.split(",")]
        out[f[0]] = tuple(float(v) * AU for v in f[2:5])
    return out


def tide(gm, d_km):
    return 2 * gm * R / d_km**3 * 1e3     # m/s^2 (km^3/s^2 * km / km^3 = km/s^2 -> x1000)


def main():
    hdir = pathlib.Path(sys.argv[1])
    years = sorted({p.name.split("_")[1][:4] for p in hdir.glob("399_*_1h.txt")})
    dist = {k: [math.inf, 0.0] for k in ["Sun"] + list(CODE)}
    moon_d = {}                          # Earth-Moon distance by Julian date, for the month-by-month swing
    for f in sorted(hdir.glob("399_*_1h.txt")):
        earth = positions(f)
        for t, e in earth.items():
            d = math.dist(e, (0, 0, 0)); dist["Sun"][0] = min(dist["Sun"][0], d); dist["Sun"][1] = max(dist["Sun"][1], d)
        for name, code in CODE.items():
            g = hdir / f.name.replace("399_", code + "_", 1)
            if not g.exists():
                continue
            other = positions(g)
            for t, e in earth.items():
                d = math.dist(other[t], e)
                dist[name][0] = min(dist[name][0], d); dist[name][1] = max(dist[name][1], d)
                if name == "Moon":
                    moon_d[float(t)] = d
    span = f"{int(min(years)) + 1}-{int(max(years)) + 1}" if years else "?"   # file for year Y starts Y-1-12-31
    moon_mean = tide(GM["Moon"], 384400.0)
    print(f"Distances from Earth's centre over {span} (hourly), and peak tidal acceleration 2GMR/d^3")
    print(f"Moon at its mean distance (384,400 km): {moon_mean:.3e} m/s^2 = 1 (reference)")
    for name, (dmin, dmax) in dist.items():
        hi, lo = tide(GM[name], dmin), tide(GM[name], dmax)
        print(f"{name:8s} d = {dmin/AU:.4f}-{dmax/AU:.4f} au ({dmin:,.0f}-{dmax:,.0f} km): "
              f"tide {lo:.3e}-{hi:.3e} m/s^2 = {lo/moon_mean:.3g}-{hi/moon_mean:.3g} x Moon(mean)")
    planets = sum(tide(GM[n], dist[n][0]) for n in ("Mercury", "Venus", "Mars", "Jupiter"))
    moon_swing = tide(GM["Moon"], dist["Moon"][0]) - tide(GM["Moon"], dist["Moon"][1])
    print(f"All four planets at their closest at once: {planets:.3e} m/s^2 = {planets/moon_mean:.2e} x Moon(mean)")
    print(f"Swing of the Moon's own tide, perigee to apogee: {moon_swing:.3e} m/s^2, "
          f"{moon_swing/planets:,.0f} times the planets' combined maximum")
    # month by month: each perigee (local minimum of the distance) against the apogee that follows it
    ts = sorted(moon_d); ds = [moon_d[t] for t in ts]
    ext = [(i, "p" if ds[i] < ds[i - 1] and ds[i] <= ds[i + 1] else "a")
           for i in range(1, len(ds) - 1) if (ds[i] < ds[i - 1] and ds[i] <= ds[i + 1]) or (ds[i] > ds[i - 1] and ds[i] >= ds[i + 1])]
    swings = [(tide(GM["Moon"], ds[i]) - tide(GM["Moon"], ds[j])) / moon_mean
              for (i, ki), (j, kj) in zip(ext, ext[1:]) if ki == "p" and kj == "a"]
    sw = sorted(swings); med = sw[len(sw) // 2] if len(sw) % 2 else (sw[len(sw) // 2 - 1] + sw[len(sw) // 2]) / 2
    print(f"Month by month ({len(sw)} perigee-to-apogee swings): {sw[0]:.3f}-{sw[-1]:.3f} of the Moon's mean tide, "
          f"median {med:.3f}; {sw[0] * moon_mean / planets:,.0f}-{sw[-1] * moon_mean / planets:,.0f} times the planets' maximum")
    if len(sys.argv) > 2:                  # optional: the same numbers as JSON
        rows = [{"body": n, "closest_au": round(d0 / AU, 4), "farthest_au": round(d1 / AU, 4),
                 "closest_km": round(d0), "farthest_km": round(d1),
                 "peak_tide_m_s2": float(f"{tide(GM[n], d0):.4g}"), "least_tide_m_s2": float(f"{tide(GM[n], d1):.4g}"),
                 "peak_vs_moon_mean": float(f"{tide(GM[n], d0) / moon_mean:.3g}"),
                 "least_vs_moon_mean": float(f"{tide(GM[n], d1) / moon_mean:.3g}")} for n, (d0, d1) in dist.items()]
        out = {"span": span, "moon_mean_tide_m_s2": float(f"{moon_mean:.4g}"), "bodies": rows,
               "four_planets_closest_m_s2": float(f"{planets:.4g}"),
               "four_planets_vs_moon_mean": float(f"{planets / moon_mean:.3g}"),
               "moon_swing_m_s2": float(f"{moon_swing:.4g}"),
               "moon_swing_vs_moon_mean": float(f"{moon_swing / moon_mean:.3g}"),
               "moon_swing_over_four_planets": round(moon_swing / planets),
               "monthly_swings": {"months": len(sw), "least_vs_moon_mean": round(sw[0], 3),
                                  "median_vs_moon_mean": round(med, 3), "largest_vs_moon_mean": round(sw[-1], 3),
                                  "least_over_four_planets": round(sw[0] * moon_mean / planets),
                                  "largest_over_four_planets": round(sw[-1] * moon_mean / planets)}}
        with open(sys.argv[2], "w", encoding="utf-8", newline="\n") as f:
            json.dump(out, f, indent=1); f.write("\n")


if __name__ == "__main__":
    main()