#!/usr/bin/env python3
"""Independent full-force and geometric-root validation; optional REBOUND5.2.2.

Input is public Cartesian initial states, never executable Engine code. Roots
reintegrate copies of an IAS15 bracket; no RK4/Hermite interpolation is reused.
"""
from pathlib import Path
import argparse
import json
import sys
import math


def main():
    parser = argparse.ArgumentParser()
    parser.add_argument("--input", type=Path, required=True)
    parser.add_argument("--out", type=Path, required=True)
    parser.add_argument("--rebound-path", type=Path)
    arguments = parser.parse_args()
    if arguments.rebound_path:
        sys.path.insert(0, str(arguments.rebound_path))
    import rebound

    assert rebound.__version__ == "5.2.2"
    source = json.loads(arguments.input.read_text())
    output = {
        "reference": "REBOUND5.2.2 IAS15 epsilon1e-12; independent bracketed minimum projected-distance roots",
        "coordinate": "observer positive x, sky y; astrocentric differences; geometric time",
        "endpoint_policy": "compare interior roots at least0.001 initial period from the requested endpoint",
        "cases": [],
    }
    for case in source["cases"]:
        simulation = rebound.Simulation()
        simulation.G = 1
        simulation.integrator = "ias15"
        simulation.integrator.epsilon = 1e-12
        for body in case["initial"]:
            simulation.add(
                m=body["mass"], x=body["x"], y=body["y"], vx=body["vx"], vy=body["vy"]
            )
        events = [[], []]

        def derivative(state, planet):
            star, body = state.particles[0], state.particles[planet]
            return (body.y - star.y) * (body.vy - star.vy)

        for index in range(8000):
            before = simulation.copy()
            target = case["duration"] * (index + 1) / 8000
            simulation.integrate(target)
            for planet in [1, 2]:
                star, body = simulation.particles[0], simulation.particles[planet]
                if not (
                    derivative(before, planet) < 0 <= derivative(simulation, planet)
                    and body.x > star.x
                ):
                    continue
                lower, upper = before.t, simulation.t
                for iteration in range(35):
                    middle = (lower + upper) / 2
                    trial = before.copy()
                    trial.integrate(middle)
                    if derivative(trial, planet) < 0:
                        lower = middle
                    else:
                        upper = middle
                value = (lower + upper) / 2
                if value < case["duration"] - case["period"] * 0.001:
                    events[planet - 1].append(value)
        differences = []
        for reference, original in zip(events, case["events"]):
            kept = [
                value
                for value in original
                if value < case["duration"] - case["period"] * 0.001
            ]
            assert len(reference) == len(kept), (
                case["name"],
                len(reference),
                len(kept),
            )
            differences.append(
                max(abs(a - b) for a, b in zip(reference, kept))
                * 12
                / case["period"]
                * 86400
            )
        output["cases"].append(
            {
                "name": case["name"],
                "settings": case["settings"],
                "events": events,
                "maximum_timing_error_seconds": differences,
                "rk4_relative_energy_error": case["energyError"],
                "rk4_relative_angular_momentum_error": case["angularMomentumError"],
            }
        )
        print(case["name"], differences, flush=True)
    arguments.out.write_text(json.dumps(output, indent=2) + "\n")


if __name__ == "__main__":
    main()