Reference-Frame Integration for Ships (#960)

Supersedes: midpoint-body-integration.md (#958 midpoint fix, confirmed no-op)

Problem

Ships orbiting moons lose orbital energy due to numerical tidal drift. The leapfrog integrator operates in ICRF where the dominant acceleration (from the parent planet) is almost entirely common-mode — it accelerates both ship and moon equally. The tidal residual that determines orbital evolution is a small difference of two large accelerations, amplifying numerical error.

Investigation (#958 comment) confirmed:

Midpoint body positions are equivalent to Strang splitting (no-op fix)
2-body integration is perfectly symplectic (zero drift)
3-body drift is physical tidal perturbation, not a timestep bug
Eccentric orbits show mild dt-dependence (2.4x at dt=0.02 vs dt=1.0) because smaller dt better resolves perigee tidal effects

Solution

Integrate ships in the reference body’s frame (Encke’s method). Instead of tracking absolute ICRF position/velocity, integrate the relative state (delta_r, delta_v) where:

delta_r = r_ship - r_ref
delta_v = v_ship - v_ref

The equation of motion in the relative frame is:

delta_a = -mu_ref / |delta_r|^3 * delta_r + a_tidal

where a_tidal is the tidal acceleration from all other bodies:

a_tidal = sum_i [ G*M_i * (r_i - r_ship)/|r_i - r_ship|^3
                - G*M_i * (r_i - r_ref)/|r_i - r_ref|^3 ]

This is mathematically exact — no terms dropped. The numerical advantage is that leapfrog integrates the dominant force (reference body gravity) directly, while the tidal perturbation enters as a small correction.

Design

Reference body selection

Reuse the existing find_reference_body() from attitude.py (Hill sphere SOI containment). This is already computed per tick using ref_bodies (pre-update positions). If no reference body is found (ship in deep space), fall back to standard ICRF integration.

Tidal acceleration

New function compute_tidal_acceleration() in nbody.py:

def compute_tidal_acceleration(
    ship_pos: np.ndarray,
    ref_body_pos: np.ndarray,
    ref_body_mass: float,
    bodies: list[CelestialBody],
    ref_body_name: str,
) -> np.ndarray:

For each body that is NOT the reference body:

Compute gravity at ship position: G*M*(r_body - r_ship) / |r_body - r_ship|^3
Compute gravity at ref body position: G*M*(r_body - r_ref) / |r_body - r_ref|^3
Tidal contribution = (1) - (2)

Sum all tidal contributions. This is the perturbation that leapfrog must track, and it is orders of magnitude smaller than the direct terms.

Ship integration (_update_ship)

Find reference body from ref_bodies (pre-update positions for Hill sphere consistency)
If reference body found: a. Compute relative state: delta_r = ship.pos - ref.pos, delta_v = ship.vel - ref.vel b. First kick: a1 = -mu/|delta_r|^3 * delta_r + a_tidal(ship_pos) + a_thrust + a_drag c. delta_v_half = delta_v + 0.5 * dt * a1 d. delta_r_new = delta_r + dt * delta_v_half e. Second kick: recompute gravity and tidal at new relative position f. delta_v_new = delta_v_half + 0.5 * dt * a2 g. Recover ICRF: ship.pos = ref_body_post.pos + delta_r_new, ship.vel = ref_body_post.vel + delta_v_new
If no reference body: use current ICRF integration (fallback)

Reference body positions for recovery

Step (g) uses the post-update reference body position (from bodies, which have been integrated by update_bodies_compute before ships). This ensures the ship’s ICRF position is consistent with the current body positions at end-of-tick.

The tidal acceleration in steps (b) and (e) uses the body positions from mid_bodies (or ref_bodies/bodies as fallback). The choice between mid/ref/post for the tidal term has negligible impact since the tidal term itself is small.

Thrust, drag, RCS

These are proper accelerations (not gravitational). They are the same in any frame. Add them to the relative acceleration just as in ICRF.

Attitude control

Unchanged. Attitude control already uses the reference body for direction computation. No frame dependence.

Station/jumpgate integration

Stations and jump gates also benefit from reference-frame integration. Apply the same approach in update_station() in nbody.py.

Cleanup

Remove midpoint body computation from process_tick() since it is a confirmed no-op. Remove mid_bodies parameters from _update_ship() and update_station(). The tidal acceleration computation replaces the midpoint approach with a physically correct solution.

Keep ref_bodies (pre-update snapshot) — still needed for:

Hill sphere reference body lookup (attitude + integration)
ICRF recovery (post-update bodies used directly)

Files modified

services/game-engine/src/physics/nbody.py — add compute_tidal_acceleration(), update update_station()
services/game-engine/src/physics/simulation.py — refactor _update_ship() for relative-frame integration, remove midpoint computation from process_tick()
services/game-engine/tests/physics/test_simulation_extra.py — update #958 regression tests
services/game-engine/tests/physics/test_958_diag.py — update diagnostic tests

Verification

2-body energy conservation: 0% drift (unchanged)
3-body eccentric Io orbit: drift < 0.005%/orbit at dt=0.02 (>7x improvement over current 0.036%/orbit)
All existing tests pass
External interfaces unchanged (Redis, API, rendering all use ICRF)