Devlog #87 – Wave 66: Fourier Epicycles, Magnetic Pendulum, Surface Wetting, Gear Train, Knapsack & Huffman

New Simulations

🌀 Fourier Epicycles — Drawing with the DFT

From drawn path to complex DFT

When you draw a closed curve on the canvas the simulator samples it at N equally-spaced parameter values, treating each sample point as a complex number z_n = x_n + i·y_n. The complex Discrete Fourier Transform is computed:

X[k] = Σ_{n=0}^{N-1}  z_n · exp(−2πi·k·n / N)    for k = 0, 1, …, N−1

Each bin k represents one epicycle:
  centre frequency : k  rotations per full drawing period
  amplitude        : |X[k]| / N
  initial phase    : arg(X[k])

Reconstruction (partial sum of M epicycles):
  z(t) = (1/N) · Σ_{k=0}^{M-1}  X[k] · exp(2πi·k·t / N)

The simulation sorts epicycles by decreasing amplitude so that the largest circles are drawn first. As M increases from 1 to N/2, the partial sum path converges to the original drawing. A slider lets you watch the approximation improve circle by circle.

Why this differs from the Fourier Series simulator

The platform's existing fourier-series simulator applies a 1D DFT to periodic waveforms (square, sawtooth, triangle) and visualises harmonic amplitudes on a bar chart. The fourier-epicycles simulator uses a 2D complex DFT on an arbitrary drawn path — both x(t) and y(t) are encoded together as a single complex signal. This means every frequency bin contributes a 2D circular motion (an epicycle), not just a vertical oscillation. The set of rotating circles is a visual proof that any continuous closed curve can be decomposed into circular motions.

Built-in presets and their spectra

The five presets were chosen to illustrate different spectral profiles:

• Circle — single non-zero bin at k = 1; one epicycle suffices for perfect reconstruction.
• Star (5-pointed) — dominant harmonics at k = 1, 4, 6, 9, 11, … (multiples of 5 ± 1); produces the characteristic pentagonal starbursts.
• Heart — asymmetric shape requires many harmonics; slow convergence reveals the cusp at the bottom.
• Figure-8 — symmetric; only odd harmonics present in the x-component, giving a clean spectrum.
• Lemniscate (∞) — parametrised as x = cos t / (1 + sin²t), y = sin t cos t / (1 + sin²t); dominated by k = 2 due to the double loop.

Rendering and performance

Computing a full DFT at N = 512 costs O(N²) multiplications. The simulator uses the Cooley-Tukey radix-2 FFT algorithm (padding to the next power of two) to reduce this to O(N log N):

// Iterative Cooley-Tukey in-place FFT
function fft(x):   // x: complex array, length must be power of 2
    N = x.length
    bit-reverse permutation of x
    for s = 1 to log2(N):
        m = 2^s
        w_m = exp(-2πi / m)
        for k = 0 to N-1 step m:
            w = 1
            for j = 0 to m/2 - 1:
                t = w * x[k + j + m/2]
                u = x[k + j]
                x[k + j]       = u + t
                x[k + j + m/2] = u - t
                w = w * w_m

The epicycle animation renders at 60 fps by recomputing only the tip positions of the rotating arms each frame (no FFT per frame — the DFT is computed once on mouse-up). The trailing path is drawn to an off-screen canvas and composited, keeping the main loop light.

🧲 Magnetic Pendulum — Chaos and Fractal Basins

Equations of motion

The pendulum bob is modelled as a point mass m suspended on a string of length L, constrained to move in the horizontal plane directly above the three magnets. The height z of the bob above the magnet plane is kept fixed at h (the rest height). The equations of motion in Cartesian coordinates (x, y) are:

Forces on the bob:
  Gravity restoring:  F_g = −(mg/L) · (x, y)   (small-angle approximation)
  Damping:            F_d = −b · (ẋ, ẏ)
  Magnetic (magnet i at position r_i = (x_i, y_i)):
    d_i = sqrt((x − x_i)² + (y − y_i)² + h²)
    F_m_i = μ · (r_i − r) / d_i³

Full ODE:
  ẍ = −(g/L)·x  −  b·ẋ  +  Σ_i  μ·(x_i − x) / d_i³
  ÿ = −(g/L)·y  −  b·ẏ  +  Σ_i  μ·(y_i − y) / d_i³

Default parameters: L = 0.5 m, h = 0.05 m, b = 0.2 (damping), μ = 1.0 (magnetic strength). The three magnets sit at angles 0°, 120°, 240° on a circle of radius 0.3 m.

RK4 integration

Each trajectory is integrated using the classical 4th-order Runge-Kutta method with a fixed time step Δt = 0.01 s. The state vector is s = (x, y, ẋ, ẏ):

k1 = f(t,       s)
k2 = f(t + h/2, s + h·k1/2)
k3 = f(t + h/2, s + h·k2/2)
k4 = f(t + h,   s + h·k3)
s_{n+1} = s_n + (h/6)·(k1 + 2k2 + 2k3 + k4)

The simulation halts when the bob speed drops below 0.001 m/s for 50 consecutive steps, recording which magnet is nearest as the "captured" magnet. For the basin map a grid of starting positions (up to 400×400 = 160,000 trajectories) is computed using Web Workers for chunked rendering so the UI stays responsive.

Fractal basin boundaries

The resulting basin map is a canonical example of a fractal. At every apparent boundary between two colour regions there is a third region nested infinitely deep — the boundaries have Hausdorff dimension strictly greater than 1. This is a direct consequence of sensitivity to initial conditions (positive Lyapunov exponents) near the separatrices. Zooming into the boundary reveals self-similar fine structure at every scale, making this one of the most visually accessible chaos demonstrations on the platform.

Two interactive controls are provided: a "magnetic strength" slider that weakens or strengthens the magnets relative to gravity, and a "damping" slider. At very low damping the bob wanders much longer before settling, producing sharper and deeper fractal structure. At very high damping it falls almost immediately into the nearest magnet, simplifying the basin boundaries toward straight Voronoi lines.

💧 Surface Wetting — From Young's Equation to the Lotus Effect

Young's equation (ideal flat surface)

When a liquid droplet rests on a flat solid surface in a vapour environment, three interfacial tensions act at the contact line: solid-vapour (γ_SV), solid-liquid (γ_SL), and liquid-vapour (γ_LV). Mechanical equilibrium along the surface requires:

Young's equation:
  γ_SV = γ_SL + γ_LV · cos(θ_Y)

  → cos(θ_Y) = (γ_SV − γ_SL) / γ_LV

Contact angle regimes:
  θ_Y = 0°         : complete wetting (liquid spreads indefinitely)
  0° < θ_Y < 90°   : hydrophilic (water wets the surface)
  90° < θ_Y < 150° : hydrophobic (water beads up)
  θ_Y > 150°        : superhydrophobic (lotus effect)

The animated droplet cross-section morphs continuously with the contact angle slider. The three interface arrows update direction and label to show the force balance.

Wenzel model (rough surface)

Real surfaces are never perfectly flat. Wenzel (1936) introduced a roughness factor r (ratio of actual surface area to projected area, always r ≥ 1) into the energy balance:

Wenzel equation:
  cos(θ_W) = r · cos(θ_Y)

  Consequence:
    If θ_Y < 90°  (hydrophilic): roughness makes it more wetting  → θ_W < θ_Y
    If θ_Y > 90°  (hydrophobic): roughness makes it less wetting  → θ_W > θ_Y

Roughness factor examples:
  Polished metal:  r ≈ 1.05
  Sand-blasted:    r ≈ 1.5 – 2.0
  Lotus leaf:      r ≈ 2.5 (microscale + nanoscale dual texture)

The surface roughness visualisation shows a magnified cross-section of the substrate with adjustable bump amplitude and spatial frequency. The droplet contact angle updates in real time as both θ_Y and r are varied.

Cassie-Baxter model and the lotus effect

When a surface has trapped air pockets (composite interface), the droplet rests partly on solid and partly on air. Cassie and Baxter (1944) derived the contact angle for this heterogeneous surface:

Cassie-Baxter equation:
  cos(θ_CB) = f₁ · cos(θ₁) + f₂ · cos(θ₂)

For solid-air composite (f₂ = air fraction, θ₂ = 180°):
  cos(θ_CB) = f_s · cos(θ_Y) − (1 − f_s)

  where f_s = fractional solid contact area

Lotus leaf: f_s ≈ 0.05 (only 5% solid contact)
  cos(θ_CB) ≈ 0.05 · cos(105°) − 0.95 ≈ −0.963
  θ_CB ≈ 164°  (superhydrophobic)

The simulator provides a toggle between Wenzel state (droplet penetrates roughness grooves) and Cassie-Baxter state (droplet sits on top of roughness peaks with air below). The Cassie-to-Wenzel transition (wetting transition) can be triggered by pressing the droplet down, modelled by crossing an energy barrier.

⚙️ Gear Train — Ratios, Torque and Compound Stages

Gear ratio fundamentals

Two meshing external spur gears with tooth counts N₁ and N₂ rotate at angular velocities ω₁ and ω₂. Because their pitch circles roll without slipping, the tangential velocity at the mesh point is equal:

Pitch circle radii: r₁ = N₁·m/2,  r₂ = N₂·m/2  (m = module = tooth size)

No-slip mesh condition: r₁·ω₁ = r₂·ω₂

  → gear ratio: ω₂/ω₁ = N₁/N₂  (magnitude)

Direction: external mesh reverses rotation → ω₂/ω₁ = −N₁/N₂
           internal mesh (ring gear) preserves rotation

Torque ratio (ideal, no friction):
  T₂/T₁ = N₂/N₁   (torque multiplied by same factor as speed is divided)

Idler gears

An idler gear (also called a parasite gear) is inserted between a driver and driven gear. It meshes with both but is not connected to any load or power source. Because it introduces two mesh reversals, the net direction of the driven gear is the same as the driver — as if no idler were present. The overall gear ratio between first and last gear is unaffected by the idler:

Driver (N₁) → Idler (N_i) → Driven (N₂):
  ω_i / ω₁ = −N₁/N_i       (first mesh, reversed)
  ω₂ / ω_i = −N_i/N₂       (second mesh, reversed again)
  ω₂ / ω₁  = N₁/N₂         (idler tooth count cancels)

Practical use: reverse output shaft direction without changing ratio

Compound gear trains

A compound gear train has two or more gears sharing a shaft (keyed together, rotating at the same angular velocity). The total ratio is the product of the individual stage ratios:

3-stage compound train:
  Stage 1: driver N₁ meshes with gear N₂ (on shaft A)
  Stage 2: gear N₃ (on shaft A) meshes with gear N₄ (on shaft B)
  Stage 3: gear N₅ (on shaft B) meshes with output N₆

  Total ratio: ω₆/ω₁ = (N₁·N₃·N₅) / (N₂·N₄·N₆)

Example: N₁=10, N₂=40, N₃=12, N₄=60, N₅=15, N₆=75
  Ratio = (10·12·15)/(40·60·75) = 1800/180000 = 1/100
  Input speed 1000 rpm → output 10 rpm, torque ×100

The animation renders each gear with a simplified involute tooth profile. Tooth engagement is tracked by monitoring the contact point position along the line of action, and the teeth animate smoothly in and out of mesh. An efficiency panel shows the approximate power loss from gear-mesh friction (typically 1–3% per stage for well-lubricated spur gears).

🎒 Knapsack — The 0/1 DP Table Exposed

Problem definition and brute-force complexity

Given n items, each with weight w_i and value v_i, and a knapsack of capacity W, find the subset of items with maximum total value subject to total weight ≤ W. Each item can be included at most once (the "0/1" constraint). The brute-force approach tries all 2ⁿ subsets — feasible only for n ≤ 20 or so.

Bottom-up dynamic programming

The DP recurrence fills an (n+1) × (W+1) table where K[i][w] = maximum value achievable using the first i items in a knapsack of capacity w:

Base case:
  K[0][w] = 0   for all w   (no items → zero value)

Recurrence (for i = 1…n, w = 0…W):
  if w_i > w:
      K[i][w] = K[i-1][w]            // item i too heavy: skip it
  else:
      K[i][w] = max(
          K[i-1][w],                 // don't take item i
          v_i + K[i-1][w - w_i]     // take item i
      )

Answer: K[n][W]
Complexity: O(n·W) time, O(n·W) space (reducible to O(W) with rolling array)

The simulation renders the complete table, colour-coded by value magnitude. As the animation fills each cell it highlights the two cells being compared (K[i-1][w] and K[i-1][w−w_i]) with a brief flash, making the recurrence structure immediately visible.

Backtracking to find selected items

Once the table is complete, the actual selected items are recovered by backtracking from K[n][W]:

i = n, w = W
while i > 0 and w > 0:
    if K[i][w] != K[i-1][w]:
        include item i           // cell changed → item i was taken
        w = w − w_i
    i = i − 1

Selected items are highlighted with a gold border in the item list. A side-by-side comparison panel runs the greedy algorithm (sort by value/weight ratio, greedily fill) and marks any discrepancy — for 0/1 knapsack, greedy frequently misses the optimum.

Why DP isn't truly polynomial

The O(nW) complexity is sometimes called "pseudo-polynomial" because W can be exponentially large in the number of bits used to encode it. If W = 2^k, the table has n·2^k cells — exponential in the input size. True polynomial-time algorithms for 0/1 knapsack are unknown (it is NP-hard). The simulation makes this concrete with a "bit-length" display: slide W from 100 to 100,000 and watch the table size (and render time) scale linearly with W but exponentially with log W.

🌳 Huffman Coding — Optimal Prefix-Free Compression

Building the tree with a min-heap

Given a symbol alphabet with known frequencies (or probabilities), Huffman's algorithm builds an optimal prefix-free binary code:

Input: symbols s_1…s_k with frequencies f_1…f_k

Algorithm:
1. Insert all symbols into a min-priority queue (min-heap), keyed by frequency.
2. While heap size > 1:
     a. Extract two nodes with minimum frequency: left (f_L), right (f_R)
     b. Create internal node with frequency f_L + f_R
     c. Insert internal node back into the heap
3. The remaining node is the root of the Huffman tree.

Code assignment:
  Traverse tree from root: going left appends '0', right appends '1'.
  Each leaf's code = the bit string on the path from root to that leaf.

The simulation animates each extraction and merge step. The growing tree is drawn with frequency labels on every node, and the two extracted nodes glow as they are merged.

Entropy and average code length

Shannon's source coding theorem states that no uniquely decodable code can achieve an average code length shorter than the entropy H of the source:

Shannon entropy:
  H = −Σ_{i=1}^{k}  p_i · log₂(p_i)    (bits per symbol)

Huffman average code length:
  L̄ = Σ_{i=1}^{k}  p_i · len(code_i)

Guarantee:  H ≤ L̄ < H + 1

For the optimal block code of length n symbols: L̄ → H as n → ∞

The live graph on the right panel plots H (horizontal reference line) and L̄ (bar) side by side. As you edit the input text, both values update and the gap between them is labelled in millibits. For natural English text H ≈ 4.1 bits/symbol and Huffman achieves L̄ ≈ 4.2 bits/symbol — within about 2% of the theoretical minimum using single-symbol codes.

Prefix-free property proof

Any code produced by reading tree paths is automatically prefix-free: no codeword is a prefix of another, because codewords only appear at leaves and no leaf is an ancestor of another leaf. This guarantees unambiguous decoding without any separator characters. The simulator includes a decoder panel — type in a bit string and it parses it symbol by symbol, highlighting the tree path taken for each decoded symbol.

Up Next

Wave 67 is already built and brings the library to a landmark milestone: the 600-simulation mark. The six additions span chaos theory (logistic map bifurcation diagram with Feigenbaum constant), optics (moiré patterns and Nyquist aliasing), computational geometry (Delaunay triangulation), transport modelling (Nagel-Schreckenberg traffic cellular automaton), algorithms (quadtree spatial partitioning), and orbital mechanics (tidal forces with spring/neap tide cycles). Read all about it in Devlog #88.