Each leaf, idealized as a thin disk, responds to wind and gravity. The air's velocity field is a space-time varying sinusoid combined additively with randomly-spawning Gaussian-enveloped vortices (referred to as swirls).

Given the relative velocity between each leaf and the surrounding air, we can compute the forces and moments that act upon each leaf. From there, we can apply Newton's second law and integrate to propagate the state of each leaf over time.

When a leaf wanders too far out of bounds, it fades and despawns, increasing the probability that a new leaf will spawn near the top of the domain.

Except where stated otherwise, everything runs on the $\left[\text{kg}\cdot\text{m}\cdot\text{s}\right]$ unit system.
$\uvec{i}, \uvec{j}, \uvec{k}$ are unit vectors in the $x, y, z$ directions, respectively
In website space, $\uvec{i}, \uvec{j}, \uvec{k}$ point right, up, and into the page, respectively
Bold symbols such as $\vec{v}$ are vectors in $\mathbb{R}^3$
Hatted symbols such as $\uvec{v}$ are unit vectors
$t$ is time
$\vec{x} = [x,y,z]$ denotes spatial coordinates
$\omega$ is a frequency in $[1/s]$
$\gamma$ is a wavevector (spatial version of frequency) in $[1/m]$
$\propto$ means "proportional to" (i.e. I'm omitting a scaling factor)
$\rho$ signifies density (could be an area density $\left[\text{kg}/\text{m}^2\right]$ or volumetric density $\left[\text{kg}/\text{m}^3\right]$)

Linear kinematics are pretty simple in this simulation. Each leaf has a position, velocity, and acceleration: $\vec{x}_{\text{leaf}}$, $\vec{v}_{\text{leaf}}$, and $\vec{a}_{\text{leaf}}.$ It's easy to propagate these quantities in time because they live in $\mathbb{R}^3$ — meaning that they compose through addition. My simulation uses symplectic Euler for this.

Rotational kinematics are not as simple to model because rotations do not compose through addition (i.e. $\,\boldsymbol{\theta}_{t+\Delta t}\leftarrow \boldsymbol{\theta}_t+\boldsymbol{\omega}\Delta t\,$ doesn't work). Instead, we need to represent rotations and track orientations with quaternions, which live in $\mathbb{R}\times\mathbb{R}^3$. It may seem like a lot, but seeing the leaves rotate is so worth it!

My sim uses quaternions for rotation math.

A quaternion is an ordered pair \(q=(w,\vec{q})\in\mathbb{R}\times\mathbb{R}^3\) with scalar part \(w\) and vector part \(\vec{q}\). We write the quaternion product as \(\otimes\).

A unit quaternion satisfies $q \otimes q^{*} = (\|q\|^{2},\,\vec{0}) = (1,\,\vec{0})$. Unit quaternions parametrize 3D rotations. To rotate a vector $\vec{v}$ by quaternion $\vec{q}$, we embed \(\vec{v}\in\mathbb{R}^3\) as a pure quaternion \((0,\vec{v})\) and sandwich multiply: $$\mathcal{R}_q(\vec{v}) \;=\; \big(q \otimes (0,\vec{v}) \otimes q^{-1}\big)$$ where, since the result is a quaternion in $\mathbb{R}\times\mathbb{R}^3$, we take only the vector part as the result of the rotation.

We express angular acceleration $\alpha$ and angular velocity $\omega$ in the inertial frame (as opposed to the local frame), and our quaternion \(q\) maps local→world.

We integrate angular velocity by treating it as constant over each small step. The kinematic ODE \[ \dot q=\tfrac12\,(0,\boldsymbol{\omega})\otimes q \] has the per-step solution \[ q_{n+1}=\underbrace{\exp\!\big(\tfrac12(0,\boldsymbol{\omega}_{n+1})\Delta t\big)}_{\text{unit quaternion with axis }\hat{\boldsymbol{\omega}}_{n+1}\text{, and angle }\|\boldsymbol{\omega}_{n+1}\|\Delta t}\ \otimes\ q_n. \] This update rule can be read more simply as $q\leftarrow dq\otimes q$.

The background field is a space-time varying sinusoid:
$$\vec{v}_{\text{background}}(\mathbf{x},t)\propto\sin(2\pi\gamma\,y)\sin(2\pi\omega\,t)\uvec{i}.$$

Each swirl is a Gaussian vortex about the $\uvec{k}$-axis, with on-off behavior modulated by smoothed Heaviside envelopes: $$\vec{v}_{\text{swirl}}(\mathbf{x},t)\propto e^{-\tfrac12(r/R)^2}\operatorname{env}(t)\symuvec{\phi}, \quad\symuvec{\phi}=\uvec{k}\times\uvec{r}.$$ where $\operatorname{env}(t,t_{0},t_{1}) = \operatorname{H_{s,2}}(t-t_{0})\operatorname{H_{s,2}}(t_{1} - t)$ turns the vortex on and off with a smoothed Heaviside function at times $t_{0}$ and $t_{1}$ respectively. Here, position vector $\vec{r} = \vec{x}-\vec{x_{\text{swirl}}}$, where $\vec{x_{\text{swirl}}}$ is the center of the swirl. Distance from the swirl center is $r = \norm{\vec{r}}$. Radial size variable $R$ scales the spatial size of the swirl.

The total velocity field is the background field, plus any existing swirls: $$\vec{v}_{\text{air}}(\vec{x},t) = \vec{v}_{\text{background}}(\vec{x},t) + \sum_{\text{swirls}}{\vec{v}_{\text{swirl}}(\vec{x},t)}.$$

The relative velocity between the air and each leaf is computed as:

By idealizing each leaf as a thin disk, we can compute the aerodynamic forces and moments analytically. The aerodynamic force is computed as:

The coefficients of drag depend on the angle $\theta$ between the relative velocity $(\uvec{v}_{\text{rel}})$ and the leaf normal $(\uvec{n}_{\text{leaf}})$:

We offset the center of pressure a small fraction of the radius along the in-plane flow direction, then compute the moment by a lever-arm $\times$ force cross-product. $$\vec{r} = -a_{\mathrm{CoP}}R\uvec{v}_{\text{rel},\parallel}$$

Let $\vec{x}_\star$ denote the camera's current target position (computed directly from scroll location). To move the camera to that location smoothly, we apply $$\vec{x}_{t+\Delta t} = \vec{x}_\star + (\vec{x}_{t} - \vec{x}_\star)e^{-\tau\Delta t}$$ where damping rate $\tau > 0$ controls how snappy/laggy the transition feels.

Initially, the scene is seeded with $N_{0}$ leaves. Over time, aerodynamics or gravity eventually cause them to leave the simulated domain and are culled.

New leaves are probabilistically spawned via a Poisson distribution. The number of leaves that will spawn during a given time-step is drawn as: $$N_{\text{new}} \sim \operatorname{Poisson}(\lambda \Delta t)$$ $$\lambda = \lambda_0 \operatorname{max}(0,N_{\text{max}} - N)/N_{\text{max}}$$ By this construction, the spawn rate $\lambda$ linearly interpolates between $0$ and the base spawn rate $\lambda_{0}$ based on the current population $N$ and the maximum population $N_{\text{max}}$. If $N \geq N_{\text{max}}$, then $\lambda = 0$ (no new leaves will spawn). If $N = 0$, then $\lambda = \lambda_{0}$ (leaves will spawn at the base spawn rate). This spawn logic is frame-rate independent.

Swirls follow the same spawning logic, but are culled differently from leaves. Each swirl has a set time duration, after which it fades out.