Devlog #84 – Wave 63: Stokes Flow, Sphere Packing, Hilbert Curve, Ulam Spiral, Max-Flow & Ice Halo

New Simulations

🌊 Stokes Flow — The Physics of Creeping Motion

The Stokes equations (Re → 0)

At very low Reynolds number Re = ρUL/μ ≪ 1, inertial forces are negligible compared to viscous forces and the Navier-Stokes equations reduce to the linear Stokes equations:

∇p = μ ∇²u       (momentum: pressure gradient balances viscous diffusion)
∇ · u = 0        (continuity: incompressible flow)

In 2D stream-function form (u = ∂ψ/∂y, v = −∂ψ/∂x):
  ∇⁴ψ = 0         (biharmonic equation)

Because the equations are linear and contain no time derivatives (time appears only through the boundary conditions), Stokes flow is quasi-static: the velocity field adjusts instantaneously to the boundary conditions. This has a profound consequence — reversibility.

Streamfunction for flow past a sphere

For uniform flow U past a sphere of radius R, the exact Stokes streamfunction in spherical coordinates (r, θ) is:

ψ(r, θ) = U · sin²θ · [r²/2  -  3R·r/4  +  R³/(4r)]

Velocity components:
  u_r = (1/r² sinθ) · ∂ψ/∂θ = U cosθ · [1 - 3R/(2r) + R³/(2r³)]
  u_θ = -(1/r sinθ) · ∂ψ/∂r = -U sinθ · [1 - 3R/(4r) - R³/(4r³)]

Drag force (Stokes drag):
  F_drag = 6π μ R U

The streamlines are symmetric fore and aft of the sphere — a property unique to Stokes flow. In contrast, at high Re the wake breaks the symmetry and drag increases dramatically.

The reversibility demonstration

The simulation's most striking feature is the reversibility demo. Coloured dye blobs are injected and advected forward by the Stokes velocity field. When the flow is reversed (U → −U), the dye particles retrace their exact paths and reconverge to their original positions — the system has no effective "mixing" because the linearity of Stokes flow means there is no chaotic stretching and folding. This is Taylor's famous corn-syrup experiment from 1966. The same demo with a Navier-Stokes flow at Re = 100 would show the blobs irreversibly smeared by turbulent mixing.

The comb theorem

A single horizontal dye line is divided into N segments, each forced through a separate parallel cylinder. After passing through the N cylinders, the N separate dye threads rejoin into a single horizontal line — identical to the original. This is a consequence of the linearity and reversibility of Stokes flow: each fluid particle's path through a Stokes-flow obstacle is uniquely determined by its entry position and reversible, so the exit positions are guaranteed to reassemble in order.

🔵 Sphere Packing — Kepler's Conjecture and Packing Fractions

2D packing fractions

In 2D, two regular (lattice) packings are known. The square packing places circle centres on a unit square grid; the hexagonal (triangular) packing staggers alternate rows by half a diameter:

Square packing:
  centres at (i, j) for i,j ∈ ℤ; radius r = 0.5
  area of circle / area of unit cell = π r² / 1² = π/4 ≈ 0.7854

Hexagonal packing:
  centres at (i + j/2, j·√3/2); radius r = 0.5
  area of unit cell = √3/2 (rhombus)
  packing fraction = π r² / (√3/2) = π/(2√3) ≈ 0.9069

Random sequential addition (RSA):
  Place disks one at a time at random positions if no overlap.
  Jammed packing fraction ≈ 0.547 (Renyi's parking constant)

The hexagonal packing is provably optimal in 2D — no arrangement of equal circles can achieve a higher density. This was proved by Thue in 1910. The 2D simulation shows live packing fraction and the characteristic hexagonal void pattern between circles.

3D crystal lattices and Kepler's conjecture

The 3D problem is far harder. The four classical lattice packings implemented in the simulation are:

Simple Cubic (SC):
  a = 2r;  packing fraction = π/6 ≈ 0.5236;  contact number Z = 6

Body-Centred Cubic (BCC):
  a = 4r/√3;  packing fraction = π√3/8 ≈ 0.6802;  Z = 8

Face-Centred Cubic (FCC):
  a = r√8;  packing fraction = π/(3√2) ≈ 0.7405;  Z = 12

Hexagonal Close-Packing (HCP):
  Same packing fraction as FCC = π/(3√2) ≈ 0.7405;  Z = 12
  Differs from FCC only in the stacking sequence (ABAB vs ABCABC)

Kepler conjectured in 1611 that π/(3√2) ≈ 74.05% is the maximum possible packing fraction for equal spheres. This was formally proved by Thomas Hales in 1998 using a computer-assisted proof of 250 pages plus 3 GB of code — one of the longest proofs in mathematics. Both FCC and HCP achieve this maximum.

3D projection and contact visualisation

The right panel shows an animated 3D projection of the selected lattice. Each sphere is rendered semi-transparent so the internal structure is visible. Contact bonds are drawn between touching sphere pairs (distance = 2r within floating-point tolerance). The contact number Z — number of nearest neighbours touching a given sphere — visually confirms the lattice type and equals the kissing number in the respective geometry.

🌀 Hilbert Curve — Space-Filling Curves and Locality

Construction by recursive substitution

The Hilbert curve is generated by a substitution grammar applied to a seed pattern. At order 1 it is a ∪-shape connecting 4 points. At order n, each of the 4 sub-squares from order n−1 is replaced by a rotated/reflected copy of the order n−1 curve, connected by bridging segments. The order-n curve visits 4ⁿ points and has total length (4ⁿ − 1) × (unit cell size), which diverges as n → ∞ while the curve remains continuous and contained in the unit square.

// d-to-(x,y) mapping (bit manipulation, after Skilling 2004):
function d2xy(n, d):
    rx, ry, s = 0, 0, 1
    x, y = 0, 0
    while s < n:
        rx = (d >> 1) & 1
        ry = (d ^ rx) & 1
        rotate(s, x, y, rx, ry)
        x += s * rx
        y += s * ry
        d >>= 2
        s <<= 1
    return x, y

// Inverse (x,y)-to-d:
function xy2d(n, x, y):
    rx, ry, s = 0, 0, n >> 1
    d = 0
    while s > 0:
        rx = (x & s) > 0 ? 1 : 0
        ry = (y & s) > 0 ? 1 : 0
        d += s * s * ((3 * rx) ^ ry)
        rotate(s, x, y, rx, ry)
        s >>= 1
    return d

Locality preservation

The key property that makes the Hilbert curve practically useful (in CPU cache management, database R-trees, GIS spatial indexing) is that points near each other on the 1D curve are also near each other in 2D space, and vice versa. This locality preservation can be quantified: for two 2D points with Euclidean distance d, their Hilbert indices differ by at most O(d). In contrast, the row-major ("raster scan") ordering has worst-case pairs at opposite ends of the index space that are just one row apart in 2D, and Z-order (Morton) is intermediate.

Peano and Morton curves

The Peano curve (1890, the first space-filling curve ever discovered) subdivides each square into 9 sub-squares in a serpentine pattern, visiting 9ⁿ points at order n. The Morton Z-order curve interleaves the bits of the x and y coordinates:

// Morton encoding (Z-order):
morton(x, y) = spread_bits(x) | (spread_bits(y) << 1)

spread_bits(v):  // insert a 0 bit between each bit of v
    v = (v | (v << 8)) & 0x00FF00FF
    v = (v | (v << 4)) & 0x0F0F0F0F
    v = (v | (v << 2)) & 0x33333333
    v = (v | (v << 1)) & 0x55555555
    return v

The hover tooltip shows both the forward mapping (curve distance d → grid position (x,y)) and the inverse, letting you see at a glance how locality is preserved or violated for any chosen point.

🔢 Ulam Spiral — Prime Patterns and Polynomial Ridges

The Ulam spiral construction

In 1963, Stanislaw Ulam, bored during a conference lecture, began writing integers in a square spiral and circling the primes. He noticed that the primes tended to cluster on diagonals — far more than would be expected by chance. The spiral assigns integer n to a grid position by winding outward from the centre:

n=1 → (0,0)
n=2 → (1,0)  n=3 → (1,1)  n=4 → (0,1)  n=5 → (-1,1)
n=6 → (-1,0) n=7 → (-1,-1) n=8 → (0,-1) n=9 → (1,-1) ...

General winding (move right 1, up 1, left 2, down 2, right 3, ...):
  segment lengths: 1,1,2,2,3,3,4,4,...
  directions: E,N,W,S,E,N,W,S,...

The simulation computes spiral positions for all N up to ~90,000 using a precomputed winding table — fast enough to render the full grid at 60 fps with interactive zoom and pan.

Why diagonals? Quadratic polynomials

The diagonals correspond to quadratic polynomials in n. Along one main diagonal: the values are f(k) = 4k² + 2k + 1 for k = 0, 1, 2, … Along the other: f(k) = 4k² − 2k + 1. Any quadratic that produces many primes will light up as a bright diagonal on the Ulam spiral. The most famous example is Euler's polynomial f(n) = n² + n + 41, which is prime for n = 0, 1, 2, …, 39 (40 consecutive values) — a record that stood for centuries. On the spiral this appears as an unmistakably bright diagonal ridge, toggle-able in the simulation.

The Sacks spiral

The Sacks spiral (Robert Sacks, 1994) maps integer n to polar coordinates by placing it at angle θ = 2π√n and radius r = √n — so that perfect squares lie on the positive x-axis and the spiral makes exactly one revolution per perfect square:

θ(n) = 2π √n
r(n) = √n
x(n) = r cosθ,  y(n) = r sinθ

Primes on the Sacks spiral form curved arcs rather than straight diagonals, but the clustering is even more visually striking. Each arc corresponds to a different quadratic residue class modulo small primes. The simulation renders both spirals side by side with synchronised zoom.

📦 Max-Flow — Edmonds-Karp and the Min-Cut Theorem

Flow networks and the max-flow problem

A flow network is a directed graph G = (V, E) where each edge (u, v) has a non-negative capacity c(u, v). Flow f(u, v) on each edge must satisfy: (1) capacity constraint f(u, v) ≤ c(u, v), and (2) flow conservation: flow into any internal node equals flow out. The max-flow problem asks for the maximum total flow from source s to sink t.

Max-flow = min-cut theorem (Ford-Fulkerson, 1956):
  The maximum s-t flow equals the minimum s-t cut capacity.

s-t cut: a partition (S, T) of V with s ∈ S, t ∈ T
Cut capacity = Σ c(u,v) for edges (u,v) with u ∈ S, v ∈ T

Edmonds-Karp algorithm

The Edmonds-Karp algorithm (1972) is Ford-Fulkerson with a specific augmenting path choice: always use the shortest path (fewest edges) from s to t in the residual graph, found by BFS. This choice guarantees polynomial time complexity O(VE²), whereas general Ford-Fulkerson can take O(E × max_flow) steps with unlucky path choices.

Edmonds-Karp:
1. Build residual graph G_f: for each edge (u,v) with capacity c and flow f,
       add forward edge residual capacity c - f
       add backward edge residual capacity f
2. Find shortest augmenting path s→t via BFS in G_f
3. Bottleneck Δ = min residual capacity along path
4. Augment: increase flow on each forward edge by Δ, decrease backward by Δ
5. Repeat until no augmenting path exists (BFS finds no path s→t)

Complexity: O(VE²) — at most VE/2 augmentations, each BFS in O(E)

Min-cut identification and animation

After the algorithm terminates, the min-cut is identified by BFS in the residual graph starting from s: all vertices reachable from s form the set S; all others form T. Every edge crossing from S to T at full capacity is a cut edge — drawn in red. The simulation shows each BFS augmenting step one-by-one: the chosen path is highlighted in yellow, the flow values on each edge update numerically, and the residual capacities animate. After convergence, the min-cut animation draws the partition boundary and labels it with its capacity (equal to the max-flow).

• Preset 1 — Classic textbook example; max-flow = min-cut = 23.
• Preset 2 — Diamond network; demonstrates a cut of size 3 separating source from sink.
• Preset 3 — Pipeline with a bottleneck edge; single cut edge is trivially the min-cut.
• Preset 4 — Bipartite matching formulation; max-flow gives maximum matching size.
• Preset 5 — Large random network; demonstrates performance with 20 nodes and 40 edges.

🌈 Ice Halo — Monte Carlo Optics of Atmospheric Ice Crystals

Physics of ice halos

Atmospheric ice halos are caused by refraction and reflection of sunlight through the faces of hexagonal ice prisms (Ih crystal structure). A hexagonal prism has two hexagonal basal faces and six rectangular prism faces. The key refraction angle depends on the prism angle α between two faces and the refractive index n of ice for visible light (n ≈ 1.31 at 550 nm, dispersive from ~1.307 at red to ~1.317 at violet).

Snell's law and minimum deviation

For a ray entering face 1 and exiting face 2 of a prism with wedge angle α, the deviation angle δ is:

At entry:   n_air · sin(θ_i) = n_ice · sin(θ_r)
At exit:    n_ice · sin(θ_r') = n_air · sin(θ_e)

Geometry:   θ_r + θ_r' = α   (angle between prism faces)
Total deviation: δ = θ_i + θ_e - α

Minimum deviation (δ_min) occurs when θ_i = θ_e (symmetric ray):
  sin((δ_min + α)/2) = n · sin(α/2)

For the 22° halo: α = 60° (alternate prism faces), n = 1.31
  sin((δ_min + 60°)/2) = 1.31 · sin(30°) = 0.655
  (δ_min + 60°)/2 = arcsin(0.655) ≈ 40.9°
  δ_min ≈ 21.8° ≈ 22°

The minimum deviation angle defines the inner edge of the halo ring: no ray can exit with deviation less than δ_min for that prism geometry. Rays pile up near δ_min (stationary phase), producing a bright ring at that angular radius.

Monte Carlo ray tracing

Rather than integrating the deviation distribution analytically, the simulation fires N = 100,000 rays per frame, each through a randomly oriented (or constrained) ice crystal:

for each ray:
    1. Sample crystal orientation:
         random → all 22° halo rays
         plate (tilt σ ≈ 1°) → sun dog / parhelion
         column (tilt σ ≈ 1°) → circumzenithal arc
    2. Compute entry face hit: intersect ray with prism geometry
    3. Apply Snell's law at entry → refracted direction inside crystal
    4. Propagate to exit face; apply Snell's law again
    5. Compute (azimuth, elevation) of exit ray
    6. Accumulate in 2D angular histogram (sky projection)
    7. Apply wavelength dispersion: repeat for λ = 400, 550, 700 nm

The resulting histogram is normalised and rendered as a sky image with an angular coordinate system. The 22° halo appears as a diffuse bright ring at 22° angular radius from the simulated sun, brighter on the inside edge (minimum deviation). Red is on the inside of the ring (smaller n → larger δ_min), blue on the outside — the opposite of a rainbow.

Sun dogs and the circumzenithal arc

Sun dogs (parhelia) are caused by plate crystals that settle with their c-axis vertical due to aerodynamic drag. Rays enter through a vertical prism face and exit through another vertical face at 60°, but now only crystals with the correct azimuthal angle contribute, concentrating light at ±22° from the sun along the horizon:

Parhelion condition: light enters/exits vertical prism faces
  Crystal tilt from horizontal: σ ≈ 0°–2° (plate settling)
  Angular position: at solar altitude h, sun dog is displaced
    Δaz = arctan(tan(22°)/cos(h))
  At h = 0°: Δaz = 22°;  at h = 30°: Δaz ≈ 26°

The circumzenithal arc (CZA) forms when light enters through the top basal face and exits through a vertical prism face of plate crystals. The geometry gives a curved arc 46° above the sun, with red on the outside (bottom) and violet on the inside — producing the most intensely coloured arc in atmospheric optics, often more saturated than a rainbow.

Up Next

Wave 64 is in development with candidates spanning drag coefficient aerodynamics, geodesic dome geometry, stable matching (Gale-Shapley), hologram interference patterns, the Louvain community detection algorithm, and grain structure / Voronoi crystallography. The library is on track to exceed 600 simulations by summer.