Numerical Insights: Integration, Drift, and Timestep Trade-offs

This document covers three topics that came out of analysing the Sun-Venus simulation:

1. How Venus's position is actually calculated — step by step from browser frame to orbital mechanics
2. What "Kepler drift" and "Horizons drift" measure — and why they are fundamentally different things
3. What happens when you change the substep size — accuracy vs performance trade-off table and reasoning

1. How Venus's Position Is Updated Each Frame

There are four distinct layers between a browser animation frame and a new planetary position. Understanding each layer is what allows the Python drift analysis script to reproduce the simulation results exactly.

Layer 1 — loop(now) — frame-rate dependent

const realDt = Math.min((now - lastRAF) / 1000, 0.1);  // real seconds since last frame
const simDt  = realDt * REAL_TO_SIM * timeWarp;         // convert to sim-seconds
advanceSim(simDt);

At 60 fps with Warp 1x: realDt ≈ 0.0167 s → simDt ≈ 10,080 sim-seconds (~2.8 sim-hours).

The frame rate does not affect physics accuracy. It only controls how many substeps are batched per render call. A slow machine that renders at 10 fps still runs the same physics as a 120 fps machine — it just runs more substeps per frame.

Layer 2 — advanceSim(totalDt) — substep batching

const steps = Math.ceil(totalDt / DT_SUB);   // DT_SUB = 3600 s
const dt    = totalDt / steps;                // always ≤ 3600 s per substep

This is the accuracy-controlling layer. Regardless of frame rate, time warp, or how much sim-time has accumulated, each individual physics step is at most 3600 sim-seconds (1 sim-hour). At Warp 50x, one frame advances ~500,000 sim-seconds, which becomes ~139 substeps of equal size.

Layer 3 — Kick-Drift-Kick Leapfrog (per substep)

A.  v += a · (dt/2)          ← half-kick: velocity update with current acceleration
B.  x += v · dt              ← full drift: position update with new velocity
    simTime += dt
    record position to history buffer if interval elapsed
C.  computeForces()          ← recompute acceleration at the new position
D.  v += a · (dt/2)          ← second half-kick with new acceleration

Why leapfrog and not Euler? Euler integration does not conserve energy — orbits spiral outward or collapse over timescales of years. The leapfrog is symplectic: it conserves a slightly modified total energy exactly, so orbits remain bounded over arbitrarily long timescales. The only accumulating error is a slow phase drift (the orbit period is slightly wrong), not an amplitude drift (the orbit size stays correct).

Layer 4 — computeForces() — light-delayed gravity

d0         = |r_Venus − r_Sun|      ← current separation (~108 million km for Venus)
lightDelay = d0 / c                  ← ~360 seconds (6 minutes)
t_ret      = simTime − lightDelay   ← when the signal left the Sun

sPast      = posHist interpolated at t_ret   ← Sun's position 360 s ago

F_on_Venus = G · M_Sun / |r_Venus − sPast|²  ← Newton's law with retarded source
             directed toward sPast

The position history buffer stores one entry every HIST_INTERVAL = 3600 s of sim time, keeping the last MAX_HIST_DURATION = 7200 s. The retarded time lookup is a binary search with linear interpolation.

Summary. Venus's position at any moment is the result of pure numerical integration starting from the Horizons or J2000 initial state vectors. There is no orbit equation. The ellipse emerges from Newton's law — applied to the retarded gravitational field — iterated one hour at a time.

2. Venus's Orbit: Elliptical, Not Circular

The orbit is genuinely elliptical (eccentricity e ≈ 0.0067). It looks nearly circular because Venus has one of the most circular orbits in the solar system, but the simulation does not hardcode a circle anywhere.

The J2000 built-in initial conditions are derived from e = 0.00677323 (IAU mean elements). The Kepler propagator explicitly computes the eccentricity vector from the initial state and propagates along the correct conic section. The leapfrog integrator produces whatever orbit shape naturally follows from Newton's law applied to the initial position and velocity.

3. Kepler Drift vs Horizons Drift — Two Different Things

The simulation's Reference Ephemeris panel shows two drift values. They measure completely different sources of error. Conflating them leads to wrong conclusions about the simulation's accuracy.

Kepler drift

Definition: Distance between Venus's simulated position and where a perfect 2-body Newtonian simulation would put it, given the same starting conditions.

Source: The light-delayed gravity perturbation. When gravity travels at c rather than instantaneously, the force direction is slightly off from the current Sun-Venus line. Over many orbits, this introduces a slow precessional drift relative to the ideal Keplerian ellipse. The leapfrog's finite-timestep phase error also contributes.

Sensitive to substep size: Yes. Smaller substeps give less integration error and smaller Kepler drift. This is the metric that changes meaningfully between dt = 100 s and dt = 20,000 s.

Horizons drift

Definition: Distance between Venus's simulated position and where JPL Horizons says Venus actually is on that date.

Source: Primarily the missing planetary perturbations. The simulation models only the Sun and Venus. Jupiter exerts a gravitational acceleration on Venus of roughly G × M_Jupiter / r² ≈ 2×10⁻⁸ km/s² — tiny per second, but accumulating continuously. Saturn, Earth, Mercury, and Mars also contribute. After decades, the orbit phase drifts measurably relative to where Venus actually is.

Not sensitive to substep size: The Horizons drift curve looks essentially identical whether you run at dt = 100 s or dt = 20,000 s, because the missing-perturbation error is orders of magnitude larger than the integration error at any reasonable substep size.

Summary

This asymmetry is why the drift analysis uses both metrics: the Horizons drift shows the hard physical limitation of the 2-body model, while the Kepler drift reveals the numerical precision of the integration itself.

4. Substep Size Trade-offs

The default substep is DT_SUB = 3600 s (1 sim-hour). The drift analysis script benchmarks six values: 100, 1000, 3600, 7200, 14000, and 20,000 seconds.

The leapfrog's per-orbit phase error scales as O(dt²). For Venus's 224.7-day orbital period:

T_Venus = 224.7 days = 19,414,080 s

dt =   100 s  →  dt/T = 1/194,141   ← sub-km error per orbit, negligible
dt =  1000 s  →  dt/T = 1/19,414    ← very small
dt =  3600 s  →  dt/T = 1/5,393     ← default; small, well-characterised
dt =  7200 s  →  dt/T = 1/2,696     ← 4× more error than default
dt = 14000 s  →  dt/T = 1/1,387     ← measurable Kepler drift over decades
dt = 20000 s  →  dt/T = 1/971       ← significant Kepler drift

Keep dt = 3600 s for general use. It provides ~1.75 million steps, Kepler drift that genuinely reflects the light-delay physics, and position-history entries spaced tightly enough (3600 s) to accurately resolve the 360 s light delay.

5. Parallelisation vs Larger dt — Big-O

The leapfrog is inherently serial — step N depends on step N-1's positions and velocities. You cannot compute step 1000 without first computing steps 1–999. For a single run, increasing dt by a factor k gives O(1/k) speedup with no parallelisation alternative that preserves accuracy.

When running N simulations at different dt values — as in the drift analysis — each run is completely independent. This is embarrassingly parallel: no communication, no coordination. Running all N simultaneously on N cores gives O(1/N) wall time.

Larger dt          →  O(1/k) speedup  but  accuracy cost per run  (Kepler drift increases)
Parallelisation    →  O(1/P) speedup  and  no accuracy cost        (each run uses its own dt)

They are big-O equivalent only when comparing "one coarse run" against "P parallel fine runs". For independent multi-run benchmarks, parallelisation wins because it preserves accuracy while reducing wall time.

6. Why the Drift Analysis Uses Python, Not the Browser Simulation

The browser simulation is unsuitable for long-baseline logging: the NASA API requires a CORS proxy, there is no mechanism to export drift data to a file from JavaScript in a browser tab, and the browser cannot run a 200-year simulation without user interaction to keep it alive.

The Python script reproduces the JS physics numerically exactly because all physical constants are identical, the coordinate transform uses the same obliquity and rotation formula, the leapfrog algorithm is the same kick-drift-kick sequence, and the Horizons API is called with the same parameters.