Q-law / Lyapunov Steering for Rendezvous

Overview

Q-law is a Lyapunov feedback control law that targets orbital element differences directly, computing optimal thrust direction each tick via Gauss variational equations (GVE). Unlike Lambert+VTG which assumes impulsive thrust, Q-law naturally handles continuous-thrust vehicles without burn/coast heuristics.

The rendezvous system uses the "manual" strategy exclusively. The Q-law math (GVE, effectivity, periapsis barrier) is used internally by the manual strategy’s steering phases. The sequential element steering functions serve as internal helpers for the manual strategy’s plane_change and adjust_orbit phases.

Proximity Quotient

The proximity quotient measures distance from current to target orbit in orbital element space:

Q = Σ W_j * ((oe_j - oe_target_j) / S_j)²

Over orbital elements {a, e, i, Ω, ω} where:

W_j — weight for element j
S_j — scaling factor: max(|oe_j - oe_target_j|, floor) where floor normalizes contributions
Angular elements (i, Ω, ω) use shortest-arc wrapping

Scaling Floors

Element	Floor
a	0.01 × a
e	0.01
i, Ω, ω	0.01 rad

Default Weights

Element	Weight	Rationale
a	1.0	Semi-major axis matching
e	1.0	Eccentricity matching
i	1.0	Inclination matching
Ω	1.0	RAAN matching
ω	0.0	Ignored for near-circular orbits (game default)

Gauss Variational Equations

GVE matrix gives d(oe_j)/dt per unit acceleration in radial (R), transverse (T), normal (N) frame.

With p = a(1-e²), h = √(μp), r = p/(1 + e·cos ν), θ = ω + ν:

da/dt = (2a²/h)[e·sin(ν)·f_r + (p/r)·f_θ]
de/dt = (1/h)[p·sin(ν)·f_r + ((p+r)·cos(ν) + r·e)·f_θ]
di/dt = (r·cos(θ)/h)·f_h
dΩ/dt = (r·sin(θ)/(h·sin(i)))·f_h
dω/dt = (1/(h·e))[-p·cos(ν)·f_r + (p+r)·sin(ν)·f_θ] - (r·sin(θ)·cos(i))/(h·sin(i))·f_h

Singularity Handling

e < 1e-4: clamp in denominators, set ω = 0 (use argument of latitude)
sin(i) < 1e-4: clamp in denominators, set Ω = 0

Optimal Thrust Direction

The Q-gradient in RTN frame: D_k = Σ_j(∂Q/∂oe_j × G_jk) for k ∈ {R, T, N}.

Optimal thrust: f_hat = -normalize(D) (steepest descent on Q).

Periapsis Barrier

Prevents unsafe low-periapsis orbits:

P(rp) = k / (rp - min_rp)²    when rp < rp_safe

Parameters scale with the reference body:

safe_alt = SAFE_ORBIT_ALTITUDES[ref_body]        # per-body minimum altitude (m)
min_rp = body_radius + safe_alt
margin = max(PERIAPSIS_MARGIN_FLOOR, PERIAPSIS_MARGIN_FRACTION × body_radius)
rp_safe = min_rp + margin
k = PERIAPSIS_K_REF × (PERIAPSIS_REF_MARGIN / margin)²

The barrier strength k is scaled so the “repulsion feel” is consistent across bodies. A wider margin (large body) gets a softer barrier; a narrow margin (small body like Phobos) gets a stronger barrier.

Body	radius	safe_alt	min_rp	margin	rp_safe	k
Earth	6,371 km	200 km	6,571 km	95 km	6,666 km	110
Phobos	11 km	2 km	13 km	5 km	18 km	40,000
Jupiter	69,911 km	500 km	70,411 km	1,049 km	71,460 km	9

Gradients added to ∂Q/∂a and ∂Q/∂e:

∂P/∂a = ∂P/∂rp × (1 - e)
∂P/∂e = ∂P/∂rp × (-a)

GVE-Normalized Q-Gradient

In-plane GVE coefficients (a, e) are ~O(a²/h) while out-of-plane (i, Ω) are ~O(r/(h·sin(i))), differing by ~10⁷ in LEO. Without correction, tiny Δa residuals dominate the Q-gradient direction over large ΔΩ errors, causing the ship to chase the rotating transverse direction instead of thrusting along the stable orbit normal.

Per-element normalization: Each element’s contribution to the Q-gradient is divided by the orbit-max GVE magnitude for that element:

D_k = Σ_j (∂Q/∂oe_j × G_jk / max_ν ||G_j(ν)||)

This makes the gradient direction depend on Q-function partials (which properly reflect error magnitudes and scaling), not on GVE asymmetry. Orbit-max GVE norms are computed by sampling at QLAW_EFFECTIVITY_SAMPLES true anomalies — same samples used for effectivity.

Effectivity (Orbit-Max Normalized)

Effectivity is normalized to [0, 1] by dividing the GVE-normalized |D| by the orbit-max:

eff_normalized = |D_norm(ν_current)| / max_ν |D_norm(ν)|

With per-element normalization, effectivity correctly reflects how favorable the current position is for the dominant correction (e.g., peaks at nodes for RAAN maneuvers), rather than being dominated by elements with large GVE coefficients.

Effectivity is included in Q-law status texts so the user can monitor thrust effectiveness in real time. Status texts also show ΔRAAN alongside Δa, Δe, Δi for complete orbital element error visibility.

Proportional Thrust

Instead of binary burn/coast, Q-law always thrusts in the optimal direction with thrust level scaled by normalized effectivity:

thrust_level = min(1.0, eff_normalized / eff_threshold)

The effectivity threshold and attitude gate are per-ship-class to accommodate different inertia/torque characteristics:

Ship Class	eff_threshold	eff_gate	Rationale
fast_frigate	0.4	0.2	Higher thresholds — agile ship can afford to coast during low-effectivity arcs
cargo_hauler	0.2	0.1	Lower thresholds — slow ship needs all thrust it can get, direction is more stable (sequential strategy)

Per-class parameters are read from SHIP_CLASS_QLAW config and stored in the maneuver Redis hash at maneuver start.

At full effectivity (>= threshold): 100% thrust. At low effectivity: proportionally reduced thrust. This ensures continuous progress even when orbital geometry is unfavorable (e.g. large plane changes from near-circular orbits), at the cost of some fuel efficiency.

RTN to ICRF Conversion

RTN basis from position and velocity vectors:

R = r̂ (radial)
H = (r × v) / r × v (orbit normal)
T = H × R (transverse, in-track direction)

Phase Machine

ORBIT_MATCH ──(converged + dist>phase_entry_dist)──→ PHASE ──(dist<approach_entry_dist)──→ APPROACH ──(d<1km, Δv<1m/s)──→ COMPLETE
     │                                                                                        ↑
     └──────────(distance < approach_entry_dist at any time)──────────────────────────────────┘

Phase Transition Distances

Transition distances scale with orbit circumference so phase matching starts at an appropriate range for any body/orbit combination. Computed at maneuver start and stored in the maneuver Redis hash:

phase_entry_dist = max(PHASE_ENTRY_FLOOR, PHASE_ENTRY_FRACTION × 2π × target_a)
approach_entry_dist = max(APPROACH_ENTRY_FLOOR, APPROACH_ENTRY_FRACTION × 2π × target_a)

Body	target_a	phase_entry_dist	approach_entry_dist
Phobos (3km alt)	~14 km	~10 km (floor)	~1 km (floor)
Earth LEO (400km)	~6,771 km	~213 km	~43 km
Earth GEO	~42,164 km	~1,324 km	~265 km

Completion criteria remain absolute: distance < 1 km and relative velocity < 1 m/s.

ORBIT_MATCH Phase

Two steering modes selected automatically based on which orbital element errors dominate:

Mode Selection

Compute the out-of-plane dominance ratio: (Q_i + Q_raan) / Q_total. Mode transitions use hysteresis to prevent oscillation at the boundary:

Enter orbit-normal: oop_ratio > QLAW_RAAN_DOMINANCE_ENTER (0.8) and ΔΩ > QLAW_ORBIT_MATCH_RAAN_TOL (2°)
Exit orbit-normal: oop_ratio < QLAW_RAAN_DOMINANCE_EXIT (0.6) or ΔΩ ≤ QLAW_ORBIT_MATCH_RAAN_TOL

The current steering mode ("orbit_normal" or "gradient") is tracked in the maneuver Redis hash as steering_mode.

Equatorial guard: Orbit-normal steering is disabled when i < QLAW_MIN_INCLINATION_FOR_HNORM (5°), because RAAN is poorly defined for near-equatorial orbits (sin(i) → 0 makes the node vector degenerate).

This prevents H-norm activation when RAAN is already matched but inclination remains — H-thrust causes symmetric per-orbit Δi oscillation with zero net change in circular orbits, so gradient mode must handle residual inclination.

Orbit-Normal Steering (RAAN-dominated)

When out-of-plane errors (Δi + ΔΩ) dominate the Q-function (>80%), the Q-law gradient direction wanders at anti-nodes because tiny in-plane residuals dominate the normalized gradient. Instead of tracking this unstable direction, lock attitude to ±orbit normal:

Compute orbit normal: H = normalize(r × v) in ICRF — fixed in inertial space
Compute H-component of GVE-normalized Q-gradient: D_h = Σ_j (∂Q/∂oe_j × G_jh / norm_j)
- Naturally blends RAAN and inclination corrections: dΩ/dt ∝ sin(θ)·f_h, di/dt ∝ cos(θ)·f_h
- Sign follows steepest descent: h_sign = -sign(D_h)
- Direction flips 180° between orbit positions where the dominant correction changes sign
Out-of-plane effectivity: eff_h = |D_h| / max_ν |D_h|
- Peaks where out-of-plane thrust is most effective
- Near zero where thrust would not change i or Ω
Thrust level: eff_h × cos(alignment_angle)
Attitude is commanded every tick to ±H — approximately constant for half-orbit, then a single 180° flip

Why this works: H is inertially fixed, so the ship holds steady attitude for ~50 min (half orbit), executes one planned 180° flip (~60s), then holds again. No wheel saturation, no direction chasing.

Deactivation: Returns to gradient steering when out-of-plane dominance drops below threshold (in-plane elements become significant).

Gradient Steering (default)

Continuous burn with cosine-scaled thrust:

Compute GVE, effectivity, optimal thrust direction in RTN, convert to ICRF
Effectivity gate on attitude: Only update attitude target when eff >= QLAW_ATTITUDE_EFF_GATE (0.15). Below this threshold the Q-law direction is unreliable (dominated by tiny residuals from non-target elements that rotate with orbital motion). Holding attitude at low effectivity lets reaction wheels desaturate and prevents angular momentum buildup from chasing wandering directions.
Spin-damping slew limiter: When alignment_angle > QLAW_SPIN_DAMP_ALIGNMENT (45°), the attitude target is set to an intermediate direction at most QLAW_SPIN_DAMP_SLEW_RATE (20°) ahead of the ship’s current forward toward the Q-law optimal direction (f_icrf). This avoids commanding a 90°+ slew that would saturate the PD controller’s torque limit and trigger the spin-trap (kp_scale → 0 at moderate angular velocity), which causes oscillatory overshoot on cargo haulers. The bounded intermediate target generates enough torque (~22 kN·m at the saturation edge) to achieve ~5°/s slew rate while preventing overshoot. When alignment is within 45°, the target is set to f_icrf directly (normal tracking).
Thrust scaled by alignment: thrust *= max(0, cos(alignment_angle))
- Aligned (0°): full thrust
- 45° off: ~71% thrust
- 90°+: zero thrust (cosine negative, clamped)
Combined with effectivity scaling: actual_thrust = min(1.0, eff / threshold) * cos_factor

SOI guard: If apoapsis > SOI radius → steer retrograde (same as Lambert strategy).

PHASE Phase

Close along-track distance after orbit shape matches:

Compute phase angle from true longitude difference
Target ahead → retrograde (lower orbit → faster angular rate → catch up)
Target behind → prograde (raise orbit → slower → let target catch up)
Proportional throttle: max(0.01, min(0.2, |Δφ| / 90°)) — caps at 20%, gentle scaling
Alignment gate: No thrust until ship is aligned within 10° of the desired direction. Prevents thrust in the wrong direction during attitude transitions (e.g., retrograde → prograde).
SMA headroom limiter: Throttle is linearly reduced from full at Δa=0 to zero at Δa=a_tol, preventing phasing burns from exceeding the orbit match tolerance. da_factor = max(0, 1 − |Δa| / a_tol). Only limits thrust in the direction that increases divergence (prograde when Δa>0, retrograde when Δa<0). This self-regulating feedback keeps orbit shape within tolerance while allowing phasing to make gradual progress over many orbits.
SOI guard: If apoapsis > SOI radius → steer retrograde at full thrust (same as orbit_match phase)
SOI cap: If target behind (normally prograde) but apoapsis > 90% SOI radius → force retrograde instead (lower orbit → faster → gain phase the long way around). Prevents oscillation between prograde phasing and SOI guard at the SOI boundary.
Orbit divergence fallback: Each tick, recompute ship and target Keplerian elements. Check in-plane elements only (a, e) against 2× orbit_match tolerances — phasing burns tangentially and cannot fix inclination or RAAN, so out-of-plane elements are excluded from the divergence check. Separately, catastrophic OOP drift (i or Ω exceeding 5× tolerance) triggers fallback to plane_change. The 2× in-plane factor catches orbit shape drift early before it diverges far from what adjust_orbit established.
Status text: All phasing status messages include orbital element errors (Δa, Δe, Δi, ΔΩ) alongside phase angle and distance, matching orbit_match visibility.

APPROACH Phase

Two-mode approach: closing to reduce distance, then braking to kill relative velocity.

Closing Mode

When relative velocity is low but distance remains large (Δv < 20.0 m/s AND distance > 1000m), the ship steers toward the target to build closing velocity along a deceleration profile:

target_closing = min(50.0, sqrt(2 × 0.1 × distance))   # deceleration budget 0.1 m/s², cap 50 m/s
closing_vel = dot(d_hat, rel_v)                           # positive = approaching

If closing_vel < target_closing × 0.8: steer toward target (ATTITUDE_TARGET, proto enum 10) and thrust proportionally (vel_deficit / dv_per_tick)
If closing_vel >= target_closing × 0.8: switch to braking mode (existing TARGET_RETROGRADE)
The 0.8× hysteresis prevents oscillation at the mode boundary

The deceleration profile produces a natural brake-to-zero curve:

50 m/s at 250km → 14 m/s at 1km → 4.5 m/s at 100m → 1.4 m/s at 10m

Constants:

APPROACH_DECEL_BUDGET = 0.1 m/s² — deceleration budget for closing profile
MAX_CLOSING_VEL = 50.0 m/s — cap closing velocity
CLOSING_VEL_THRESHOLD = 20.0 m/s — switch to closing mode when Δv below this

Braking Mode

When relative velocity is high (≥ 20.0 m/s) or distance is small (≤ 1000m), kill relative velocity using TARGET_RETROGRADE attitude mode (original behavior).

Proportional Throttle

Thrust is scaled proportionally to prevent overshoot when the ship’s per-tick Δv exceeds the remaining relative velocity:

accel = max_thrust / total_mass          # ship acceleration (m/s²)
dv_per_tick = accel × time_scale         # Δv delivered per tick at full thrust
thrust = min(1.0, rel_vel / dv_per_tick)

No hardcoded floor — thrust goes as low as needed for sub-m/s corrections. This prevents oscillation when a high-TWR ship (e.g., fast frigate: 600 kN, ~20t → 30 m/s² accel) tries to cancel only 10-20 m/s of relative velocity.

Completion Criteria

Distance < 1 km AND relative velocity < 1.0 m/s.

Sequential Strategy

The "sequential" strategy corrects one orbital element at a time instead of blending all corrections into a single thrust direction via the Q-gradient. Each single-element thrust direction is a GVE row — stable and physically intuitive. This avoids the direction oscillation (~2°/s) that the blended Q-gradient produces, which high-inertia ships (cargo hauler: 6.4M kg·m² inertia, 22 kN·m max torque) cannot track.

When to Use

Use "sequential" for ships with high inertia-to-torque ratios where the PD controller cannot track the Q-law gradient’s rapid directional changes. Ships with low inertia or high torque authority can use "qlaw" (gradient) for faster convergence.

Phase Machine

Sequential shares the same phase machine as Q-law: ORBIT_MATCH → PHASE → APPROACH → COMPLETE. Only the orbit_match steering differs.

Element Priority

Elements are corrected one at a time in priority order. Out-of-plane elements first (±H direction, inertially fixed), then in-plane:

Priority	Element	GVE direction	ICRF behavior	Effectivity peaks
1	RAAN (Ω)	±H (normal)	Inertially fixed	θ = 90°, 270°
2	Inclination (i)	±H (normal)	Inertially fixed	θ = 0°, 180° (nodes)
3	SMA (a)	±T (transverse)	Rotates at orbital rate (~0.07°/s)	~Constant (circular)
4	Eccentricity (e)	R+T mix	Varies with ν	ν-dependent

Element Selection

Each tick, scan the priority list and pick the first element with |error| > tolerance. Re-evaluated each tick for robustness (handles perturbations, SOI guard disruptions that alter multiple elements).

Hysteresis: When the active element transitions from a higher-priority element (e.g., RAAN) to a lower-priority one (e.g., a+e), re-activating the higher-priority element requires exceeding a re-entry threshold of 2.0 × tolerance instead of the base tolerance. This prevents chattering at tolerance boundaries, which causes rapid attitude target switching (~90° between orbit-normal and transverse) that high-inertia ships cannot track. The 2.0× factor (increased from 1.5×) provides sufficient margin for cross-coupling perturbations — e.g., SMA burns perturb inclination by ~0.5° per orbit on near-circular orbits, which at 1.5× repeatedly crossed the re-entry threshold for a 0.5° inclination tolerance.

The prev_element is stored in the maneuver hash and passed to sequential_element_steering(). When checking elements with higher priority than prev_element, the re-entry factor (2.0×) is applied. When prev_element has itself converged (below base tolerance), normal tolerances apply for all elements.

Combined Out-of-Plane Steering

When both RAAN and inclination exceed their effective tolerances simultaneously, they are corrected together as a single “oop” (out-of-plane) element using a tolerance-normalized blend, exactly analogous to the a+e in-plane blend:

w_raan = (Ω_ship - Ω_target) / tol_raan
w_i    = (i_ship - i_target) / tol_i

D_k = w_raan × GVE_raan_hat[k] + w_i × GVE_i_hat[k]
f_rtn = -normalize(D)

Why combined: Correcting RAAN and i separately causes oscillation — each correction cross-couples to the other through orbit-normal thrust, and the priority system bounces between them. During the oscillation, in-plane elements (eccentricity) grow unchecked because the ship never reaches a+e priority. Combined steering corrects both simultaneously: RAAN-dominant orbital positions get RAAN correction, i-dominant positions get i correction, within a single orbit pass.

When only one of RAAN or i exceeds tolerance, single-element steering is used instead.

Periapsis Safety Override

When the current periapsis radius rp = a × (1 - e) drops below rp_safe, all out-of-plane corrections (combined “oop”, single “raan”, single “i”) are skipped entirely, and the element scanner falls through to in-plane a+e correction.

Why this is necessary: Orbit-normal thrust (±H) has zero first-order effect on SMA or eccentricity. When the orbit is degraded to a dangerously low perigee, OOP corrections cannot raise it — only in-plane (transverse/radial) thrust can increase SMA and circularize the orbit. Without this override, the priority system keeps the ship on OOP while perigee drops toward atmospheric entry.

rp_current = a × (1 - e)
periapsis_critical = rp_safe > 0 and rp_current < rp_safe

if periapsis_critical:
    skip combined OOP block
    skip single-element RAAN and i in priority scan
    → falls through to a+e correction

The override is transient: once in-plane thrust raises perigee above rp_safe, normal OOP priority resumes. No hysteresis is applied to the safety override — it activates immediately when perigee is critical.

Single-Element Steering

For active element j:

f_rtn = -sign(error_j) × normalize(GVE_row_j(ν))
eff_j = ||GVE_row_j(ν_current)|| / max_ν ||GVE_row_j(ν)||

The sign is determined by the element error direction (e.g., RAAN too high → thrust in -H direction to reduce it). The GVE row gives the physically correct RTN direction for changing that element.

Combined a+e Steering

When the active element is SMA (“a”) or eccentricity (“e”), thrust direction uses a tolerance-normalized blend that simultaneously corrects both coupled in-plane elements. Transverse thrust for SMA also changes eccentricity (GVE coupling), and correcting them independently creates costly correction cycles.

For active element ∈ {a, e}:

w_a = (a_ship - a_target) / tol_a     # signed, tolerance-normalized
w_e = (e_ship - e_target) / tol_e     # signed, tolerance-normalized

GVE_a_hat = normalize(GVE_a_row)       # unit direction for SMA change
GVE_e_hat = normalize(GVE_e_row)       # unit direction for eccentricity change

D_k = w_a × GVE_a_hat[k] + w_e × GVE_e_hat[k]   for k ∈ {r, θ, h}
f_rtn = -normalize(D)

Why tolerance-normalized, not Q-gradient: The Q-function’s inverse-square scaling (dQ/dx ∝ 1/|Δx|) creates pathological imbalance between elements with different units. With Δa = 5000 km and Δe = 0.04, the ratio dQ_de / dQ_da ≈ 130 million — the blend becomes pure eccentricity correction. Since eccentricity GVE flips sign twice per orbit, this produces retrograde thrust for half the orbit, canceling SMA progress. Tolerance normalization weights elements by how many tolerances they are from convergence: w_a/w_e = (5000km/10km)/(0.04/0.01) = 125:1, correctly prioritizing SMA.

Natural convergence behavior: When an element is within tolerance, its weight is < 1.0 and is dominated by the unconverged element’s weight (typically hundreds). No explicit suppression is needed — the tolerance normalization naturally reduces converged elements’ contribution to negligible levels.

Out-of-plane elements (RAAN, i) use combined “oop” steering when both exceed tolerance, or single-element GVE rows when only one needs correction
Effectivity = D_current / max_ν( D ) — blended orbit-max, sampled at 12 true anomaly points
Status text shows “a+e” as element name when blending

Effectivity

Single-element effectivity: ratio of current GVE row magnitude to orbit-max GVE row magnitude, sampled at 12 true anomaly points. Maps to [0, 1].

Thrust scaling same as Q-law: thrust = min(1.0, eff / QLAW_EFFECTIVITY_THRESHOLD) × cos(alignment_angle).

Attitude Handling

Uses the same slew limiter as Q-law gradient mode (_intermediate_direction at 20° max step) for initial alignment recovery. Once aligned, single-element directions are inherently stable — the ship stays aligned without chasing a wandering target.

Ship Acceleration

The tick-engine needs ship acceleration (m/s²) for dynamic effectivity gating. Ship class configs (max_thrust) are defined in the physics service. The tick-engine maintains a matching lookup table (SHIP_CLASS_THRUST) keyed by ship class name. Total mass = dry_mass + fuel (both available from ship Redis hash).

Status Text

Format: "Seq. {element} {thrust%}, eff {eff}"

Manual Strategy

The "manual" strategy emulates a human pilot’s committed-phase approach to rendezvous. Unlike sequential (which re-evaluates element priority every tick, causing thrashing), manual locks into discrete phases and completes each before advancing.

Motivation

Sequential’s per-tick priority re-evaluation causes thrashing between OOP and in-plane corrections. A human pilot does rendezvous in committed phases: plane change first, then orbit shaping (Hohmann-like Ap/Pe adjustments), then phase match, then approach.

Phase Machine

PLANE_CHANGE ──→ ADJUST_ORBIT ──→ PHASE ──→ APPROACH ──→ COMPLETE
     ^                │  ^           │
     └────────────────┘  └───────────┘
     (OOP diverged >5× tol)  (Δe diverged >5× tol)

Phases already converged at start are skipped.

Apoapsis/Periapsis GVE Derivation

The key insight is that apoapsis and periapsis can be controlled independently via decomposed GVE rows:

Ap = a × (1 + e)  →  d(Ap)/dt = (1+e)·d(a)/dt + a·d(e)/dt
Pe = a × (1 - e)  →  d(Pe)/dt = (1-e)·d(a)/dt - a·d(e)/dt

As GVE row vectors in RTN:

GVE_Ap = (1+e) × GVE_a + a × GVE_e
GVE_Pe = (1-e) × GVE_a - a × GVE_e

Natural decoupling:

At periapsis (ν = 0): GVE_Ap peaks, GVE_Pe ≈ 0 → burn here to adjust Ap
At apoapsis (ν = π): GVE_Pe peaks, GVE_Ap ≈ 0 → burn here to adjust Pe

Since Ap and Pe effectivity are complementary (one peaks when the other is near zero), the ship can work on both in the same phase by picking whichever element has higher effectivity at the current orbital position. This gives nearly continuous productive burning with minimal cross-coupling.

Effectivity-Weighted Thrust

Both plane_change and adjust_orbit use the same thrust pattern: thrust_pct=1.0 with effectivity weighting.

actual_thrust = 1.0 × eff × cos(alignment_angle)

The steering function computes max_gve (peak GVE over orbit) for effectivity normalization: eff = |GVE(ν)| / max_gve.

Effectivity weighting (× eff): naturally concentrates thrust near apsides where GVE magnitude peaks (eccentric orbits) while allowing continuous thrust for near-circular orbits where eff is ~uniform. No dead zones from binary gating.
No proportional throttle: previous versions used throttle = min(1.0, |error| / correction_per_orbit) which over-throttled for small errors (e.g., 0.12% throttle for a 1km error on a near-circular orbit). Effectivity weighting alone provides sufficient modulation — it concentrates burns where they are most effective, matching the plane_change pattern which converges reliably without proportional throttle.

Attitude always updates (direction is prograde/retrograde, always well-defined). Thrust is continuous, not gated.

Phase Descriptions

PLANE_CHANGE:

Combined RAAN+i steering (tolerance-normalized raw GVE blend, same as sequential “oop”)
Single-element fallback when only one OOP element exceeds tolerance
Effectivity-weighted thrust: actual_thrust = 1.0 × eff × cos(alignment_angle). Attitude always updates (orbit-normal direction is always well-defined). No binary eff_gate.
Periapsis safety: if rp < rp_safe, temporarily switch to in-plane a+e steering
Transition: when both Δi ≤ tol_i and ΔΩ ≤ tol_raan

ADJUST_ORBIT:

Each tick: compute eff_ap, eff_pe at current position. Pick element with higher effectivity (both need correction) or the one that needs correction.
Attitude always updates toward the selected element’s GVE direction (prograde/retrograde)
Effectivity-weighted thrust: actual_thrust = 1.0 × eff × cos(alignment_angle). Same pattern as plane_change — continuous modulation, no binary gate. Thrust naturally concentrates near apsides for eccentric orbits (high eff) and distributes broadly for near-circular orbits (uniform eff).
Convergence: both ΔAp ≤ a_tol and ΔPe ≤ a_tol
Burn ETA: Estimated time from |error| / peak_rate where peak_rate = max_gve × accel. Uses peak (not mean) rate since effectivity weighting already modulates thrust. Displayed as minutes or hours.
Catastrophic OOP fallback: regress to PLANE_CHANGE if Δi > 5×tol_i or ΔΩ > 5×tol_raan

PHASE:

Existing phase-matching logic (pro/retrograde based on phase angle)
Transition to approach: uses approach_entry_dist (not phase_entry_dist) — phasing continues until the ship is within approach range (~75km for LEO), giving it time to close the along-track gap naturally
SMA headroom limiter: da_factor = max(0, 1 − |Δa| / (PHASE_SMA_HEADROOM_FACTOR × a_tol)) — linearly ramps throttle from full at Δa=0 to zero at Δa=PHASE_SMA_HEADROOM_FACTOR × a_tol. Allows phasing to use a significantly different orbit for faster phase closure while still limiting divergence. Self-regulating: orbit raises gradually, phasing progresses over fewer orbits with larger angular rate difference.
ETA: Estimated time remaining from orbital mean motion difference: eta = |Δφ| / |n_target − n_ship| where n = √(μ/a³). Updated each tick. Displayed in status_text as “Xm” (< 1 hour) or “X.Yh” (≥ 1 hour).

Catastrophic eccentricity check: if Δe > MANUAL_CATASTROPHIC_FALLBACK_FACTOR × e_tol, fall back to ADJUST_ORBIT. The SMA headroom limiter controls Δa but not Δe — sustained prograde/retrograde phasing burns can increase eccentricity over many orbits, creating a mismatched orbit shape that makes phasing ineffective. This check detects the divergence and re-enters ADJUST_ORBIT to correct Ap/Pe before resuming phasing.

Catastrophic OOP check: if Δi > 5×tol_i or ΔΩ > 5×tol_raan, fall back to PLANE_CHANGE.

APPROACH:

Existing approach logic (TARGET_RETROGRADE)

Constants

Constant	Value	Description
`MANUAL_CATASTROPHIC_FALLBACK_FACTOR`	5.0	Regress to earlier phase if OOP elements exceed this × tolerance
GVE orbit samples	360	Samples around orbit (1° resolution) for mean_gve computation

Constants

Constant	Value	Description
`QLAW_ORBIT_MATCH_A_TOL_RELATIVE`	0.001	SMA tolerance as fraction of target SMA (0.1%)
`QLAW_ORBIT_MATCH_A_TOL_FLOOR`	100,000 m	Minimum SMA tolerance floor (100 km)
`QLAW_ORBIT_MATCH_E_TOL`	0.01	Eccentricity tolerance
`QLAW_ORBIT_MATCH_I_TOL`	0.5°	Inclination tolerance
`QLAW_ORBIT_MATCH_RAAN_TOL`	2.0°	RAAN tolerance (weighted by `sin(max(i_ship, i_target))` — RAAN is undefined for equatorial orbits)
`QLAW_EFFECTIVITY_THRESHOLD`	0.3	Default full-thrust effectivity (per-class override via `SHIP_CLASS_QLAW`)
`QLAW_ATTITUDE_EFF_GATE`	0.15	Default attitude gate (per-class override via `SHIP_CLASS_QLAW`)
`QLAW_SPIN_DAMP_ALIGNMENT`	45.0°	Freeze attitude target above this alignment angle (let controller damp spin)
`QLAW_SPIN_DAMP_SLEW_RATE`	20°	Max angle ahead of ship forward for intermediate attitude target (keeps PD controller error bounded)
`QLAW_EFFECTIVITY_SAMPLES`	12	True anomaly samples for orbit-max normalization
`QLAW_RAAN_DOMINANCE_ENTER`	0.8	Enter orbit-normal steering when oop_ratio exceeds this
`QLAW_RAAN_DOMINANCE_EXIT`	0.6	Exit orbit-normal steering when oop_ratio drops below this
`QLAW_MIN_INCLINATION_FOR_HNORM`	5°	Minimum inclination for orbit-normal steering (equatorial guard)
`PHASE_SMA_HEADROOM_FACTOR`	5.0	Allow 5× a_tol SMA deviation during phasing for faster phase closure
`PHASE_ENTRY_FRACTION`	0.005	Phase transition distance as fraction of orbit circumference (0.5%)
`PHASE_ENTRY_FLOOR`	10,000 m	Minimum phase transition distance (10 km)
`APPROACH_ENTRY_FRACTION`	0.001	Approach transition distance as fraction of orbit circumference (0.1%)
`APPROACH_ENTRY_FLOOR`	1,000 m	Minimum approach transition distance (1 km)
`PERIAPSIS_MARGIN_FRACTION`	0.015	Barrier margin as fraction of body radius (1.5%)
`PERIAPSIS_MARGIN_FLOOR`	5,000 m	Minimum barrier margin (5 km)
`PERIAPSIS_K_REF`	100.0	Reference barrier strength
`PERIAPSIS_REF_MARGIN`	100,000 m	Reference margin for barrier scaling

Effective SMA tolerance = max(floor, relative × target_a). The 100 km floor ensures near-circular orbits don’t require impractically precise Ap/Pe corrections (which have strong cross-coupling at low eccentricity). For larger orbits (a > 100,000 km), the 0.1% relative tolerance takes over. Since the phasing phase follows ADJUST_ORBIT, residual Ap/Pe errors within tolerance are acceptable — the approach phase handles final velocity matching.

Burn Approach Alert

During automated maneuvers, the ship spends long periods coasting between burn windows (e.g., waiting to reach apoapsis/periapsis in adjust_orbit, or waiting to reach orbital nodes in plane_change). An escalating audio alert sounds 60 seconds before each burn, letting players monitor maneuvers without watching the screen.

Backend: `burn_eta_seconds` Field

The burn_eta_seconds field is set in the maneuver Redis hash during coasting periods:

adjust_orbit: Computes dt (time to next Pe/Ap), then subtracts the burn lead time — the time from when the effectivity gate opens to when the ship reaches the apsis. The burn starts when eff ≥ dynamic_gate, which corresponds to a true anomaly offset of arccos(√gate) from the apsis. Lead time = T_orbit × arccos(√gate) / (2π). Set burn_eta_seconds = int(dt - lead_time) (clamped ≥ 0). Always computed (not gated on thrust level), so the countdown reaches the ≤60s alert window before burning begins.

plane_change: Computes time to next orbital node (where out-of-plane effectivity peaks). Argument of latitude u = ω + ν. Nodes at u = 0 and u = π. Find nearest node’s true anomaly, convert via eccentric anomaly → mean anomaly → time. Always computed (not gated on thrust level), since plane_change uses effectivity-weighted thrust that ramps up gradually near nodes — burn_eta_seconds reflects time-to-node regardless of current thrust state.

phase: Set burn_eta_seconds = int(eta_s) where eta_s is the estimated time for the phase angle to close (computed from mean motion difference). This triggers the burn approach alert as the ship nears the approach transition. Set to 0 when eta_s is not computable (e.g., SMA difference too small).

approach: Do not set the field (continuous terminal guidance, no discrete burn windows). When transitioning from phase to approach, delete the stale burn_eta_seconds field via remove_maneuver_fields.

Client: Escalating Beep Schedule

When burn_eta_seconds is present, ≤ 60, and > 0, the client plays escalating beeps:

ETA Range	Beep Interval	Urgency
60–30s	Every 10s	0.0–0.5 (low pitch, long tone)
30–10s	Every 5s	0.5–0.8 (medium pitch)
10–0s	Every 2s	0.8–1.0 (high pitch, short tone)

Audio characteristics:

Sine oscillator: frequency = 440 + urgency × 440 Hz (A4 → A5)
Duration: 200 − urgency × 100 ms (shorter = more urgent)
Gain: 0.2 (quieter than targeted alert)

The client maintains a local 1-second countdown timer synced to the server’s burn_eta_seconds value (resync each maneuver_status update, ~5s interval). This provides smooth countdown without requiring per-second server messages.

Settings

burnAlertEnabled (default: true) — toggle in Audio settings tab. Disabling stops all burn approach beeps immediately.

High-Energy Transfer Planning

Pre-computed two-burn transfer orbits for fast/express/parabolic strategies. Used by the adjust_orbit phase when the strategy is not Hohmann.

Transfer Orbit Geometry

Given ship at circular orbit radius r1 and target at circular orbit radius r2:

Transfer apogee: r_a = factor × max(r1, r2) where factor depends on strategy
Transfer SMA: a_t = (r1 + r_a) / 2
Transfer eccentricity: e_t = (r_a - r1) / (r_a + r1)

Strategy	Factor	Notes
fast	2.0	Moderate energy transfer
express	5.0	High energy transfer
parabolic	∞	e_t = 1, departure at escape velocity

Departure Δv (prograde, at periapsis of transfer)

v_transfer = sqrt(mu × (2/r1 - 1/a_t))
v_circular = sqrt(mu / r1)
dv1 = v_transfer - v_circular

For parabolic: dv1 = sqrt(2 × mu / r1) - sqrt(mu / r1)

Arrival Δv (at target orbit crossing)

At radius r2, the transfer orbit velocity and flight path angle:

h = r1 × v_transfer           # angular momentum (periapsis)
v_t = sqrt(mu × (2/r2 - 1/a_t))  # speed at r2
v_tangential = h / r2          # tangential component
v_radial = sqrt(v_t² - v_tangential²)  # radial component
gamma = atan2(v_radial, v_tangential)   # flight path angle

v_target = sqrt(mu / r2)      # target circular velocity
dv2 = sqrt(v_t² + v_target² - 2 × v_t × v_target × cos(gamma))

For parabolic: v_t = sqrt(2 × mu / r2) (vis-viva for e=1, a→∞).

Transfer Time

Elliptic (fast, express): Kepler’s equation from periapsis to target orbit crossing.

cos_nu = (a_t × (1 - e_t²) / r2 - 1) / e_t    # true anomaly at r2
E = 2 × atan2(sqrt(1 - e_t) × sin(nu/2), sqrt(1 + e_t) × cos(nu/2))
M = E - e_t × sin(E)
t = M / sqrt(mu / a_t³)

Parabolic: Barker’s equation.

D = tan(nu / 2)
t = sqrt(2 × r1³ / mu) × (D + D³/3) / 2

Where nu is the true anomaly at r2 for the parabolic orbit: cos_nu = (p/r2 - 1) with p = 2 × r1.

Estimate Function

transfer_estimate(r_ship, r_target, mu, strategy) computes {total_dv, transfer_time} for UI display. Covers all four strategies including Hohmann:

Hohmann estimate:

dv = |sqrt(mu/r1) × (sqrt(2r2/(r1+r2)) - 1)| + |sqrt(mu/r2) × (1 - sqrt(2r1/(r1+r2)))|
t = π × sqrt(((r1+r2)/2)³ / mu)

Lambert Solver

Solves Lambert’s problem: given two position vectors and a time of flight, find the velocity vectors for the transfer orbit. Uses universal variable formulation with Stumpff functions.

Algorithm

Compute transfer angle from position vectors. Direction (prograde/retrograde) from cross product sign.
Compute parameter A from geometry.
Newton-Raphson iteration on universal variable ψ with bisection safeguard:
- ψ > 0: elliptic transfer
- ψ = 0: parabolic
- ψ < 0: hyperbolic
Compute f, g, gdot Lagrange coefficients from converged ψ.
Return v1, v2 velocity vectors.

Functions

lambert_solve(r1_vec, r2_vec, tof, mu, prograde) → (v1, v2) or None
lambert_transfer(r1_vec, v1_current, r2_vec, tof, mu) → {dv1, dv2, dv1_vec, v1_transfer, v2_transfer, transfer_time} or None

Integration

Added as "lambert" strategy in STRATEGY_LABELS. Can be used by transfer_plan when strategy == "lambert" to compute transfer from actual position vectors instead of radii-based Hohmann. Lambert needs the target’s predicted future position — propagate target mean anomaly forward by estimated TOF.

Optimal Departure Timing

Scans candidate departure positions around the ship’s orbit to find the minimum total Δv point (Oberth effect optimization).

Algorithm

Sample 12 true anomaly positions over 2 orbits
At each candidate: compute ship radius, Hohmann Δv from that radius via compute_full_transfer
Return best candidate (minimum total Δv)
Only activate for eccentric orbits (e > 0.01) — circular orbits have uniform departure cost
Require > 5% savings over immediate departure

Function

optimal_departure_scan(oe_ship, oe_target, mu, strategy, n_samples) → {total_dv, nu_departure, wait_time, r_departure} or None

Bi-Elliptic Transfers

For large radius ratios (> 11.94:1), a three-burn bi-elliptic transfer can be more Δv-efficient than Hohmann.

Algorithm

Three burns:

At r1: raise apoapsis to r_intermediate
At r_intermediate: lower periapsis to r2
At r2: circularize

dv1 = |v_dep1 - v_circ1|    (r1 → transfer ellipse 1)
dv2 = |v_dep2 - v_arr1|     (transfer ellipse 1 → 2 at r_intermediate)
dv3 = |v_circ2 - v_arr2|    (transfer ellipse 2 → circular at r2)

Constants

Constant	Value	Description
`BI_ELLIPTIC_RATIO_THRESHOLD`	11.94	Minimum r_max/r_min ratio for bi-elliptic consideration
`BI_ELLIPTIC_FACTOR`	2.0	r_intermediate = factor × max(r1, r2)
`BI_ELLIPTIC_MIN_SAVINGS`	0.05	Minimum 5% Δv savings over Hohmann

Functions

bi_elliptic_dv(r1, r2, r_intermediate, mu) → {dv1, dv2, dv3, total_dv} or None
should_use_bi_elliptic(r1, r2, mu, min_savings) → (bool, r_intermediate)

Implementation

Bi-elliptic is implemented as two sequential Hohmann transfers via the existing state machine:

Cycle 1: r1 → r_intermediate (Hohmann)
Cycle 2: r_intermediate → r2 (Hohmann)

The transfer_cycle field tracks which cycle. _circularize_complete() routes cycle 1 completion to re-enter transfer_plan for cycle 2.

Keplerian Elements

Computed from Cartesian state (body-relative position and velocity):

Returns None for degenerate orbits (h ≈ 0, r ≈ 0, parabolic energy ≈ 0)
Accepts hyperbolic orbits (e ≥ 1.0, a < 0) — see below
Edge cases: e < 1e-4 → ω = 0; sin(i) < 1e-4 → Ω = 0
Elements: {a, e, i, raan, omega, nu, h_mag, r_mag, p}

Hyperbolic Orbit Support

keplerian_from_cartesian accepts hyperbolic trajectories (e ≥ 1.0, a < 0). This is needed because cross-SOI capture can leave ships on hyperbolic orbits around the target body. Previously, the function returned None for e ≥ 1.0 or a ≤ 0, which caused downstream phases (plane_change, phase, phase_coast) to enter unrecoverable “degenerate orbit, coasting” states.

Why the math works:

Semi-latus rectum: p = a(1 − e²) is positive for hyperbolic orbits because both a < 0 and (1 − e²) < 0, so their product is positive
Specific angular momentum: h = √(μp) is real since p > 0
GVE coefficients are valid: all GVE formulas use p, h, and trigonometric functions of ν — none require a > 0 or e < 1
True anomaly via atan2 works for all eccentricities

Derived quantities for hyperbolic orbits:

Periapsis: a(1 − e) — positive (both a < 0 and (1 − e) < 0)
Apoapsis: undefined (∞) — hyperbolic orbits do not close
ap_pe_converged returns False when e ≥ 1.0 (no meaningful apoapsis convergence)

Phase fallback list: Phases that can proceed without valid Keplerian elements (using fallback radial estimates) include: transfer_plan, transfer_burn, transfer_coast, circularize, departure_wait, capture, plane_change, phase, phase_coast.

Hyperbolic recovery in orbit-matching phases: If keplerian_from_cartesian returns elements with e ≥ 1.0 in plane_change, phase, or phase_coast, the phase transitions to transfer_plan for re-planning. This allows the ship to recover from unexpected hyperbolic capture orbits.

Post-Capture Routing

After cross-SOI capture (both from the capture phase and from already-bound SOI entry in the interplanetary phase), the next phase depends on the maneuver strategy:

Strategy	Next Phase	Rationale
`brachistochrone`	`brachistochrone`	High-thrust ships skip orbit-matching; direct point-to-point guidance is faster
All others (`hohmann`, `fast`, `express`, `parabolic`, `manual`, `sequential`, `qlaw`)	`circularize`	Low-thrust ships need orbit-matching sequence (circularize → phase → approach)

This applies at two transition points:

Capture phase exit (line ~497): when retrograde braking achieves bound orbit with apoapsis within SOI
Interplanetary already-bound SOI entry (line ~222): when ship enters target SOI already on a bound orbit fitting within SOI

Module Structure

services/tick-engine/src/qlaw.py is a pure-math module (no async, Redis, or gRPC). Functions are organized into logical sections with # ── separator comments and a table of contents at the top of the file.

Sections

Section	Functions	Description
Orbital Element Conversion	`_angle_diff`, `keplerian_from_cartesian`, `keplerian_from_state`	Cartesian ↔ Keplerian
GVE & Q-function	`gve_coefficients`, `proximity_quotient`, `out_of_plane_dominance`, `orbit_normal_steering`, `_q_gradient_partials`, `compute_gve_norms`, `_d_magnitude`, `optimal_thrust_direction_rtn`, `rtn_to_icrf`, `orbit_match_converged`, `compute_effectivity`	Gauss variational equations, Q-function evaluation, thrust direction
Sequential Element Steering	`sequential_element_steering`, `sequential_eta`, `gve_ap_pe`	Priority-ordered single-element steering with hysteresis
Manual Strategy Steering	`oop_steering`, `ap_pe_steering`, `ap_pe_combined_steering`, `ae_combined_steering`, `aeio_combined_steering`, `oop_interleave_check`, `oop_converged`, `ap_pe_converged`, `oop_coast_steering`	Fixed-phase steering (plane_change, adjust_orbit)
Phasing	`compute_phase_angle`, `compute_phasing_orbit_sma`	Phase angle computation and phasing orbit design. `compute_phasing_orbit_sma` accepts optional `mu` to iteratively compensate for Hohmann transfer phase shift (#1066)
Hohmann & Transfer Planning	`hohmann_transfer_time`, `hohmann_departure_dv`, `hohmann_burn_plan`, `transfer_burn_plan`, `compute_full_transfer`, `transfer_estimate`, `compute_departure_phase`	Transfer orbit mechanics
Anomaly Helpers	`angle_near_zero`, `time_to_anomaly`, `angle_near_pi`, `manual_eta`	True anomaly utilities and ETA estimation
Lambert Solver	`_stumpff_c2`, `_stumpff_c3`, `lambert_solve`, `lambert_transfer`, `optimal_departure_scan`	Lambert’s problem via universal variable
Bi-elliptic & Brachistochrone	`bi_elliptic_dv`, `should_use_bi_elliptic`, `brachistochrone_time_estimate`	Alternative transfer strategies
Kepler Propagation	`kepler_propagate`	Universal variable propagation with Stumpff functions

Q-law / Lyapunov Steering for Rendezvous

Overview

Proximity Quotient

Scaling Floors

Default Weights

Gauss Variational Equations

Singularity Handling

Optimal Thrust Direction

Periapsis Barrier

GVE-Normalized Q-Gradient

Effectivity (Orbit-Max Normalized)

Proportional Thrust

RTN to ICRF Conversion

Phase Machine

Phase Transition Distances

ORBIT_MATCH Phase

Mode Selection

Orbit-Normal Steering (RAAN-dominated)

Gradient Steering (default)

PHASE Phase

APPROACH Phase

Closing Mode

Braking Mode

Proportional Throttle

Completion Criteria

Sequential Strategy

When to Use

Phase Machine

Element Priority

Element Selection

Combined Out-of-Plane Steering

Periapsis Safety Override

Single-Element Steering

Combined a+e Steering

Effectivity

Attitude Handling

Ship Acceleration

Status Text

Manual Strategy

Motivation

Phase Machine

Apoapsis/Periapsis GVE Derivation

Effectivity-Weighted Thrust

Phase Descriptions

Constants

Constants

Burn Approach Alert

Backend: burn_eta_seconds Field

Client: Escalating Beep Schedule

Settings

High-Energy Transfer Planning

Transfer Orbit Geometry

Departure Δv (prograde, at periapsis of transfer)

Arrival Δv (at target orbit crossing)

Transfer Time

Estimate Function

Lambert Solver

Algorithm

Functions

Integration

Optimal Departure Timing

Algorithm

Function

Bi-Elliptic Transfers

Algorithm

Constants

Functions

Implementation

Keplerian Elements

Hyperbolic Orbit Support

Post-Capture Routing

Module Structure

Sections

Backend: `burn_eta_seconds` Field