DEV Community

Cover image for I Ported a Peer-Reviewed Astrodynamics Paper to TypeScript. Here's Why.
Sarkis M
Sarkis M

Posted on

I Ported a Peer-Reviewed Astrodynamics Paper to TypeScript. Here's Why.

Spacecraft proximity operations math is world-class. Papers in JGCD present elegant, validated formulations — analytical State Transition Matrices, density-model-free drag representations, quasi-nonsingular element sets that avoid coordinate singularities. Decades of theoretical refinement.

The tooling is plot3() in a MATLAB figure window. Scripts with magic numbers. Variable names like aa, bb, dd_lambda because MATLAB doesn't have types. Visualization is a non-interactive figure window. Collaboration is a shared drive. Every new team reimplements the same papers from scratch, introduces the same bugs, builds the same bad UIs.

The web has solved this problem for virtually every other domain. Finance has interactive dashboards. Data science has Observable notebooks. Mapping has Mapbox. Astrodynamics has MATLAB figure windows and STK if you can afford the license.

I decided to change that for one particular domain. I took Koenig, Guffanti, and D'Amico's 2017 paper on State Transition Matrices for perturbed relative motion — from the Journal of Guidance, Control, and Dynamics — and implemented it as a TypeScript library with a 3D interactive mission planner that runs entirely in the browser. Adam Koenig, the paper's lead author and a current colleague, has validated the implementation. This is a production-quality port of peer-reviewed math, not a weekend approximation.

Try the live demo →

No MATLAB license. No install. No backend. Click to place waypoints in 3D space, toggle perturbation models, watch a deputy spacecraft fly the computed trajectory. This is the same math that flies real missions.

What I Actually Built

The project has two layers:

1. An orbital mechanics library (src/orbital/) — pure TypeScript, zero framework dependencies. Implements three fidelity levels of State Transition Matrices from the Koenig et al. paper:

Model Matrix Size What It Captures
Keplerian 6×6 Two-body motion only
J2-Perturbed 6×6 Earth's equatorial bulge (nodal regression, apsidal precession)
J2 + Drag 7×7 / 9×9 Atmospheric drag with density-model-free formulation

Plus a targeting solver (Newton-Raphson iterative root-finding for rendezvous), a mission planner that chains multi-waypoint transfers, Gauss Variational Equations (how impulsive burns map to orbit changes), and ROE ↔ RIC coordinate transforms.

2. A 3D mission planner — React 19, React Three Fiber, Zustand, Tailwind. Interactive waypoint placement, real-time trajectory recomputation, simulation playback with telemetry HUD, and data export to JSON/CSV.

The library is framework-agnostic. You could use it in a Node script, a Bun project, a WASM module — the React app is just one consumer.

Why This Paper

Not all astrodynamics papers are created equal for the purpose of porting to a web environment.

Koenig et al.'s approach uses analytical State Transition Matrices instead of numerical integration. This is the key architectural decision that makes browser-native implementation viable.

Numerical integration (Runge-Kutta, etc.) requires hundreds or thousands of timesteps, each with expensive trigonometric evaluations. It's inherently iterative and hard to make fast in JavaScript. Analytical STMs reduce propagation to a single matrix-vector multiply — for a 6×6 matrix, that's 36 multiplies and 30 adds. Done. The state at any future time falls out directly.

This means the targeting solver — which calls propagation on every Newton-Raphson iteration, and again 6 times per Jacobian evaluation — runs in microseconds instead of milliseconds. Which means the mission planner can recompute trajectories on every frame while the user drags a waypoint in 3D space. Which means the tool feels interactive rather than "click compute and wait."

Drag a waypoint 200 meters radially outward. The planner recomputes departure burns, coast arcs, and arrival burns in under a millisecond. That's analytical propagation at work — not clever optimization, not caching, just math that's inherently fast enough.

The paper also uses a Relative Orbital Elements (ROE) formulation instead of Cartesian state vectors. ROE describes how two orbits differ — relative semi-major axis, relative eccentricity vector, relative inclination vector — rather than tracking two absolute positions. Perturbations separate cleanly in this space. J2 effects show up as known secular drift rates. Drag adds coupling terms with known structure. Each perturbation layer extends the STM while preserving the analytical form.

Porting PhD-Level Math to TypeScript: What That Actually Involves

This is the part that's hard to convey in a README.

The paper is 20+ pages of derivations. The STMs reference orbital factors (κ, η, P, Q, and others) that are themselves functions of eccentricity, inclination, and Earth's J2 coefficient. The drag model introduces a "density-model-free" formulation where you don't need an atmospheric density profile — instead, you parameterize the drag effect as rates of change of ROE components. The 9×9 arbitrary-eccentricity drag matrix has elements that are polynomials in these factors.

Translating this means:

Reading the paper deeply enough to implement it, not just cite it. You need to understand what each matrix element represents physically, because when your output is wrong, the debugger doesn't tell you "sign error in row 3, column 5 of the J2 STM." It tells you the trajectory curves the wrong way. Knowing that row 3 corresponds to δex (eccentricity x-component) and that a sign flip there would reverse the apsidal precession is how you find the bug.

Handling the notation translation. Academic papers use compact notation that maps poorly to code. Subscripts, superscripts, overlines, primes — all of these become variable names. The paper's Ω̇ becomes something like raanDrift and you need to make sure every consumer of that value knows it's a rate in rad/s, not rad/orbit.

Here's the moment it stops being abstract and becomes software — the J2 STM in TypeScript. The paper defines a 6×6 matrix where each row captures how one ROE component evolves under J2 perturbations. Orbital factors like kappa, E, F, G, P, Q, S, T encode how Earth's oblateness affects the orbit — and they compose into matrix elements that track apsidal precession (eccentricity vector rotation) and nodal regression (orbit plane drift):

// From src/orbital/stm/j2.ts — the J2 STM, Equation A6
// Orbital factors (kappa, E, F, G, P, Q, S, T) are precomputed
// from eccentricity, inclination, and J2

return [
  // Row 1: delta-a is constant (no J2 secular effect on semi-major axis)
  [1, 0, 0, 0, 0, 0],

  // Row 2: delta-lambda evolution (Keplerian drift + J2 corrections)
  [
    -(1.5 * n + 3.5 * kappa * E * P) * tau,
    1,
    kappa * ex_i * F * G * P * tau,
    kappa * ey_i * F * G * P * tau,
    -kappa * F * S * tau,
    0,
  ],

  // Row 3: delta-ex evolution (apsidal precession rotates eccentricity vector)
  [
    3.5 * kappa * ey_f * Q * tau,
    0,
    cos_wt - 4 * kappa * ex_i * ey_f * G * Q * tau,
    -sin_wt - 4 * kappa * ey_i * ey_f * G * Q * tau,
    5 * kappa * ey_f * S * tau,
    0,
  ],

  // Row 4: delta-ey evolution (apsidal precession, conjugate to row 3)
  [
    -3.5 * kappa * ex_f * Q * tau,
    0,
    sin_wt + 4 * kappa * ex_i * ex_f * G * Q * tau,
    cos_wt + 4 * kappa * ey_i * ex_f * G * Q * tau,
    -5 * kappa * ex_f * S * tau,
    0,
  ],

  // Row 5: delta-ix is constant
  [0, 0, 0, 0, 1, 0],

  // Row 6: delta-iy evolution (nodal regression)
  [3.5 * kappa * S * tau, 0, -4 * kappa * ex_i * G * S * tau,
   -4 * kappa * ey_i * G * S * tau, 2 * kappa * T * tau, 1],
];
Enter fullscreen mode Exit fullscreen mode

This is a direct translation of the paper's Equation A6. Each element is a product of orbital factors and time — no numerical integration, no timestep loops, just algebra. That's what makes it fast enough to call hundreds of times per second in the targeting solver.

Getting the constants and units right. The J2 coefficient, Earth's gravitational parameter, equatorial radius — these are well-known values, but the paper may normalize differently than your implementation. A factor-of-two error in how you apply J2 produces trajectories that look plausible but are quantitatively wrong. There's no compiler warning for "your orbital factor η should use 1 - e², not 1 - e."

Making the API not terrible. The paper's formulation is beautiful math. The TypeScript API needs to be something a developer can use without reading the paper. That means clear type definitions, sensible defaults, and function signatures that guide you toward correct usage:

export type ManeuverLeg = {
  readonly from: Vector3;
  readonly to: Vector3;
  readonly targetVelocity: Vector3;
  readonly tof: number;
  readonly burn1: Maneuver;
  readonly burn2: Maneuver;
  readonly totalDeltaV: number;
  readonly converged: boolean;
  readonly iterations: number;
  readonly positionError: number;
};
Enter fullscreen mode Exit fullscreen mode

Every field is readonly. The types are self-documenting. You can't accidentally mutate a maneuver leg and corrupt downstream computations. This is the kind of thing MATLAB doesn't even have a concept for.

Validating correctness is the scariest part. Getting a matrix "mostly right" is worse than getting it obviously wrong. A 6×6 STM with one wrong element out of 36 will often produce trajectories that converge — the Newton-Raphson solver is robust enough to compensate — but the trajectories will be physically wrong. The solver finds a minimum. It's just not the minimum. You get a smooth, plausible-looking trajectory to the wrong destination. No error, no divergence, no warning. A wrong index throws. A null pointer crashes. A wrong sign in element (3,4) of a State Transition Matrix gives you a spacecraft that calmly spirals in the wrong direction over six orbits.

The validation approach: implement the same dynamics as an RK4 numerical integrator, propagate using both the STM and the brute-force integrator, and verify they agree within numerical tolerance. This was done in the original Bun implementation with full test coverage. Adam Koenig's direct review of the matrix implementations closes the loop — the person who derived the math checking the port.

Making It Interactive: The 3D Planner

The math is the foundation, but the tool is the point.

The planner uses React Three Fiber for 3D rendering. The chief spacecraft sits at the origin, the deputy moves along computed trajectories, and the user places waypoints by clicking in 3D space. Dragging a waypoint triggers incremental replanning — only the legs affected by the change are recomputed, keeping the interaction fluid even with complex missions.

Three Zustand stores with cross-store subscriptions. Mission state (orbit parameters, waypoints, computed plan), simulation state (playback time, speed), and UI state (sidebar, HUD) are separate stores, with a sync layer that bridges them. The mission store uses subscribeWithSelector so the 3D scene only re-renders when trajectory data changes, not when the user toggles the sidebar.

Synchronous solver on the main thread. The solver doesn't use Web Workers. For 6×6 matrix multiplies, the structured clone serialization overhead of posting data to a worker would likely exceed the computation cost itself. The entire multi-waypoint planning loop — including TOF optimization via golden section search — completes fast enough that it runs between frames.

The solver needs to handle every waypoint a user can click, not just textbook cases. A pure Newton-Raphson solver diverges on many real inputs. Here's what that looks like in practice:

// From src/orbital/targeting/rendezvous.ts — adaptive damping
let dv1Correction: Vector3;
try {
  const jacobianInv = invert3x3(jacobian);
  dv1Correction = matMul3x3_3x1(jacobianInv, positionError);
} catch {
  // Jacobian singular — use gradient descent fallback
  dv1Correction = positionError;
}

let damping: number;
if (iter < 3) {
  damping = 0.5;
} else if (iter < 10) {
  damping = 0.8;
} else {
  damping = 1;
}
dv1 = add3(dv1, [
  damping * dv1Correction[0],
  damping * dv1Correction[1],
  damping * dv1Correction[2],
]);
Enter fullscreen mode Exit fullscreen mode

Aggressive damping early, relax toward full Newton steps as the solver converges, fall back to gradient descent when the Jacobian goes singular. None of this was derived from theory — it came from watching the solver blow up on edge cases and tuning until it converged reliably across all scenario presets.

This Is What Browser-Native Mission Design Looks Like

A mission analyst opens a URL, loads an ISS-altitude scenario, places three waypoints in 3D space, toggles J2 perturbations, and watches the delta-v budget update in real time. No MATLAB license, no installation, no waiting.

The math in these papers is public. The STMs, the Gauss Variational Equations, the targeting algorithms — all published and peer-reviewed. We don't need to wait for legacy vendors to modernize.

TypeScript gives you types that prevent the class of bugs MATLAB can't even detect. Browser deployment eliminates the install. Analytical STMs give you the performance budget for real-time interactivity. Five years ago, this would have required a desktop application and a license. Now it runs in a browser tab with zero infrastructure.

The orbital mechanics library is framework-agnostic and designed to be consumed independently. If you're building anything involving relative spacecraft motion, it might save you from reimplementing Koenig et al. from a MATLAB script on a shared drive.

GitHub: sakobu/koenig-guffanti-damico-roe-stm


What domain-specific math are you sitting on that deserves a serious frontend? If you've hit the MATLAB-to-browser gap in your own field, I'd love to hear what paper you'd port first.

Top comments (0)