orbital-test/soi_calculator.py

"""
SOI (Sphere of Influence) crossing-time calculator.

Given two bodies in Keplerian orbits around the same parent,
find all times when the distance between them equals a given
SOI radius — i.e., when one body enters or exits the other's
gravitational sphere of influence.

Algorithm:
  1. Quick radial-overlap check (rejects impossible cases).
  2. Coarse time-scan across the search window.
  3. Newton-Raphson refinement of each sign-change bracket.
"""

import math
from typing import List, Optional, Tuple

from orbital_elements import OrbitalElements


def find_soi_crossings(
    body: OrbitalElements,
    planet: OrbitalElements,
    soi_radius: float,
    window_start: float,
    window_end: float,
    n_steps: int = 500,
) -> List[Tuple[float, float]]:
    """
    Find all (enter_time, exit_time) SOI crossing pairs.

    Parameters
    ----------
    body, planet : OrbitalElements
        Two bodies orbiting the same central body (same mu).
    soi_radius : float
        Radius of the planet's sphere of influence.
    window_start, window_end : float
        Search window in simulation time.
    n_steps : int
        Number of coarse-scan steps (default 500).

    Returns
    -------
    List of (enter, exit) time pairs, sorted chronologically.
    Empty list if no crossings found.
    """
    # ── Quick rejection: orbital radius ranges must overlap ──
    b_min, b_max = body.radius_range()
    p_min, p_max = planet.radius_range()

    if b_max + soi_radius < p_min or p_max + soi_radius < b_min:
        return []  # orbits never come close enough

    # ── Coarse time scan ────────────────────────────────────
    duration = window_end - window_start
    dt = duration / n_steps

    crossings: List[Tuple[float, float]] = []
    prev_dist: Optional[float] = None
    enter_bracket: Optional[Tuple[float, float]] = None

    for i in range(n_steps + 1):
        t = window_start + i * dt

        bx, by = body.position_at(t)
        px, py = planet.position_at(t)
        dist = math.hypot(bx - px, by - py)

        if prev_dist is not None:
            prev_t = t - dt
            was_outside = prev_dist > soi_radius
            is_inside = dist <= soi_radius

            if was_outside and is_inside:
                # ── Entering SOI — save narrow bracket ──
                enter_bracket = (prev_t, t)

            elif not was_outside and not is_inside and enter_bracket is not None:
                # ── Exiting SOI — refine both crossings ──
                exit_bracket = (prev_t, t)
                enter_t = _refine(
                    body=body, planet=planet, soi_radius=soi_radius,
                    bracket=enter_bracket, entering=True,
                )
                exit_t = _refine(
                    body=body, planet=planet, soi_radius=soi_radius,
                    bracket=exit_bracket, entering=False,
                )
                crossings.append((enter_t, exit_t))
                enter_bracket = None

        prev_dist = dist

    # ── Still inside SOI at window end? ─────────────────────
    if enter_bracket is not None:
        t_end = window_end
        exit_bracket = (window_end - dt, window_end)
        enter_t = _refine(
            body=body, planet=planet, soi_radius=soi_radius,
            bracket=enter_bracket, entering=True,
        )
        exit_t = _refine(
            body=body, planet=planet, soi_radius=soi_radius,
            bracket=exit_bracket, entering=False,
        )
        crossings.append((enter_t, exit_t))

    return crossings


# ── Newton-Raphson refinement ───────────────────────────────────


def _refine(
    body: OrbitalElements,
    planet: OrbitalElements,
    soi_radius: float,
    bracket: Tuple[float, float],
    entering: bool,
) -> float:
    """
    Newton-Raphson refinement of an SOI crossing time.

    Solves  |r_body(t) − r_planet(t)| − soi_radius = 0.

    Parameters
    ----------
    bracket : (t_low, t_high)
        Narrow bracket around the crossing (ideally one scan-step wide).
    entering : bool
        True  → distance *drops* to SOI (outer → inner).
        False → distance *rises*  to SOI (inner → outer).
    """
    t_low, t_high = bracket

    # Start from the side closest to the expected root.
    # For entering: outside edge (t_low, where dist > soi).
    # For exiting:  outside edge (t_high, where dist > soi).
    t = t_low if entering else t_high

    newton_tol = 1e-8
    max_iter = 30

    for _ in range(max_iter):
        bx, by, bvx, bvy = body.compute_state_at(t)
        px, py, pvx, pvy = planet.compute_state_at(t)

        dx = bx - px
        dy = by - py
        dvx = bvx - pvx
        dvy = bvy - pvy

        dist = math.hypot(dx, dy)
        f_val = dist - soi_radius

        if dist < 1e-14:
            break  # co-located — degenerate case

        fprime = (dx * dvx + dy * dvy) / dist

        if abs(fprime) < 1e-14:
            break  # stationary w.r.t. SOI boundary

        dt_correction = -f_val / fprime
        t += dt_correction

        # Stay inside bracket
        t = max(t_low, min(t_high, t))

        if abs(dt_correction) < newton_tol:
            break

    return t


# ── Convenience: next SOI event ──────────────────────────────────


def find_next_soi_event(
    body: OrbitalElements,
    planet: OrbitalElements,
    soi_radius: float,
    from_time: float,
    max_orbits: int = 3,
) -> Optional[Tuple[float, float]]:
    """
    Find the first (enter, exit) SOI crossing after *from_time*.

    Searches *max_orbits* × the longer of the two orbital periods.
    Returns None if no crossing found in that window.
    """
    window = max_orbits * max(body.period, planet.period)
    crossings = find_soi_crossings(
        body, planet, soi_radius, from_time, from_time + window,
    )
    return crossings[0] if crossings else None