Devlog #89 – Wave 68: Barnes-Hut, Marching Squares, Poisson Disk, Bell Inequality, Chua Circuit & Pedestrian Crowd

New Simulations

Barnes-Hut N-Body: Quadtree Gravity

Direct N-body gravity is beautiful but brutal: every particle must interact with every other, giving O(n²) force evaluations per time step — n(n−1)/2 pairs for n=2000 means nearly 2 million calculations per frame at 60 fps. Barnes-Hut slashes this to O(n log n) using a recursive quadtree (octree in 3D).

The θ-Criterion

Each internal node stores the total mass and centre of mass of all contained particles. When computing the force on a particle, instead of recursing into a node you replace it with its aggregate if the ratio s/d satisfies the condition s² < θ²d² — where s is the node's side length and d is the distance from particle to node centre of mass. Smaller θ means more accuracy but more evaluations; θ = 0.5 is the standard trade-off.

Velocity-Verlet Integration

Positions and velocities are advanced with the Velocity-Verlet scheme: x(t+dt) = x(t) + v(t)dt + ½a(t)dt², then forces are re-evaluated, and finally v(t+dt) = v(t) + ½(a(t)+a(t+dt))dt. This is a symplectic integrator — it conserves a shadow Hamiltonian, preventing artificial energy drift over long runs. The tree is rebuilt every frame since particles move continuously.

• 4 presets: Gaussian cloud, two-galaxy collision, ring + central mass, central attractor + cloud
• Click any particle to highlight its quadtree decomposition path
• Live counter: BH evaluations vs naive n(n−1)/2
• N slider up to 2,000 particles

Marching Squares: Iso-Contour Lookup Table

Marching Squares extracts iso-contours from a 2D scalar field by processing each grid cell independently. Each cell has four corners — each either above or below the iso-value — giving 2⁴ = 16 distinct topological cases.

The 16-Case Table and Linear Interpolation

A precomputed lookup table maps the 4-bit corner bitmask to edge intersection pairs. For each active edge, the intersection is placed using linear interpolation: t = (iso − v₀) / (v₁ − v₀) along the edge, producing smooth contours rather than jagged pixel boundaries.

Saddle Point Disambiguation

Cases 5 and 10 are ambiguous — the diagonal connectivity is undefined. The simulation resolves these by sampling the bilinear centre value: if the centre is above threshold the contour connects one way; below, the other. This prevents topological artefacts like phantom disconnections in metaball fields.

• Four scalar fields: metaballs (sum of 1/(r²+ε)), fBm noise (6 octaves), animated sine, user-paintable
• Multi-level contouring: draw up to 8 simultaneous iso-values with distinct colours
• Bilinear heatmap overlay shows the raw field beneath the contours

Poisson Disk Sampling: Bridson's Algorithm

Poisson disk sampling produces point sets where no two points are closer than a minimum distance r, while filling space as densely as possible. The result is "blue noise" — a power spectrum that is flat in the mid-frequency range, making it perceptually ideal for dithering, stippling, and terrain generation.

O(n) Background Grid

Bridson's 2007 algorithm achieves O(n) (linear in output size) via a background grid with cell size r/√2. Each cell holds at most one sample; neighbour lookups require checking only a fixed 5×5 window. The active list grows as each accepted sample spawns up to k=30 candidate annular tries in the range [r, 2r].

Blue Noise Verification

The radial power spectrum of the sample set is shown alongside random uniform sampling and a regular grid. Poisson disk displays a characteristic flat plateau in the mid-range — this is the defining blue-noise property, indicating absence of low-frequency clumping and high-frequency aliasing.

• Animate generation step-by-step or batch-complete instantly
• Adjust minimum radius r from 8 to 40 pixels
• Up to ~2000 samples in a 800×600 canvas
• Toggle radial power spectrum overlay for noise quality comparison

Bell Inequality: Quantum Correlations Beat Classical Limits

John Bell's 1964 theorem proved that no local hidden-variable theory can reproduce all quantum mechanical predictions. The CHSH (Clauser-Horne-Shimony-Holt) inequality makes this experimentally testable: measure the combination S = |E(a,b) − E(a,b') + E(a',b) + E(a',b')|.

Quantum vs Classical Bounds

For a singlet state |ψ⁻⟩ = (|01⟩ − |10⟩)/√2, the quantum correlation function is E(a,b) = −cos(a − b). With optimal measurement angles a=0°, a'=45°, b=22.5°, b'=67.5°, this gives S = 2√2 ≈ 2.828 — the Tsirelson bound, the maximum quantum mechanics allows. Any classical hidden-variable model is constrained to |S| ≤ 2.

Live Experiment

The simulation samples random λ values for the classical model and quantum random rotations for the quantum model. Both S values are tracked over N trials and converge toward their theoretical limits as N grows, visually demonstrating the gap that no classical model can bridge.

• Slider to vary measurement angles and watch S change in real time
• Histogram of per-pair products E = A·B accumulates over trials
• Classical and quantum S values displayed with running average

Chua Circuit: The Simplest Chaotic Circuit

Leon Chua's 1983 circuit is the canonical laboratory for electronic chaos. It produces a double-scroll strange attractor — two interleaved spirals in phase space — from just three passive components and one nonlinear resistor (Chua's diode).

The ODEs and Chua's Diode

The dynamics are governed by three coupled ODEs:

dV_C1/dt = (1/C1) · [G(V_C2 − V_C1) − f(V_C1)]
dV_C2/dt = (1/C2) · [G(V_C1 − V_C2) + I_L]
dI_L/dt  = −(1/L) · V_C2

where f(V) is the piecewise-linear I-V of Chua's diode: f(V) = m₁V + ½(m₀−m₁)(|V+1|−|V−1|). The simulation integrates with RK4 at a fixed timestep, rendering trajectories directly onto the canvas.

Lyapunov Exponent and Bifurcation

The largest Lyapunov exponent λ₁ is estimated via the two-trajectory Benettin method: a second trajectory starts infinitesimally offset, and its log-divergence is averaged and renormalised periodically. λ₁ > 0 confirms chaos. As the parameter C1 is swept, the simulation shows period-1 → period-2 → period-4 → chaos (period-doubling cascade).

• Phase portraits: (V_C1, V_C2) plane and (V_C1, I_L) plane side by side
• Real-time Lyapunov exponent readout
• Parameter sweep slider to observe period-doubling

Pedestrian Crowd: Social Force Model

The Helbing-Molnar social force model (1995) treats pedestrian motion as Newtonian dynamics with psychosocial forces. Each agent has a desired velocity toward a goal and experiences repulsion from other pedestrians and walls.

Force Components

The total force on agent i is the sum of three terms:

• Motivation: F_motivation = m·(v₀·ê − v) / τ — drives agent toward desired speed v₀ in desired direction ê with relaxation time τ
• Social repulsion: exponential decay A·exp((r_ij − d_ij)/B)·n̂_ij — soft repulsion with contact hard-core correction when agents overlap
• Wall repulsion: same exponential form applied to the nearest wall point

Emergent Phenomena

The magic of the social force model is what emerges without explicit rules. In a bidirectional corridor, opposing streams spontaneously segregate into stable lanes — minimising head-on encounters. At a narrow bottleneck, a pressure-induced arch forms at the exit, intermittently releasing bursts. In panic mode, raising desired speed (faster-is-slower effect) decreases throughput as arching becomes severe and agents block each other harder.

• Up to 200 pedestrians; spatial hash grid for O(n) neighbour queries
• Three scenario presets: corridor, bottleneck, open room
• Toggle panic mode to demonstrate faster-is-slower
• Flow-rate counter: agents per second through exit

Up Next

Wave 69 arrives with convex hull algorithms (Graham scan vs Jarvis march vs Quickhull compared side by side), Grover's quantum search with geometric amplitude visualisation, river meander dynamics with oxbow lake formation, Verlet rope physics, the classic Rosenblatt perceptron, and Kepler orbital resonance with conjunction clustering.