def rhs(t: float, y: np.ndarray) -> np.ndarray:
d = np.zeros(dim)
# ── Per-target coordinated-turn dynamics ─────────────────
positions = np.empty((N_targets, 3))
for i in range(N_targets):
b = i * _STATES_PER_TARGET
vx, vy, vz = y[b + 3], y[b + 4], y[b + 5]
omega = y[b + 6]
d[b + 0] = vx
d[b + 1] = vy
d[b + 2] = vz
d[b + 3] = -omega * vy
d[b + 4] = omega * vx
d[b + 5] = 0.0
d[b + 6] = 0.0
positions[i, 0] = y[b + 0]
positions[i, 1] = y[b + 1]
positions[i, 2] = y[b + 2]
# Process noise (interpolated)
ch = 4 * i
d[b + 3] += np.interp(t, noise_ts, noise_table[ch + 0])
d[b + 4] += np.interp(t, noise_ts, noise_table[ch + 1])
d[b + 5] += np.interp(t, noise_ts, noise_table[ch + 2])
d[b + 6] += sigma_omega * np.interp(t, noise_ts, noise_table[ch + 3])
# ── Proximity-based dynamic coupling ──────────────────────
for i in range(N_targets):
bi = i * _STATES_PER_TARGET
for j in range(i + 1, N_targets):
bj = j * _STATES_PER_TARGET
dx = positions[j, 0] - positions[i, 0]
dy = positions[j, 1] - positions[i, 1]
dz = positions[j, 2] - positions[i, 2]
dist_sq = dx * dx + dy * dy + dz * dz
if dist_sq < gate_sq:
# Soft association weight (Gaussian decay within gate)
w = coupling_strength * np.exp(-0.5 * dist_sq / (gate_sq * 0.25))
inv_dist = 1.0 / np.sqrt(dist_sq + _EPS)
ux, uy, uz = dx * inv_dist, dy * inv_dist, dz * inv_dist
# Symmetric velocity perturbation (proximity-induced track confusion)
d[bi + 3] += w * ux
d[bi + 4] += w * uy
d[bi + 5] += w * uz
d[bj + 3] -= w * ux
d[bj + 4] -= w * uy
d[bj + 5] -= w * uz
# Cross-coupling bleeds into turn-rate estimates
bearing_ij = np.arctan2(dy, dx + _EPS)
d[bi + 6] += 0.1 * w * np.sin(bearing_ij)
d[bj + 6] -= 0.1 * w * np.sin(bearing_ij)
return d