N-body simulator

Explore a few interacting particles. Start with the two-body circular orbit, inspect its energy and momentum, then change a mass or velocity. Start/resume continues a paused run. One step pauses and advances exactly one Δt; Reset rebuilds the state from the editable initial values. A physics edit stops the run. Selecting a scene replaces all initial values, force constants, custom formulas and overrides. Switching dimensions preserves hidden coordinates, velocities and formula drafts; only the selected dimensions evolve. Custom pair formulas survive switching to another force preset and back.

Units and forces. Choose consistent mass, length, time and charge units. These are model units, not automatically SI. Positions use length, velocities length/time, forces mass·length/time² and energy mass·length²/time². G has units length³/(mass·time²). Coulomb k has units energy·length/charge²; spring k has units energy/length²; Lennard–Jones k is an energy and σ is a length. ε is the Gravity/Coulomb softening length. A 1D or 2D run constrains motion to a line or plane while retaining the same point-force law.

For a pair, dx=x−xj and similarly for y/z. r²=dx²+dy²+dz² and r=√r² are geometric distances. Gravity uses F=−G·m·mj·Δr/(r²+ε²)^(3/2), Coulomb uses F=k·q·qj·Δr/(r²+ε²)^(3/2), and pair springs use F=−k·Δr with zero rest length. Lennard–Jones uses U=4k[(σ/r)^12−(σ/r)^6] and F=24k[2(σ/r)^12−(σ/r)^6]·Δr/r²: repulsion at short range, attraction beyond r=2^(1/6)σ. It uses no softening or mass factor. Singular overlaps or undefined expressions stop the run instead of producing a zero force.

Custom expressions. Components are forces, divided by the receiving mass to obtain acceleration. Pair variables: t, m, mj, q, qj, dx, dy, dz, r, r2, x, y, z, xj, yj, zj, G, k, eps, sigma. Global/local field variables: t, m, q, x, y, z, G, k. An override replaces all pair forces received by that particle; it does not replace the force received by its partner. Local external forces add to global forces. Use Fx; Fy; Fz in override editors; omitted components are zero and a wholly blank pair override keeps the shared law. Inactive dimensions are retained but ignored. For example a global force 2*m; −m; 3*m gives constant acceleration (2; −1; 3) in 3D; enter ordinary ASCII minus signs in formulas.

Supported operators are +, -, *, /, %, ^, **, comparisons, &&, ||, ! and condition ? yes : no. Both powers associate to the right: 2^3^2=512 and −2^2=−4. Conditional branches short-circuit, so r>0 ? 1/r : 0 is valid. Pure mathematical functions include abs, sqrt, cbrt, pow, exp, log/log10/log2, trigonometric and hyperbolic functions and their inverses, atan2, hypot, min, max, floor, ceil, round, trunc, sign, expm1 and log1p. Math.pow and other allowed Math-qualified names also work; constants include PI and E. There are no globals, assignments, object access or random(). Each component allows 1,024 characters, 256 syntax nodes and nesting depth 32. Nonfinite intermediate calculations are errors.

Reading results. Kinetic energy, momentum and centre of mass are always available. Total energy and E−E₀ are shown only for a known pair preset without active additional fields or overrides; arbitrary custom formulas have no inferred potential. For directed overrides or external fields, momentum need not be conserved. The numerical method is fixed-step velocity Verlet. Smaller Δt generally improves accuracy; close encounters and stiff repulsion may need much smaller steps. This is an educational few-particle model, without collision merging or an adaptive solver.

The circular scene has G=m=1, radius 0.5 and speed √0.5. The triangle uses the rotation speed appropriate to its softening. Figure-eight data are rounded published initial conditions, so exact long-term recurrence is not promised; the perturbed scene is an illustration, not a chaos classifier. The plot offers orthographic XY/XZ/YZ projections in 3D and equal scales. Fit, projection and trail visibility do not change model time. Trails retain physical xyz and time, at most 500 samples per particle, spaced by at least 0.02 model time units.

Limits and saving. There are 2–12 particles, Δt=0.0001–0.1 and playback 0.01–10 model time units/second. Masses are 10⁻⁶–10⁶; initial coordinates, velocities and charges are bounded by ±10⁶; G and k by 0–10⁶; ε by 0–1,000; σ by 10⁻⁶–1,000. Up to 32 integration steps are performed per frame, with at most 0.05 seconds of wall time accepted; the run slows and reports the limit instead of accumulating a backlog. Hidden tabs/tools stop the model clock. With reduced motion, the first Start prepares a paused state; a second Start explicitly animates. Only input drafts and view settings are serialized, not trajectories. Drafts up to the platform’s 1,500-character state limit can survive reload; larger drafts stay in module memory during navigation and are lost on reload, with a visible notice.

References: LAMMPS 12–6 potential, velocity Verlet, Plummer softening, Montgomery’s figure-eight data.