Devlog #86 – Wave 65: Capillary Action, A*, Parametric Curves, Vortex Beam, Colloids & Internet Routing

New Simulations

💧 Capillary Action — Jurin's Law and the Meniscus

Surface tension and contact angle

Surface tension γ (N/m) arises from the asymmetric intermolecular forces at a liquid-air interface: molecules at the surface have fewer neighbours than bulk molecules and are in a higher energy state. When a liquid contacts a solid, three interface energies compete: liquid-air (γ_LA), solid-liquid (γ_SL), and solid-air (γ_SA). Their balance determines the contact angle θ via Young's equation:

Young's equation (force balance at contact line):
  γ_SA = γ_SL + γ_LA · cos θ
  cos θ = (γ_SA − γ_SL) / γ_LA

Contact angle classes:
  θ = 0°:    perfect wetting (liquid spreads completely)
  0 < θ < 90°: wetting (concave meniscus, capillary rise)
  θ = 90°:   neutral
  90° < θ < 180°: non-wetting (convex meniscus, capillary depression)
  θ = 180°:  perfect non-wetting (Cassie-Baxter superhydrophobic state)

Water on glass:  θ ≈ 20°  (wetting)
Mercury on glass: θ ≈ 140° (non-wetting)

Jurin's law: capillary rise

In a tube of radius r, the meniscus creates a pressure difference across the curved surface (Young-Laplace equation). For a spherical meniscus in a circular tube, this pressure difference supports a column of liquid of height h:

Young-Laplace pressure across curved interface:
  ΔP = 2γ / R_c
  where R_c = r / cos θ   (radius of curvature of meniscus in tube)

Pressure balance (capillary pressure = hydrostatic pressure):
  2γ cos θ / r = ρ g h

Jurin's law:
  h = 2γ cos θ / (ρ g r)

Example — water in glass (r = 0.5 mm):
  γ = 0.073 N/m, θ = 20°, ρ = 1000 kg/m³, g = 9.81 m/s²
  h = 2 × 0.073 × cos(20°) / (1000 × 9.81 × 0.0005) ≈ 28 mm

The simulation renders side-by-side tubes of different radii with animated liquid rise. The meniscus shape is drawn as a circular arc with curvature 1/R_c. A live plot of h vs 1/r confirms the linear relationship predicted by Jurin's law. For mercury, the contact angle θ > 90° gives cosθ < 0, so h is negative — mercury is depressed below the reservoir level.

Temperature dependence of surface tension

Surface tension decreases approximately linearly with temperature for most liquids. For water, γ(T) can be approximated by Eötvös rule and is tabulated from 0°C (0.0756 N/m) to 100°C (0.0589 N/m). The simulation's temperature slider adjusts γ using real water data and recomputes the rise height in all tubes simultaneously, demonstrating that heated water climbs lower — relevant to capillary phenomena in boiling and microfluidics.

🗺️ A* Pathfinding — Heuristics, Optimality, and the Open Set

The A* algorithm

A* (Hart, Nilsson, Raphael, 1968) searches a weighted graph for the shortest path from start s to goal g by evaluating each node n with a cost estimate f(n) = g(n) + h(n):

f(n) = g(n) + h(n)

  g(n) = exact cost of shortest path found so far from s to n
  h(n) = heuristic estimate of cost from n to goal g

A* maintains two sets:
  Open set:   nodes discovered but not yet expanded (min-heap by f value)
  Closed set: nodes already expanded

Main loop:
  while open set is not empty:
    n = node in open set with lowest f(n)   // O(log N) with binary heap
    if n == goal: return path
    add n to closed set
    for each neighbour m of n:
      tentative_g = g(n) + edge_cost(n, m)
      if m in closed set and tentative_g ≥ g(m): continue
      if tentative_g < g(m):
        g(m) = tentative_g
        f(m) = g(m) + h(m)
        parent(m) = n
        add m to open set (or decrease key if already present)

Heuristic admissibility and optimality

A* is guaranteed to find the optimal path if and only if the heuristic h(n) is admissible: it never overestimates the true cost to reach the goal. The three heuristics implemented are:

Grid coordinates (x, y) for a cell, goal at (gx, gy):

Manhattan distance (L1 norm):
  h(n) = |x − gx| + |y − gy|
  Admissible for 4-connected grids (no diagonal movement)
  Inadmissible for 8-connected grids (underestimates diagonal shortcuts)

Euclidean distance (L2 norm):
  h(n) = √((x − gx)² + (y − gy)²)
  Admissible for any movement; tends to expand fewer nodes than Manhattan
  on 8-connected grids but requires float comparisons

Diagonal (Chebyshev) distance:
  h(n) = max(|x − gx|, |y − gy|)
  Admissible for 8-connected grids with unit diagonal cost
  Exact for uniform-cost diagonal movement; expands fewest nodes

Comparison with Dijkstra and greedy best-first

The algorithm comparison panel runs all three algorithms simultaneously on the same maze and displays node expansion counts side by side:

• Dijkstra (h = 0): Expands in concentric rings; guaranteed optimal; most nodes expanded.
• Greedy best-first (f = h only, ignores g): Rushes toward the goal; fewest nodes expanded; not optimal — may find a longer path.
• A* with Manhattan/Euclidean/Diagonal: Optimal and typically expands far fewer nodes than Dijkstra; Diagonal heuristic is closest to perfect on 8-connected grids.

The binary heap implementation keeps each iteration O(log N) where N is the open-set size, making A* practical on grids up to 500×500 (250,000 nodes) at interactive speeds. The node expansion counter updates live on each frame of the step-by-step animation.

✏️ Parametric Curves — Bézier, B-Splines, and NURBS

Bézier curves and de Casteljau construction

A degree-n Bézier curve is defined by n+1 control points P₀, …, P_n and parameterised by t ∈ [0, 1]. The de Casteljau algorithm evaluates any point on the curve by repeated linear interpolation:

de Casteljau algorithm:
  P_i^(0) = Pᵢ             (initialise with control points)
  P_i^(r) = (1−t) · P_i^(r−1) + t · P_{i+1}^(r−1)   for r = 1, …, n

  Curve point at parameter t: B(t) = P_0^(n)

Bernstein basis representation (equivalent):
  B(t) = Σᵢ₌₀ⁿ  C(n,i) · tⁱ · (1−t)^(n−i) · Pᵢ
  where C(n,i) = n! / (i! (n−i)!)  (binomial coefficient)

Degree 3 (cubic) Bézier — most common:
  B(t) = (1−t)³P₀ + 3t(1−t)²P₁ + 3t²(1−t)P₂ + t³P₃

The simulation animates de Casteljau construction frame by frame as t sweeps 0→1: each intermediate interpolation segment is drawn in a progressively lighter colour, and the final point traces the curve. Dragging control points updates the animation in real time.

B-splines and the Cox-de Boor recurrence

B-splines generalise Bézier curves with local control: moving one control point only affects the curve in a local neighbourhood rather than the entire curve. They are defined by a knot vector T = {t₀, …, t_m} and basis functions N_i,k(t) computed by Cox-de Boor:

Cox-de Boor recurrence:
  N_{i,0}(t) = 1  if tᵢ ≤ t < t_{i+1},  else 0
  N_{i,k}(t) = (t − tᵢ)/(t_{i+k} − tᵢ) · N_{i,k−1}(t)
             + (t_{i+k+1} − t)/(t_{i+k+1} − t_{i+1}) · N_{i+1,k−1}(t)

B-spline curve of degree k with n+1 control points:
  C(t) = Σᵢ₌₀ⁿ  N_{i,k}(t) · Pᵢ

Properties:
  Local support: N_{i,k}(t) ≠ 0 only on [tᵢ, t_{i+k+1}]
  Partition of unity: Σᵢ N_{i,k}(t) = 1 for all t
  C^(k−1) continuity everywhere (C^(k−2) at repeated knots)

The knot vector editor lets you drag individual knot values or insert/remove knots. Repeated knots reduce continuity locally — a knot repeated k times produces a corner at that parameter value, useful for cusps and sharp features. The simulation colour-codes each basis function N_i,k(t) individually and plots all of them in a strip below the main canvas.

NURBS and the exact circle

Non-Uniform Rational B-Splines (NURBS) add a weight w_i to each control point, enabling exact representation of conics — including circles and arcs — which polynomial splines can only approximate:

NURBS curve:
  C(t) = Σᵢ N_{i,k}(t) · wᵢ · Pᵢ
         ─────────────────────────
         Σᵢ N_{i,k}(t) · wᵢ

Exact quarter-circle with 3 control points (degree 2):
  P₀ = (1, 0),  w₀ = 1
  P₁ = (1, 1),  w₁ = √2/2 ≈ 0.7071   (weight of intermediate point)
  P₂ = (0, 1),  w₂ = 1

  Knot vector: [0, 0, 0, 1, 1, 1]

Full circle: 9 control points (three quarter-arcs joined C¹ continuous)
             Weights alternate between 1 and √2/2

The NURBS panel shows the control polygon with per-point weight sliders. Reducing a weight below 1 pulls the curve away from that control point; increasing it above 1 draws the curve toward it. The exact circle is the default NURBS example — changing the central weight from √2/2 to 1 deforms it into a parabola, demonstrating the role of weights.

🌀 Vortex Beam — Laguerre-Gaussian Modes and Orbital Angular Momentum

Laguerre-Gaussian beam modes

Laguerre-Gaussian (LG) modes are exact solutions to the paraxial wave equation in cylindrical coordinates (ρ, φ, z). They form an orthogonal basis for beams with rotational symmetry. An LG beam of radial index p and azimuthal index ℓ has electric field amplitude:

LG_{p,ℓ}(ρ, φ, z):

  u(ρ, φ, z) = C_{p,ℓ} · (ρ√2/w(z))^|ℓ| · L_p^|ℓ|(2ρ²/w²(z))
              · exp(−ρ²/w²(z))                     ← Gaussian envelope
              · exp(i ℓ φ)                          ← helical phase
              · exp(−i k ρ²/(2R(z)))               ← wavefront curvature
              · exp(i (2p+|ℓ|+1) ζ(z))            ← Gouy phase

where:
  w(z)  = beam waist at position z: w(z) = w₀√(1 + (z/z_R)²)
  R(z)  = radius of curvature: R(z) = z(1 + (z_R/z)²)
  ζ(z)  = Gouy phase: arctan(z/z_R)
  z_R   = Rayleigh range: πw₀²/λ
  L_p^|ℓ| = generalised Laguerre polynomial

Helical phase and orbital angular momentum

The azimuthal phase factor exp(iℓφ) winds the phase through 2πℓ radians around one full rotation in φ. This helical phase front means the beam carries orbital angular momentum (OAM) of ℓℏ per photon — distinct from the spin angular momentum ±ℏ of circularly polarised light:

OAM per photon: L_z = ℓ ħ

Intensity profile: I ∝ |u|² ∝ ρ^(2|ℓ|) · [L_p^|ℓ|(…)]² · exp(−2ρ²/w²)

For p = 0 (no radial nodes):
  Single bright ring ("doughnut") at ρ_max = w(z) · √(|ℓ|/2)
  Dark core (phase singularity) at ρ = 0

For p > 0:
  p concentric dark rings inside the outer ring
  Total of p+1 bright rings

Topological charge ℓ:
  |ℓ| = number of 2π phase windings per azimuthal circuit
  Sign determines handedness (clockwise vs anti-clockwise phase ramp)

The simulation renders the 2D intensity profile as a heat map and overlays the phase map as a colour wheel at any chosen z-plane. An interference panel shows the characteristic spiral (single-slit) or fork (double-slit) patterns used in the laboratory to measure the topological charge ℓ. When an LG beam interferes with a plane wave, the resulting pattern has |ℓ| extra or missing fringes at the fork point.

🧪 Colloids — DLVO Theory and Colloidal Stability

Brownian motion and the Stokes-Einstein equation

A colloidal particle (radius r, 1 nm–1 μm) suspended in a fluid undergoes Brownian motion because thermal fluctuations continuously kick it. The diffusion coefficient D relates the mean-squared displacement to time and is given by the Stokes-Einstein equation:

Stokes-Einstein equation:
  D = k_B T / (6π η r)

where:
  k_B = 1.38 × 10⁻²³ J/K  (Boltzmann constant)
  T   = temperature (K)
  η   = dynamic viscosity of solvent (Pa·s)
  r   = particle radius (m)

Einstein relation (diffusion–mobility):
  D = μ k_B T   where μ = mobility = 1/(6π η r)

Mean-squared displacement:
  ⟨r²⟩ = 6 D t   (3D Brownian motion)
  ⟨r²⟩ = 4 D t   (2D)

Example — 100 nm particle in water at 25°C:
  D = (1.38×10⁻²³ × 298) / (6π × 8.9×10⁻⁴ × 50×10⁻⁹) ≈ 4.9 × 10⁻¹² m²/s

DLVO theory: stability and flocculation

DLVO theory (Derjaguin, Landau, Verwey, Overbeek) describes the total interaction potential between two charged colloidal particles as the sum of electrostatic repulsion and van der Waals attraction:

Total DLVO potential:
  V_total(h) = V_EDL(h) + V_vdW(h)

Electrostatic double-layer repulsion (linearised Poisson-Boltzmann):
  V_EDL(h) = 64π ε ε₀ r (k_B T / ze)² tanh²(zeψ₀/4k_BT) · exp(−κh)

  κ = Debye-Hückel screening length:
  κ⁻¹ = √(ε ε₀ k_B T / (2 N_A e² I))
  I = ionic strength = ½ Σ cᵢ zᵢ²

van der Waals attraction (Hamaker):
  V_vdW(h) = −A_H / (12π h²) · f(r, h)
  A_H = Hamaker constant (typical value ~10⁻²⁰ J for oxide/water/oxide)

Stability barrier:
  If V_max > ~15 k_B T: stable colloid (kinetically inhibited flocculation)
  If V_max < ~5 k_B T:  rapid coagulation (Smoluchowski fast aggregation)
  Adding salt compresses Debye length → lowers barrier → flocculation

The ionic strength slider compresses the Debye length κ⁻¹ from ~100 nm (low salt) to ~1 nm (physiological), reducing the energy barrier V_max and triggering flocculation. The simulation visualises this as a live V(h) curve with the barrier height annotated in units of k_BT. When the barrier drops below ~5k_BT, particles begin aggregating in the canvas and the aggregate size distribution histogram updates in real time.

Sedimentation

Particles larger than ~1 μm settle under gravity. The sedimentation velocity for a sphere in the Stokes regime is:

Stokes sedimentation velocity:
  v_s = 2r²(ρ_p − ρ_f) g / (9η)

Péclet number (advection vs diffusion):
  Pe = v_s · r / D = v_s · 6π η r² / (k_B T)

Pe < 1: diffusion dominates → stable suspension
Pe > 1: sedimentation dominates → gravitational settling

Example — 1 μm silica sphere in water:
  ρ_p = 2200 kg/m³, ρ_f = 1000 kg/m³
  v_s = 2 × (500×10⁻⁹)² × 1200 × 9.81 / (9 × 8.9×10⁻⁴) ≈ 0.73 μm/s
  D = 4.9×10⁻¹³ m²/s  →  Pe ≈ 1.5 (borderline)

The simulation colours particles by their current aggregate cluster size (singleton = blue, dimer = green, …, large cluster = red) and overlays settling velocity vectors for particles where Pe > 1.

🌐 Internet Routing — Distance-Vector, Link-State, and BGP

Bellman-Ford distance-vector routing

Distance-vector protocols (e.g., RIP) have each router periodically broadcast its entire routing table to direct neighbours. Each router updates its own table using the Bellman-Ford recurrence:

Bellman-Ford recurrence:
  d_x(y) = min over all neighbours v of { c(x,v) + d_v(y) }

  d_x(y) = estimated cost from router x to destination y
  c(x,v) = cost of direct link x → v
  d_v(y) = neighbour v's advertised cost to y

Update rule (when x receives update from neighbour v):
  for each destination y:
    if c(x,v) + d_v(y) < d_x(y):
      d_x(y) = c(x,v) + d_v(y)
      next_hop_x(y) = v    // update routing table

Convergence: O(diameter) update rounds for correct topology

The count-to-infinity problem

When a link fails, distance-vector protocols can fall into a pathological loop where two routers alternately increase their advertised cost to a destination, bouncing between each other until the distance "counts to infinity" (typically capped at 16 in RIP):

Topology: A — B — C  (B-C link fails; only A-B link remains)

Before failure: d_B(C) = 1, d_A(C) = 2
After failure:
  B notices B-C link is down; d_B(C) = ∞
  But B's table still shows d_A(C) = 2 (stale)
  B updates: d_B(C) = 1 + d_A(C) = 3  (wrong! via A, but A routes via B)
  A updates: d_A(C) = 1 + d_B(C) = 4
  B updates: d_B(C) = 5 … and so on until ∞

Fix — Split Horizon:
  Router A never advertises to B a route it learned from B
  → A's advertisement for C (learned via B) is suppressed to B
  → B correctly learns d_B(C) = ∞ after link failure

Fix — Poison Reverse:
  More aggressive: A advertises d_A(C) = ∞ explicitly back to B
  Eliminates count-to-infinity for two-node loops immediately

The simulation lets you trigger a link failure by clicking any edge and then step through the Bellman-Ford update rounds. With split horizon off, the count-to-infinity loop is clearly visible in the router tables. Enabling split horizon shows the same failure resolved in a single round.

Dijkstra link-state (OSPF) and reconvergence

Link-state protocols (OSPF, IS-IS) take a fundamentally different approach: every router floods the entire network with its local link-state advertisement (LSA), so all routers build an identical topology map. Each router then runs Dijkstra locally to compute shortest paths:

Link-State Flooding:
  Each router generates LSA: {router ID, [neighbour, cost, ...]}
  LSA is flooded to all routers in the area (sequence numbers prevent loops)
  All routers build identical Link-State Database (LSDB)

Dijkstra on the LSDB (per router):
  Identical to Dijkstra with the router itself as source
  Result: full shortest-path tree to all destinations
  Installed in routing table as forwarding entries

Reconvergence after link failure:
  Failed router/link generates new LSA with infinite cost
  Floods network → all routers update LSDB → re-run Dijkstra
  Convergence in O(E log V) per router vs O(diameter) rounds for BF

BGP (Border Gateway Protocol):
  Path-vector protocol: advertises full AS path, not just distance
  Eliminates routing loops: reject updates containing own AS in path
  Policy-based: operators control which routes are accepted/preferred

The reconvergence animation shows LSA flood propagation as a wave across the network diagram, then highlights each router re-running Dijkstra as the flood reaches it. Router tables update visually in the panel on the right. The BGP tab demonstrates the path-vector mechanism with a simplified three-AS topology, showing how the AS-path attribute prevents the count-to-infinity problem that plagues distance-vector protocols.

Up Next

Wave 66 is in preparation with six further simulations: Fourier epicycles animating any closed curve as a sum of rotating circles, the magnetic pendulum with fractal basins of attraction, surface wetting and contact angle hysteresis, gear train kinematics, the 0/1 knapsack problem solved by dynamic programming, and Huffman coding with entropy bounds. The library approaches 600 simulations.