"""
Orbit trajectory point generation for rendering.

generate_orbit_points() — full orbit (elliptical) or truncated
  hyperbolic path within a given cut-off radius.
"""

import math
from typing import List, Optional, Tuple

from orbital_elements import OrbitalElements, TAU


def generate_orbit_points(
    orbit: OrbitalElements,
    num_points: int = 200,
    max_radius: Optional[float] = None,
) -> List[Tuple[float, float]]:
    """
    Generate *num_points* world-space positions along the orbit.

    Elliptical orbits → full closed loop, equally spaced in E.
    Hyperbolic orbits → samples H ∈ [−H_max, +H_max], truncated
      at *max_radius* (e.g. the planet's SOI edge).

    Parameters
    ----------
    max_radius : float or None
        Only used for hyperbolic orbits.  Points with r > max_radius
        are omitted.  If None, samples out to r ≈ 5·|a| from focus.
    """
    if orbit.is_hyperbolic:
        return _hyperbolic_points(orbit, num_points, max_radius)
    return _elliptical_points(orbit, num_points)


def _elliptical_points(
    orbit: OrbitalElements, num_points: int,
) -> List[Tuple[float, float]]:
    points = []
    for i in range(num_points):
        E = TAU * i / num_points
        points.append(orbit.position_from_E(E))
    return points


def _hyperbolic_points(
    orbit: OrbitalElements,
    num_points: int,
    max_radius: Optional[float],
) -> List[Tuple[float, float]]:
    """
    Sample the hyperbolic branch that passes through periapsis.

    Truncates at *max_radius* by solving for the H values where
    r = max_radius, then sampling uniformly between those limits.
    """
    r_limit = max_radius if max_radius is not None else 5.0 * abs(orbit.a)

    # Find H where r = r_limit.
    # r = |a|(e·cosh(H) − 1)  →  cosh(H) = (r/|a| + 1) / e
    abs_a = abs(orbit.a)
    cosh_Hmax = (r_limit / abs_a + 1.0) / orbit.e
    if cosh_Hmax < 1.0:
        return []  # entire orbit is outside the radius
    H_max = math.acosh(cosh_Hmax)

    points = []
    for i in range(num_points):
        frac = -1.0 + 2.0 * i / (num_points - 1) if num_points > 1 else 0.0
        H = frac * H_max
        x, y = orbit.position_from_H(H)
        r = math.hypot(x, y)
        if r <= r_limit + 1e-9:
            points.append((x, y))

    return points