We built a Keplerian solar system in 700 lines of WebGL

TL;DR — A pure-Keplerian, time-controllable solar system runs in about 700 lines of TypeScript on top of React Three Fiber. The cheat code is solving Kepler’s equation properly (11 lines of Newton–Raphson) and sharing one solver across 20+ moons, 9 planets, and 2 comets. We’ll walk through the math, the refactor that cut our 3D layer by 73%, and the gotchas (eclipse traps, GPU float precision, time injection).

Live demo: /universe · personal-birthday variant: /universe/birthday · interactive solver: /tools/kepler

---

When we set out to build Universe, we had two constraints in tension:

1. It has to be astronomically real. Saturn at its actual J2000 longitude for any date. The Moon at the right phase. Halley sweeping past on the right century.
2. It has to ship in a few thousand lines. Universe is a feature inside a larger product — not a standalone Stellarium clone.

The cheap path is cos(angle * t) per planet. That looks fine for ten seconds — until somebody asks “where was Mars the day I was born?” and the simulation is off by months.

The right path is Kepler. And it turns out Kepler is short.

This is the engineering tour. What we built, what we cut, and the refactor that took our orbital + 3D layer from ~6,500 lines down to ~1,750.

---

Why circles don’t work

Here’s the trap, in every starter “3D solar system” tutorial on the internet:

// Looks fine. Wrong by months within a year.
const angle = (time / period) * Math.PI * 2;
const x = Math.cos(angle) * radius;
const z = Math.sin(angle) * radius;

Two reasons it’s wrong:

• Orbits are ellipses. Earth swings between 147.1 and 152.1 million km from the Sun. Pluto’s eccentricity (0.25) means its periapsis is closer to the Sun than Neptune’s orbit.
• Orbital speed is not constant. Kepler’s second law — equal areas in equal times — means bodies move fastest at periapsis and slowest at apoapsis. A constant-ω model averages it out and gets the instantaneous position increasingly wrong as eccentricity grows.

For Halley (e ≈ 0.97), the circular model is laughable. For the Moon (e ≈ 0.055), it’s still enough to visibly drift after a few orbits.

The fix is Kepler’s equation:

M = E − e·sin(E)

M is the mean anomaly (linear in time). E is the eccentric anomaly (the geometric angle we want, almost). e is eccentricity. There’s no closed form for E given M. You iterate.

---

The solver, in 11 lines

function solveKepler(meanAnomaly: number, e: number, maxIterations = 8): number {
    let E = meanAnomaly;
    for (let i = 0; i < maxIterations; i++) {
        const f = E - e * Math.sin(E) - meanAnomaly;
        const fPrime = 1 - e * Math.cos(E);
        const delta = f / fPrime;
        E -= delta;
        if (Math.abs(delta) < 1e-10) break;
    }
    return E;
}

Newton–Raphson on f(E) = E − e·sin(E) − M. Converges quadratically. Measured iteration counts on real bodies:

Body	Eccentricity	Iterations to converge
Most moons	< 0.05	1
Earth	0.0167	1
Mercury	0.206	2
Pluto	0.249	2–3
Nereid	0.749	4–5
Halley	0.967	6–7

8 is the budget for the worst case, and the early-out on delta < 1e-10 means we only pay actual cost per body per frame. No body in our system has ever burned the full 8.

A heads-up. Textbooks often suggest E₀ = M + e·sin(M) as the starting estimate to speed convergence for high-e orbits. We tried it. Saves one iteration on Halley, zero on everything else. Not worth the extra line.

---

From eccentric anomaly to a 3D position

Once you have E, the rest is high-school geometry. True anomaly:

const trueAnomaly = 2 * Math.atan2(
    Math.sqrt(1 + e) * Math.sin(E / 2),
    Math.sqrt(1 - e) * Math.cos(E / 2)
);

Distance from focus (the Sun, for planets):

const r = (a * (1 - e * e)) / (1 + e * Math.cos(trueAnomaly));

Position in the orbital plane:

return new THREE.Vector3(
    r * Math.cos(trueAnomaly),
    0,                              // ecliptic plane
    r * Math.sin(trueAnomaly)
);

All orbits go in the XZ plane (the ecliptic). Real inclinations exist — Pluto’s 17°, Halley’s 162° (retrograde!) — but for our use case the silhouette of the solar system is dominated by what’s in the ecliptic. Skipping the 3×3 rotation per body per frame is a real perf win.

Axial tilts and lunar inclination get layered on top in the renderer (the orbital math itself stays planar).

---

The refactor: one solver, twenty moons

Here’s the part you can copy verbatim. Our original layout had ~20 files like this:

io.ts                # 70 lines, copy-pasted Kepler with Io constants
europa.ts            # 70 lines, copy-pasted Kepler with Europa constants
titan.ts             # 70 lines, copy-pasted Kepler with Titan constants
... 17 more ...

~1,400 lines of orbital code for the moon set alone. Almost all of it duplicated, all of it would-have-to-fix-21-times if we found a bug in the solver.

The refactor: one generic module.

// moonOrbit.ts — the solver lives here, exactly once.
export interface MoonOrbitParams {
    periodHours: number;
    semiMajorAxisKm: number;
    eccentricity: number;
    sceneRadius: number;     // visual scaling — see next section
}

export function getMoonPosition(
    date: Date,
    params: MoonOrbitParams
): THREE.Vector3 {
    // solveKepler + geometry — 25 lines total
}

Each moon collapses to a thin data file:

// rheaOrbit.ts — the entire file, all 18 lines
import { getMoonPosition, MoonOrbitParams } from "./moonOrbit";

const RHEA_PARAMS: MoonOrbitParams = {
    periodHours: 108.4,
    semiMajorAxisKm: 527_040,
    eccentricity: 0.001,
    sceneRadius: 0.9,
};

export const getRheaPosition = (date: Date) =>
    getMoonPosition(date, RHEA_PARAMS);

Per-moon files went from ~70 lines to ~18 — a 75% drop across the moon set. We did the same exercise for the mesh/texture loading layer (one generic Moon3D + per-moon data), and across the whole 3D + orbital layer:

6,470 → 1,754 lines. A 73% reduction without removing a single feature.

Why keep one file per moon instead of one big table? Because the data is per-moon, the comments are per-moon (“Io’s libration is driven by Jupiter, not its own orbit”), and the import graph stays clean — a component that only renders Rhea pulls in only Rhea’s params.

---

Scene units vs reality

We use scene units internally (~1.0 for “interesting size”) and convert at the data boundary. Rhea’s real semi-major axis is 527,040 km. Its sceneRadius is 0.9. The conversion lives inside the solver:

const sceneScale = sceneRadius / semiMajorAxisKm;
const rScene = r * sceneScale;

Two reasons you want this:

1. Editorial control over visual scaling. At true-to-scale, the moons are pinpricks next to Saturn, and Saturn is a pinprick next to its orbit. Designers pick sceneRadius; astronomers pick eccentricity and periodHours. Both win — the orbit is correctly shaped and correctly timed, just sized for the camera.

2. GPU float precision. Three.js uses 32-bit floats on the GPU. Working in scene units near 1.0 stays in the precise range. Naively dumping kilometers into Three.js (1.5e8 for Earth’s distance from the Sun) and you hit the float precision cliff — bodies start to jitter at certain camera angles, geometry z-fighting appears, the camera near plane stops being where you think it is.

You can spend a week debugging it. Or you can convert at the boundary and move on.

---

Time as an injected dependency

A working orbit module isn’t enough. Users want to scrub time, pause on their birthday, hit play and watch planets move. We have one global TimeManager:

class TimeManager {
    private current = new Date();
    private paused = false;
    private speed = 1.0;

    setDate(d: Date) { this.current = new Date(d); }
    pause() { this.paused = true; }
    resume() { this.paused = false; }
    setSpeed(s: number) { this.speed = s; }

    tick(deltaMs: number) {
        if (this.paused) return;
        this.current = new Date(this.current.getTime() + deltaMs * this.speed);
    }

    now() { return this.current; }
}

Every body’s render frame pulls timeManager.now() and feeds it to getMoonPosition / getPlanetPosition. The whole pipeline is functionally pure from Date → Vector3, so freezing on a specific date is literally timeManager.pause() after setDate(...).

This is what powers our Birthday Sky feature: type your birthday, the time manager snaps to that moment, the renderer redraws every body from Kepler’s equation, you see the solar system exactly as it was the day you were born. No special-case code path — same render loop, same now().

Try it on your own birthday →

---

Camera direction in 3D space

The orbital math is one thing. Framing it cinematically is another. Birthday Sky has two camera modes:

• Heliocentric — high oblique shot from behind Earth. Earth as foreground disk. Sun glowing in the background. Inner planets falling naturally in frame.
• Earth POV — just above Earth’s night-side surface, looking outward into deep space.

Each is a handful of vector operations:

// Helio composition.
side.copy(up).cross(sunToEarth).normalize();
outPos
    .copy(earthPos)
    .addScaledVector(sunToEarth, 5)   // 5 units behind Earth (anti-Sun axis)
    .addScaledVector(up, 6)           // 6 units above the ecliptic
    .addScaledVector(side, 2);        // small lateral — adds depth
outLook.set(0, 0, 0);                 // look at the Sun

The eclipse trap

It’s tempting to put the camera directly behind Earth on the anti-Sun line. Don’t.

Earth’s apparent disk from 5 units away (~5° angular diameter, given Earth’s render radius of 0.6 — 0.8 in scene units) is bigger than the Sun’s apparent disk at 19 units (~3–4°). Camera directly behind Earth = Earth fully eclipses the Sun, no light source visible, dark dead frame.

The 2-unit lateral offset and 6-unit vertical lift give you ~25° angular separation between Earth center and Sun center as seen from camera. Earth as foreground disk, Sun glowing past Earth’s edge. Iconic composition, intentional geometry.

We spent two debug sessions getting this right. The math:

Earth angular radius from 8 units away = atan(0.8 / 8) ≈ 5.7°
Sun angular radius from 14.5 units away = atan(1.0 / 14.5) ≈ 4.0°
Required angular separation > 5.7 + 4.0 = 9.7° to avoid overlap.

We aim for ~25° to leave breathing room.

---

Things we deliberately didn’t do

The honest part. Cutting these is what kept us at 700 lines:

• Perturbations. Real planets don’t follow pure Keplerian orbits — Jupiter tugs on Mars, the Sun’s quadrupole moment matters for Mercury. We ignore all of it. For a 100-year window centered on J2000, visual error is negligible. For asteroid orbit prediction, you’d need numerical integration with N-body forces.
• Light-time correction. When you “see” Saturn from Earth, you’re seeing where it was ~80 minutes ago. We don’t correct. Doesn’t matter for a visualization; would matter for astrometric tools.
• Precession of the equinoxes. Earth’s axis wobbles with a 26,000-year period. We hold the J2000 tilt fixed. Set the date to 12,000 AD and our Polaris is no longer where the pole should be. We’re OK with that.
• Orbital inclination. Every orbit goes in the ecliptic. Real Pluto cranks 17°. We don’t render the tilt.

The rule: the simulation should look right on the days users care about — their birthday, today, tomorrow — and degrade gracefully toward absurd dates.

---

What it all weighs

The final tally for the math + time + camera layer:

Module	Lines
moonOrbit.ts (generic solver)	82
planetOrbit.ts (generic, planets + comets)	176
earthOrbit.ts (axial tilt, ECI conversion)	302
20 moon wrappers × ~18 lines	~360
Comet orbits (Halley, Encke)	~120
TimeManager + scrubber controls	~140
Total	~1,180

About 700 of those are the core orbit code. The rest is time management and the comet special cases.

The 3D rendering layer (sphere meshes, textures, atmosphere shaders, moon meshes, ISS via SGP4, etc.) sits on top of it: another ~1,800 lines built on React Three Fiber. R3F is the right call here — it lets you describe the scene declaratively and reuse React’s diffing for object lifecycles, while still giving you useFrame for the per-tick imperative math.

---

Lessons for your past self

1. Solve Kepler’s equation properly from day one. 11 lines. The difference between a demo and a tool. Don’t ship the circular-orbit version “just to start” — you’ll spend the same time on it and have to throw it out.
2. Pick a scene-unit scale, convert at the boundary. Don’t pass real kilometers around inside Three.js. The float precision cliff is the worst kind of bug — intermittent, scene-angle dependent, and you’ll blame your math first.
3. Make time a globally injected dependency. Every position function takes a Date. Every system reads from one TimeManager. Pause / scrub / replay become free.
4. Share the solver, parameterize the bodies. The line count reduction is nice; the bug-fix reduction is huge. A correctness fix in solveKepler propagates to all 20 moons automatically.
5. Visual scaling is not a bug. It’s the difference between “scientifically accurate but invisible” and “real where it matters and gorgeous everywhere else.” Pick render scales by hand. Keep periods and eccentricities honest.
6. Camera framing is geometry, not vibes. Calculate the angular size of your subjects from the camera, calculate the required separation, then place the camera. Eyeballing “looks about right” sends you on debug spirals.

---

Try it

The live thing lives at /universe — scrub time, toggle layers, focus any body. The personal variant (/universe/birthday) snaps the simulation to a moment you care about and adds the cultural context (sun sign, moon phase, the 24 solar terms, the traditional Chinese 时辰).

If you’ve built something similar and have questions — about the solver, the time injection, the camera director, or the refactor — happy to dig in. Reach out via Contact.