"""Standalone Cartesian Newtonian three-body science model; no browser code."""

import math

INNER_MASS = 3e-5
INNER_PERIOD = 2 * math.pi / math.sqrt(1 + INNER_MASS)
DAY_SCALE = 12 / INNER_PERIOD


def initial_state(settings):
    outer_mass = settings["outerMassEarth"] / 332946.0487
    axis = (
        settings["periodRatio"] * math.sqrt((1 + outer_mass) / (1 + INNER_MASS))
    ) ** (2 / 3)
    phase = settings["outerPhaseDeg"] * math.pi / 180
    eccentricity = settings["innerEccentricity"]
    speed = math.sqrt((1 + outer_mass) / axis)
    masses = [1, INNER_MASS, outer_mass]
    points = [
        [0, 0, 0, 0],
        [
            1 - eccentricity,
            0,
            0,
            math.sqrt((1 + INNER_MASS) * (1 + eccentricity) / (1 - eccentricity)),
        ],
        [
            axis * math.cos(phase),
            axis * math.sin(phase),
            -speed * math.sin(phase),
            speed * math.cos(phase),
        ],
    ]
    center = [
        sum(mass * point[component] for mass, point in zip(masses, points))
        / sum(masses)
        for component in range(4)
    ]
    return masses, [
        value - center[component]
        for point in points
        for component, value in enumerate(point)
    ]


def derivative(state, masses):
    result = []
    for first in range(3):
        horizontal = 0.0
        vertical = 0.0
        for second in range(3):
            if first == second:
                continue
            dx = state[4 * second] - state[4 * first]
            dy = state[4 * second + 1] - state[4 * first + 1]
            factor = masses[second] / math.hypot(dx, dy) ** 3
            horizontal += factor * dx
            vertical += factor * dy
        result.extend(
            [state[4 * first + 2], state[4 * first + 3], horizontal, vertical]
        )
    return result


def rk4(state, masses, step):
    first = derivative(state, masses)
    second = derivative(
        [value + step / 2 * slope for value, slope in zip(state, first)], masses
    )
    third = derivative(
        [value + step / 2 * slope for value, slope in zip(state, second)], masses
    )
    fourth = derivative(
        [value + step * slope for value, slope in zip(state, third)], masses
    )
    return [
        value + step / 6 * (a + 2 * b + 2 * c + d)
        for value, a, b, c, d in zip(state, first, second, third, fourth)
    ]


def interpolate(position, end, velocity, end_velocity, step, fraction):
    return (
        (2 * fraction**3 - 3 * fraction**2 + 1) * position
        + (fraction**3 - 2 * fraction**2 + fraction) * step * velocity
        + (-2 * fraction**3 + 3 * fraction**2) * end
        + (fraction**3 - fraction**2) * step * end_velocity
    )


def record_orbit(settings, resolution=960, periods=200):
    masses, state = initial_state(settings)
    step = INNER_PERIOD / resolution * (1 - settings["innerEccentricity"]) ** 1.5
    time = 0.0
    step_index = 0
    events = [[], []]
    duration = periods * INNER_PERIOD
    maximum_steps = math.ceil(duration / step) + 2
    while time < duration:
        actual = min(step, duration - time)
        following = rk4(state, masses, actual)
        for planet in [1, 2]:
            offset = 4 * planet
            previous_y = state[offset + 1] - state[1]
            next_y = following[offset + 1] - following[1]
            if not previous_y < 0 <= next_y:
                continue
            lower, upper = 0.0, 1.0
            for root_iteration in range(42):
                middle = (lower + upper) / 2
                value = interpolate(
                    previous_y,
                    next_y,
                    state[offset + 3] - state[3],
                    following[offset + 3] - following[3],
                    actual,
                    middle,
                )
                if value < 0:
                    lower = middle
                else:
                    upper = middle
            fraction = (lower + upper) / 2
            front = interpolate(
                state[offset] - state[0],
                following[offset] - following[0],
                state[offset + 2] - state[2],
                following[offset + 2] - following[2],
                actual,
                fraction,
            )
            if front > 0:
                events[planet - 1].append(time + fraction * actual)
        state = following
        step_index += 1
        if step_index > maximum_steps:
            raise RuntimeError(
                "The numerical orbit clock failed to advance to the requested duration"
            )
        time = min(duration, step_index * step)
    return events


def timing_residuals(events):
    center = (len(events) - 1) / 2
    mean = sum(events) / len(events)
    slope = sum(
        (index - center) * (value - mean) for index, value in enumerate(events)
    ) / sum((index - center) ** 2 for index in range(len(events)))
    return [
        value - mean - slope * (index - center) for index, value in enumerate(events)
    ]