Glass Leaf Fall

A WebGL2 research demo combining maple leaf morphogenesis, N-BK7 glass optics, quasi-steady falling of a thin plate, incompressible turbulence, and finite Reynolds number drag on glass particles, all within a single real-time rendering engine

The capability of AI is remarkable. I fed it the papers directly and it reproduced the results almost exactly.

This work does not claim to have "computed a real glass maple leaf" faithfully. It uses realistic size and mass for the aerodynamics, but the optical thickness is intentionally exaggerated to produce visible refraction; the leaf venation is generated from a plant morphology model, and the optical surfaces treat those veins as glass ridges.

1. Model Scope and Attribution

GLSL is owned exclusively by shaders/.
In practice, when working with AI there are cases where the two get mixed together, so keeping them explicitly separated turned out to be important. For the leaf model, the outline was captured via vision OCR and handed over as-is, then implemented in .mjs form.

No shader strings or low-quality fallback shaders are placed inside JavaScript. js/shaders.js fetches explicit file URLs and synthesizes #include "..." exactly once per program. If a file or extension is missing, the system fails at startup rather than silently falling back.

index.html
style.css
js/
  main.js boot, frame loop, system wiring
  shaders.js explicit GLSL file loader and include resolver
  glass.js N-BK7 Sellmeier equation and per-wavelength refractive index
  dynamics.js rigid-body motion of a thin elliptical disk with RK4
  air.js divergence-free turbulent velocity field
  world.js leaf, fragment, and particle populations with measurement statistics
  renderer.js WebGL2 passes, instancing, GPU audit
  controls.js UI and audit table
shaders/
  leaf-outline.glsl
  leaf-venation.glsl generated file
  leaf-thickness.glsl
  background.frag
  leaf.frag
  particle.frag
  audit.frag
  post.frag
tools/
  leaf-model.mjs outline, tooth/serration, and venation generator

Because ES modules and shader fetch calls are used, the project must be opened over HTTP rather than file://.

Runtime asset policy

• No images, external textures, audio, or CDN JavaScript are used.

• The background, bokeh, leaf surfaces, and glass dust are all generated procedurally by local GLSL.

• The only external automatic request is Google Fonts for the UI, and it has no involvement in the rendering output.

• Every .vert, .frag, and .glsl file is registered in the explicit new URL(..., import.meta.url) manifest in js/shaders.js. This is a requirement for Makonea Lab to statically identify project files inside an iframe CSP.

In practice, applying this runtime asset policy successfully isolates the project within a sandboxed block. Personally, what an implementation looks like ultimately changes depending on the state of whatever you are trying to deploy.

2. Reference Constants and Default State

The values below are not rewritten file by file with similar-looking numbers. Air density and gravity are exported from dynamics.js, and it is verified that both leaves and particles import the same values.

Optical thickness and aerodynamic thickness are not the same value.

The base leaf falls as a plate with a total aerodynamic thickness of 200 µm, but the optical thickness field on screen occupies approximately 0.072–0.174 in leaf coordinates. Converted directly to a half-width of 70 mm, that corresponds to glass several millimeters thick. A 200 µm plate can move like a falling leaf but produces almost no strong caustics on screen, while glass several millimeters thick refracts beautifully but does not fall like a leaf.

3. Maple Leaf: Building Outline and Venation as a Single Structure

3.1 Lobe Profile

Each lobe is a function with axis c, length l, and angular half-width w

To simultaneously produce broad shoulders and a sharp tip, the code takes the maximum of two profiles per lobe. The actual code computes the expression directly.

float lobe(float th, float c, float len, float w) {
    return len * pow(max(0.0, 1.0 - abs(th - c) / w), LOBE_P);
}

float leafEnvelope(float th) {
    float r = 0.0;
    for (int i = 0; i < 5; i++) {
        r = max(r, lobe(th, LOBE_AXIS[i], LOBE_BL[i], LOBE_BW[i]));
        r = max(r, lobe(th, LOBE_AXIS[i], LOBE_SL[i], LOBE_SW[i]));
    }
    return r;
}

A profile whose derivative at the axis is zero, like a raised cosine, makes the tip rounded and wide. The current linear-taper power function leaves a corner at the axis, preserving the sharp tip characteristic of maple leaves.

3.2 Teeth Are Spaced by Arc Length, Not by Angular Interval

The arc length of the polar outline r(θ) is as follows.

The generator divides this s into equal intervals to place tooth vertices, which is a surprisingly elegant implementation. Dividing by angle would cause teeth to spread apart at lobe tips and cluster at indentations; the radius between two vertices is

indented inward. Because q=TOOTH_SHARP>1, V-shaped indentations form between vertices, and near the lobe axis d is sent to 0 to prevent the lobe tip from becoming a fan of teeth.

float u = (th - t0) / (t1 - t0);
float base = mix(TOOTH_R[lo], TOOTH_R[lo + 1], u);
float notch = pow(sin(PI * u), TOOTH_SHARP);
r = min(base * (1.0 - dip * notch), env);
r = max(r, floorR);

3.3 Space Colonization Venation

First, five actinodromous primary veins are placed along the lobe axes. Then, auxin sources are scattered inside the leaf blade and just inward of the serrations. Each source attracts the nearest node within an influence radius of R_I=0.13, and a node is consumed when a new node reaches within R_K=0.045; the growth step size is D=0.034 in leaf coordinates.

Because the serration vertices themselves act as sources, secondary veins and minor veins grow toward the teeth. Rather than overlaying the outline and venation as separate decorations, the craspedodromous condition is built directly into the generation process.

3.4 Murray Radius Law

Vein radii at branching points are computed back-to-front using the following conservation equation:

for (const i of order) {
    if (children[i].length === 0) { w[i] = 1; continue; }
    let s = 0;
    for (const c of children[i]) s += Math.pow(w[c], 3);
    w[i] = Math.pow(s, 1 / 3);
}

The w values up to this point are normalized radii of the biological veins, but the Gaussian ridge drawn on screen is given a minimum radius of ε=0.065 to prevent thin veins from vanishing at sub-pixel scales.

float radius = mix(VEIN_MIN_RADIUS, 1.0, w);
float widthSq = uVeinW * radius * radius;
float ridge = w * exp(-d * d / max(widthSq, 1e-7));

For a semicircular cross-section ridge, both height and half-width are proportional to the radius, so the height exponent is set to 1. An earlier implementation used w^2.35, but that made the finest branches far flatter than the Murray radius would require (and it genuinely did not look good). If a separate aesthetic height law is desired, it should be specified as an independent parameter rather than mixed in under the Murray label.

Multiple line segments are combined using a 6th-order Lp blend instead of max or a plain sum.

float ridge2 = ridge * ridge;
union6 += ridge2 * ridge2 * ridge2;
return pow(max(union6, 1e-18), 1.0 / 6.0);

max introduces curvature creases where the dominant segment changes, while a plain sum double-counts shared segments at branch points and produces bright beads. With an L6 blend, two segments of equal height meeting at a point produce a combined height of only 2^(1/6)=1.122 times the individual height.

3.5 Generation Results

The header of leaf-venation.glsl records values measured directly by the generator:

Each leaf is assigned one of three venation meshes and one of two y-axis mirror options, yielding six possible venation patterns. The outline is shared among leaves of the same species, while only the venation varies from leaf to leaf.

4. N-BK7 Dispersion: From Equation to Running Code

4.1 Sellmeier Equation

Insert wavelength λ in µm

const B = [1.03961212, 0.231792344, 1.01046945];
const C = [0.00600069867, 0.0200179144, 103.560653];

export function nBK7(um) {
    const l2 = um * um;
    let s = 1;
    for (let i = 0; i < 3; i++) s += (B[i] * l2) / (l2 - C[i]);
    return Math.sqrt(s);
}

The values derived from the coefficients are as follows.

The renderer samples three representative RGB wavelengths: 650, 550, and 450 nm.

Shorter-wavelength blue light bends the most while red bends the least; the particle shader receives the same scale factors as a uniform and applies them identically.

vec3 bend = (uIndex - 1.0) / (uIndex.g - 1.0);
vec3 bent = vec3(
    texture(uBg, vScreen + px * bend.r).r,
    texture(uBg, vScreen + px * bend.g).g,
    texture(uBg, vScreen + px * bend.b).b);

The earlier particle values 1.15 / 1.00 / 0.84 not only exaggerated the spread of dispersion but also reversed the order, displacing red the most. The current implementation has the leaf and particles share the same N-BK7 dispersion.

5. Thickness Field and Thin-Element Optics

5.1 Separating Outline from Thickness

The optical body is the sum of an ellipsoidal cap in Cartesian coordinates and a leaf-vein ridge.

float thickness(vec2 p, int variant) {
    float d = domeProfile(p);
    return uThick * (DOME_FLOOR + (1.0 - DOME_FLOOR) * d)
         + uVein * veinField(p, variant) * (0.25 + 0.75 * d);
}

The thickness field is not clipped to zero at the leaf outline, because if a finite-difference stencil crosses that step it reads the infinite curvature of the clipped boundary rather than the true surface curvature. The outline is handled at the end via coverage and grazing-angle Fresnel (honestly, I still don't know why I did it this way).

5.2 Thin-Element Mapping

Let d denote the propagation distance to the background plane; a thin thickness field t displaces a ray by

The Jacobian of this mapping is

and in the single-lobe approximation the brightness scale factor is

where regions passing through det J=0 appear as folds and caustic lines.

5.3 Gradient and Hessian

The gradient and Hessian are computed together from a total of 9 samples: one center point, four axial points, and four diagonal points.

float txx = (px - 2.0 * t0 + mx) / (e * e);
float tyy = (py - 2.0 * t0 + my) / (e * e);
float txy = (pp - pm - mp + mm) / (4.0 * e * e);

j.grad = vec2(px - mx, py - my) / (2.0 * e);
j.detJ = (1.0 + k * txx) * (1.0 + k * tyy)
       - (k * txy) * (k * txy);

An earlier implementation multiplied again by e², which the raw second-order finite difference already incorporated, reducing the Hessian contribution by a factor of e²=3.6×10⁻⁵. Conversely, an attempt that used the outline itself as the dome radius fed all high-frequency curvature from the teeth and indentations directly into the thickness, folding a large area; the current version corrects both the Hessian units and the smoothness of the thickness field independently.

5.4 Caustics for Screen Display

1/|det J| diverges at folds, and the current renderer does not gather and integrate multiple preimages. Therefore the raw determinant is preserved in the audit table, while the screen uses a log-compressed single-lobe gain.

float logGain = -log(max(abs(detJ), 0.08));
return exp(0.58 * clamp(logGain, -0.72, 0.72));

The denominator is not a value recomputed similarly on the CPU; instead, audit.frag renders the same thicknessJet() to a 128² RGBA32F target and measures the mean interior value of the leaf via readPixels (to keep computation faster).

5.5 Fresnel, Beer–Lambert, and 3D Orientation

The surface normal is constructed from the gradient of the front face at z=t/2.

This normal is rotated by the integrated leaf tilt and screen spin, rather than using the 2D trick of flattening only the silhouette while keeping the normal face-on.

Unpolarized Fresnel reflectance is the average of the s- and p-polarization reflectances.

Transmitted light uses the Beer–Lambert form T=exp(−σL). Note, however, that uTint and the per-channel absorption scale factors are visual tint parameters rather than SCHOTT spectral absorption data. The refractive index and Fresnel terms are tied to material data, but the absorption color is not claimed to be a physical material property.

6. Falling: A Quasi-Steady Model for a Thin Elliptical Disk

6.1 Mass, Inertia, and Added Mass

A 2D elliptical disk with semi-width a and semi-thickness b is interpreted with unit depth,

For the default leaf, I*=0.5824 and the broadside force-equilibrium velocity is 1.797 m/s.

6.2 Forces and Munk Moment

The relative velocity in body coordinates is written as (u,v), with angle of attack α=atan2(v,u).

The sign of the Munk moment matters.

const CD = P.CD0 + (P.CD90 - P.CD0) * sa * sa;
const CL = P.CL0 * Math.sin(2 * al);
const q = 0.5 * RHO_AIR * V * V * 2 * P.a;
const Fu = -q * (CD * ca - CL * sa);
const Fv = -q * (CD * sa + CL * ca);

const tau = -(P.m22 - P.m11) * u * v
          + P.kcp * P.a * ca * Fv
          - P.Ctau * RHO_AIR * P.a ** 4 * Math.abs(w) * w;

The state s=[x,y,θ,u,v,ω] is integrated with RK4 by subdividing each frame into four sub-steps; the out-of-plane tilt on screen comes from this θ rather than a separate sinusoidal term.

6.3 What Was Not Implemented

The quasi-steady C_L(α) and C_D(α) do not include trailing vortex shedding, which could not be implemented at this stage. This aspect will be revisited in the future after consulting additional literature.

7. Air: Divergence-Free Multi-Scale Velocity Field

The octave sum of the 2D stream function is taken as

and the velocity is constructed as

so that ∂u/∂x+∂v/∂y=0 is satisfied by construction. To match the Kolmogorov structure function δv∝ℓ^(1/3), the velocity amplitude is set to scale as k^(−1/3) and the stream-function amplitude as k^(−4/3),

const c = Math.cos(g.kx * x + g.ky * y + g.om * t + g.ph) * g.amp;
u += c * g.ky;
v += -c * g.kx;

The largest length scale across the 6 octaves is 0.3696 m and the smallest is 0.01493 m. The integral scale was capped at 0.37 m to prevent a single vortex larger than the 0.60 m viewport from pushing all leaves in the same direction. The random-phase sum is normalized by RMS rather than by amplitude sum, and the in-run measured value is 0.999.

8. Glass Particles and Shattered Dust

8.1 Handling the Regime Beyond Stokes Flow

The Stokes terminal velocity of a spherical particle is

however, larger shattered fragments do not stay within Re<1. Instead of clamping the velocity to an arbitrary upper bound, the force balance is solved by iterating the Schiller–Naumann correction.

The correction to the Stokes velocity is

and since Re depends on v, ten fixed-point iterations are performed.

for (let i = 0; i < 10; i++) {
    reynolds = Math.max(1e-8,
        (2 * RHO_AIR * radius * settle) / AIR_VISCOSITY);
    correction = 1 + 0.15 * reynolds ** 0.687;
    settle = vStokes / correction;
}

The computed results for the current radius range are as follows.

Particle velocity is updated using the exact one-frame solution dv/dt=(u_air−v)/τ as 1−exp(−Δt/τ), so even small particles with τ<Δt do not diverge numerically.

8.2 Physical Size vs. Visible Size

Particles in the 18–143 µm range are barely visible on a 0.60 m screen. Therefore, the physical radius is used only for drag and terminal velocity calculations, while the rendering half-size is scaled up to 0.45–1.45 mm for ambient particles and 0.70–3.10 mm for fragment particles. (This is purely a visual adjustment; the sizes are not rendered at true physical proportions.)

Each particle shape is rendered not as a plain circle but as a small facet with 3rd- and 5th-order angular-periodic deformations, and it samples the background through N-BK7 wavelength ratios to produce dispersion. This is a separate low-cost shader that does not evaluate the full 616-segment thickness field of the leaf for every particle.

9. Fracture Model

Glass fragments are formed not as clean triangles but as Voronoi cells from 12 seeds scattered around the centroid of the leaf.

float shardMask(vec2 p, float seed, float which) {
    float dm = length(p - shardSeed(which, seed));
    float best = 1e9;
    for (int i = 0; i < 12; i++) {
        if (float(i) == which) continue;
        best = min(best, length(p - shardSeed(float(i), seed)));
    }
    return best - dm;
}

The cell positions, fracture timing, and rupture impulses are all visual modeling choices; this is not fracture mechanics solving crack fronts, stress intensity factors, fracture energy, or impact contact.

10. Runtime Audit

10.1 GPU Optical Audit

audit.frag includes the same leaf-thickness.glsl as the rendering shader. It runs thicknessJet() directly on a 128² RGBA32F target and reads back the following values.

• Minimum and maximum values of det(I+d(n−1)H)

• Area of folded surface where det<0

• Area of caustic lines where |det|<0.15

• Central deflection d(n−1)|∇t|

• Leaf-interior mean of display gain

• Thickness minimum and maximum

Measured on 2026-08-18 with default settings in the Chrome browser, the values are as follows; this table is an observational record for regression verification, not a code constant.

FPS varies with device and window size (and therefore does not necessarily match this value exactly)

10.2 Fall statistics are not hardcoded

The actual integrated velocity of each leaf was mapped back into world coordinates to measure the instantaneous angle between the leaf's trajectory and the vertical.

That said, hardcoding this part at a later stage appears advantageous from a performance standpoint.

The statistics panel shows the median of up to 2,048 recent samples, while also displaying the net displacement summed across all leaves since the last statistics reset.

Changing the gust slider resets the statistics so that samples collected under the previous airflow are not mixed in.

11. Revision History: From Incorrect Implementation to Current Implementation

While working with AI, there were cases where the calculations in the implementation were actually wrong. Because I couldn't get the implementation right on the first try, I ended up measuring and recording incorrect values for those implementations.

The AI is far smarter than I am, but it doesn't seem to do exactly what I want on the first try. My guess is that it's caused by intermediate values vanishing when papers and formulas are fed in together.

12. Rendering Pipeline

background.frag → RGBA16F background
        │
        ├─ copy.frag → RGBA16F scene
        │
        ├─ leaf.frag
        │ ├─ N-BK7 RGB refraction
        │    ├─ Fresnel + Beer–Lambert
        │ └─ GPU-normalized caustics
        │
        ├─ particle.frag → glass dust and shatter glint
        │
        └─ post.frag → exposure / tone mapping → screen output

audit.frag → RGBA32F 128² → readPixels → audit table and caustic normalization

Leaves and fragments are instanced quads. Smaller, distant objects are drawn first, then closer objects overdraw them. Because the background is a screen-space texture, overlapping leaves do not refract one another again.

13. Controls

The panel is collapsed by default and left inert (to avoid thumbnail issues); opening it with the button reveals both the equation parameters and the current audit results.

14. Known Limitations

• Thin-element, single-lobe screen-space optics. Multiple preimages in folded regions are not path-traced.

• Lateral displacement inside the glass, multiple internal reflections, and polarization interference are not computed.

• The uTint absorption is a design parameter, not the measured spectral transmittance of N-BK7 glass.

• Overlapping leaves do not refract one another. All leaves sample the same original background.

• Optical thickness and aerodynamic thickness are intentionally decoupled.

• Aerodynamics follow a 2D quasi-steady model; there is no unsteady wake, vortex shedding, or 3D precession.

• The airflow field is a divergence-free synthetic field, not a Navier–Stokes DNS solution.

• The Schiller–Naumann particles assume spherical, isolated-particle approximations and do not represent the actual non-spherical drag of glass fragments.

• The displayed size of particles is larger than their physical size for pixel visibility (making them physically smaller looked worse visually).

• Fragmentation is a Voronoi visual model, not a stress-based fracture analysis.

• Leaf venation uses only three fixed seeds with inversions, yielding just six variants, so repetition is visible across many leaves.

• Spatial colonization is a morphogenesis algorithm, not a direct solution of the auxin transport PDE.

• WebGL2, EXT_color_buffer_float, and an HTTP-served environment are required.

16. References

1. SCHOTT, Optical Glass Datasheet — N-BK7.

2. L. J. Hickey, "Classification of the Architecture of Dicotyledonous Leaves," American Journal of Botany 60 (1973), 17–33.

3. A. Runions et al., "Modeling and visualization of leaf venation patterns," ACM Transactions on Graphics 24 (2005), 702–711, DOI: 10.1145/1073204.1073251.

4. C. D. Murray, "The Physiological Principle of Minimum Work," Journal of General Physiology 9 (1926), 835–841.

5. M. Born and E. Wolf, Principles of Optics, 7th ed., Cambridge University Press, 1999.

6. H. Lamb, Hydrodynamics, 6th ed., Cambridge University Press, 1932, §71.

7. M. M. Munk, The Aerodynamic Forces on Airship Hulls, NACA Report 184 (1924).

8. A. Andersen, U. Pesavento and Z. J. Wang, "Unsteady aerodynamics of fluttering and tumbling plates," Journal of Fluid Mechanics 541 (2005), 65–90.

9. S. B. Field et al., "Chaotic dynamics of falling disks," Nature 388 (1997), 252–254.

10. L. Schiller and A. Naumann, "A drag coefficient correlation," Zeitschrift des Vereines Deutscher Ingenieure 77 (1935), 318–320.

Personal note: the world has really come a long way. You just find a handful of similar papers, look up the relevant values from each one, plug them in, and out comes an animation. It's close to what I had in mind. I think it would be worth refining the formulas a bit more on the leaf venation side.