Planets do not move at a constant speed and their orbits are not circles. Kepler worked that out from Tycho's tables alone, decades before anyone could explain why.
Move the cursor over the canvas — it is part of the simulation.
Drag your finger across the canvas — it is part of the simulation.
What to look for
Your cursor is a second mass. Bring it near an orbit and you can pull a planet into a completely different one, or throw it out of the system entirely. That instability is the three-body problem in miniature.
The maths
r = a(1 − e²)/(1 + e·cos θ)
First law: the orbit is an ellipse with the star at one focus, not at the centre. e is the eccentricity the slider controls.
dA/dt = constant
Second law: equal areas in equal times. The shaded wedges are drawn at a fixed number of frames apart, so they are all the same area even though they look different.
T² ∝ a³
Third law. It is why the outer planets in the simulation crawl and the inner ones sprint, without that ever being programmed in.
How it is built
The three laws are not coded. Only Newtonian gravity is, integrated with leapfrog (kick-drift-kick) at 64 sub-steps per frame, and the ellipses, the varying speed and the period ratios all fall out of it. Leapfrog conserves angular momentum exactly, which is why the swept wedges, drawn from the stored path at a fixed frame spacing, have equal areas to machine precision.
The code, step by step
Cut straight from lab.js. These are the real lines that run above, not a simplified version.
01
Only Newton is coded
const acel = () => {
const dx = cx - b.x, dy = cy - b.y, q = dx * dx + dy * dy + 0.25;
const k = GM / (q * Math.sqrt(q));
let ax = dx * k, ay = dy * k;
if (M.in && P.perturb) { // el cursor es otra masa
const ux = M.x - b.x, uy = M.y - b.y, ur = Math.hypot(ux, uy) + 14;
ax += (ux / ur) * (P.perturb / (ur * ur));
ay += (uy / ur) * (P.perturb / (ur * ur));
}
return [ax, ay];
};
const H = 1 / 64;
let a = acel();
for (let s = 0; s < 64; s++) {
b.vx += a[0] * H / 2; b.vy += a[1] * H / 2;
b.x += b.vx * H; b.y += b.vy * H;
a = acel();
b.vx += a[0] * H / 2; b.vy += a[1] * H / 2; // segundo medio impulso
Inverse-square attraction toward the star with a half-pixel Plummer softening, integrated with leapfrog at 64 sub-steps per frame. Kepler's three laws are never written anywhere: the ellipses, the varying speed and the period ratios emerge from these lines. Leapfrog conserves angular momentum exactly and the orbits close; the earlier scheme added 6 px to the distance and made them precess 60° per revolution.
02
Your cursor is a second mass
if (M.in && P.perturb) {
ctx.strokeStyle = `rgba(${acc},0.4)`; ctx.lineWidth = 1;
ctx.beginPath(); ctx.arc(M.x, M.y, 4 + P.perturb / 900, 0, TAU); ctx.stroke();
}
const n = P.planets | 0;
for (let i = 0; i < n && i < ps.length; i++) {
const b = ps[i];
// Leapfrog, 64 subpasos por frame, suavizado de Plummer de 0,5 px. Con Euler
// y la distancia +6 px la órbita precesaba 60° por vuelta (suite de
// validación de BN Lab Simulations); así baja de medio grado.
const acel = () => {
const dx = cx - b.x, dy = cy - b.y, q = dx * dx + dy * dy + 0.25;
const k = GM / (q * Math.sqrt(q));
let ax = dx * k, ay = dy * k;
if (M.in && P.perturb) { // el cursor es otra masa
const ux = M.x - b.x, uy = M.y - b.y, ur = Math.hypot(ux, uy) + 14;
ax += (ux / ur) * (P.perturb / (ur * ur));
ay += (uy / ur) * (P.perturb / (ur * ur));
It contributes its own inverse square term. Two masses is solvable in closed form; adding this third one is not, which is why the orbits go unpredictable so easily.
03
Draw the swept areas
for (let k = b.p.length - 1; k > b.p.length - 60 && k > 11; k -= 12) {
ctx.beginPath(); ctx.moveTo(cx, cy);
ctx.lineTo(b.p[k][0], b.p[k][1]); ctx.lineTo(b.p[k - 11][0], b.p[k - 11][1]);
ctx.closePath(); ctx.fill();
The wedges are taken every twelve stored frames, so they cover equal times by construction. They look wildly different in shape and are the same area: that is the second law, measured rather than asserted.
Why it matters
Kepler spent years fitting circles to Mars and failing by eight arcminutes. He trusted the data over the shape everyone assumed was perfect, and that decision is where modern astronomy starts.
The whole thing
Everything above, in one piece. This is the entire experiment as it lives in lab.js.
FULL
lab.js · kepler
{ id: 'kepler',
name: L('Kepler orbits', 'Órbitas de Kepler'),
note: L('Equal areas in equal times. The law Newton later explained.',
'Áreas iguales en tiempos iguales. La ley que Newton explicó después.'),
params: [
{ key: 'ecc', label: L('Eccentricity', 'Excentricidad'), min: 0, max: 0.9, step: 0.01, def: 0.55 },
{ key: 'planets', label: L('Planets', 'Planetas'), min: 1, max: 6, step: 1, def: 3 },
{ key: 'perturb', label: L('Cursor mass', 'Masa del cursor'), min: 0, max: 4000, step: 100, def: 900 }
],
make(w, h) {
let ps = [];
const seed = (w, h, ecc) => {
ps = [];
const cx = w / 2, cy = h / 2;
for (let i = 0; i < 6; i++) {
const r = Math.min(w, h) * (0.12 + i * 0.055);
const vc = Math.sqrt(2600 / r) * Math.sqrt(1 - ecc);
ps.push({ x: cx + r, y: cy, vx: 0, vy: vc, p: [], area: [] });
}
};
let lw = w, lh = h;
seed(w, h, 0.55);
return {
resize(w, h) { lw = w; lh = h; seed(w, h, 0.55); },
reset() { seed(lw, lh, 0.55); },
step(ctx, w, h, t, acc, P, M) {
const cx = w / 2, cy = h / 2, GM = 2600;
ctx.fillStyle = `rgba(${acc},0.95)`;
ctx.beginPath(); ctx.arc(cx, cy, 6, 0, TAU); ctx.fill();
if (M.in && P.perturb) {
ctx.strokeStyle = `rgba(${acc},0.4)`; ctx.lineWidth = 1;
ctx.beginPath(); ctx.arc(M.x, M.y, 4 + P.perturb / 900, 0, TAU); ctx.stroke();
}
const n = P.planets | 0;
for (let i = 0; i < n && i < ps.length; i++) {
const b = ps[i];
// Leapfrog, 64 subpasos por frame, suavizado de Plummer de 0,5 px. Con Euler
// y la distancia +6 px la órbita precesaba 60° por vuelta (suite de
// validación de BN Lab Simulations); así baja de medio grado.
const acel = () => {
const dx = cx - b.x, dy = cy - b.y, q = dx * dx + dy * dy + 0.25;
const k = GM / (q * Math.sqrt(q));
let ax = dx * k, ay = dy * k;
if (M.in && P.perturb) { // el cursor es otra masa
const ux = M.x - b.x, uy = M.y - b.y, ur = Math.hypot(ux, uy) + 14;
ax += (ux / ur) * (P.perturb / (ur * ur));
ay += (uy / ur) * (P.perturb / (ur * ur));
}
return [ax, ay];
};
const H = 1 / 64;
let a = acel();
for (let s = 0; s < 64; s++) {
b.vx += a[0] * H / 2; b.vy += a[1] * H / 2;
b.x += b.vx * H; b.y += b.vy * H;
a = acel();
b.vx += a[0] * H / 2; b.vy += a[1] * H / 2; // segundo medio impulso
}
b.p.push([b.x, b.y]); if (b.p.length > 700) b.p.shift();
ctx.strokeStyle = `rgba(${acc},0.35)`; ctx.lineWidth = 1;
ctx.beginPath();
b.p.forEach((q, k) => k ? ctx.lineTo(q[0], q[1]) : ctx.moveTo(q[0], q[1]));
ctx.stroke();
// Segunda ley: sectores barridos en intervalos iguales
if (i === 0 && b.p.length > 12) {
ctx.fillStyle = `rgba(${acc},0.14)`;
for (let k = b.p.length - 1; k > b.p.length - 60 && k > 11; k -= 12) {
ctx.beginPath(); ctx.moveTo(cx, cy);
ctx.lineTo(b.p[k][0], b.p[k][1]); ctx.lineTo(b.p[k - 11][0], b.p[k - 11][1]);
ctx.closePath(); ctx.fill();
}
}
ctx.fillStyle = `rgba(${acc},0.95)`;
ctx.beginPath(); ctx.arc(b.x, b.y, 3, 0, TAU); ctx.fill();
}
}
};
}
},
