skills/shader-dev/reference/particle-system.md

Particle Systems — Detailed Reference

This document is a detailed supplement to SKILL.md, containing prerequisites, in-depth explanations of each step, variant details, performance optimization analysis, and complete code for combination suggestions.

NOTE: Code examples in this document primarily target the ShaderToy environment. For standalone HTML deployment, refer to the WebGL2 single-file template in SKILL.md, which includes complete HTML + JS + GLSL code.

Prerequisites

• GLSL basic syntax (uniforms, varyings, built-in functions)
• 2D/3D vector math (dot product, cross product, normalization, matrix rotation)
• ShaderToy architecture (mainImage, iTime, iResolution, iChannel, multi-Buffer passes)
• Basic physics concepts: velocity = derivative of position, acceleration = force/mass
• Usage of texelFetch (precise pixel data reading from Buffers)

Implementation Steps in Detail

Step 1: Hash Random Functions

What: Define a function that generates pseudo-random numbers from a float (particle ID). This is the foundational infrastructure for all particle systems.

Why: Each particle needs unique but deterministic attributes (color, size, initial direction, etc.); hash functions provide repeatable "randomness".

Three hash function dimensions are provided:

• hash11: 1D → 1D, for scalar randomness (lifetime, brightness, etc.)
• hash12: 1D → 2D, for 2D randomness (initial position, etc.)
• hash33: 3D → 3D, for 3D velocity perturbation

// Standard 1D -> 1D hash, returns [0, 1)
float hash11(float p) {
    return fract(sin(p * 127.1) * 43758.5453);
}

// 1D -> 2D hash, for 2D randomness
vec2 hash12(float p) {
    vec3 p3 = fract(vec3(p) * vec3(0.1031, 0.1030, 0.0973));
    p3 += dot(p3, p3.yzx + 33.33);
    return fract((p3.xx + p3.yz) * p3.zy);
}

// 3D -> 3D hash, for 3D velocity perturbation
vec3 hash33(vec3 p) {
    p = fract(p * vec3(443.897, 397.297, 491.187));
    p += dot(p.zxy, p.yxz + 19.19);
    return fract(vec3(p.x * p.y, p.z * p.x, p.y * p.z)) - 0.5;
}

Step 2: Particle Lifecycle Management

What: Compute birth time, lifespan, current age for each particle, and auto-respawn after death.

Why: Lifecycle is the core mechanism of particle systems — the cycle of birth, motion, fade-out, death, and respawn. fract or mod implements infinite cycling without additional state.

Key design:

• spawnTime: Each particle's birth time differs, generated by hash11 from the ID, spread across the [0, START_TIME] interval
• lifetime: Each particle's lifespan differs, random within the [LIFETIME_MIN, LIFETIME_MAX] interval
• mod(time - spawnTime, lifetime): Automatic cycling; the particle respawns immediately after death
• floor(...) computes the current life cycle number, used to generate different random attributes each cycle

#define NUM_PARTICLES 100   // adjustable: particle count
#define LIFETIME_MIN 1.0    // adjustable: minimum lifespan (seconds)
#define LIFETIME_MAX 3.0    // adjustable: maximum lifespan (seconds)
#define START_TIME 2.0      // adjustable: time for all particles to be born

// Returns: x = current normalized age [0,1], y = current life cycle number
vec2 particleAge(int id, float time) {
    float spawnTime = START_TIME * hash11(float(id) * 2.0);
    float lifetime = mix(LIFETIME_MIN, LIFETIME_MAX, hash11(float(id) * 3.0 - 35.0));
    float age = mod(time - spawnTime, lifetime);
    float run = floor((time - spawnTime) / lifetime);
    return vec2(age / lifetime, run);
}

Step 3: Stateless Particle Position Computation

What: Compute 2D/3D position solely from particle ID and time, without relying on any Buffer.

Why: For decorative effects (starfields, fireworks, orbiting light points), the stateless approach is simplest and most efficient. Define the main trajectory via parametric curves (e.g., Lissajous curves), then add random offset and gravity.

Position is composed of three components:

1. Main trajectory (harmonic oscillator): Multiple cosine waves superimposed to form smooth Lissajous curves, controlling the overall motion path of the particle group
2. Random drift: Each particle linearly diffuses from the main trajectory position over time; DRIFT_MAX controls the diffusion range
3. Gravity: Parabolic descent via 0.5 * g * t²; age² is the normalized form of time

#define GRAVITY vec2(0.0, -4.5)     // adjustable: gravity direction and strength
#define DRIFT_MAX vec2(0.28, 0.28)  // adjustable: maximum random drift amplitude

// Harmonic superposition for smooth main trajectory
float harmonics(vec3 freq, vec3 amp, vec3 phase, float t) {
    float val = 0.0;
    for (int h = 0; h < 3; h++)
        val += amp[h] * cos(t * freq[h] * 6.2832 + phase[h] / 360.0 * 6.2832);
    return (1.0 + val) / 2.0;
}

vec2 particlePosition(int id, float time) {
    vec2 ageInfo = particleAge(id, time);
    float age = ageInfo.x;
    float run = ageInfo.y;

    // Main trajectory (harmonic oscillator)
    float slowTime = time * 0.1; // time along main trajectory
    vec2 mainPos = vec2(
        harmonics(vec3(0.4, 0.66, 0.78), vec3(0.8, 0.24, 0.18), vec3(0.0, 45.0, 55.0), slowTime),
        harmonics(vec3(0.415, 0.61, 0.82), vec3(0.72, 0.28, 0.15), vec3(90.0, 120.0, 10.0), slowTime)
    );

    // Random drift (grows linearly with time)
    vec2 drift = DRIFT_MAX * (vec2(hash11(float(id) * 3.0 + run * 4.0),
                                    hash11(float(id) * 7.0 - run * 2.5)) - 0.5) * age;
    // Gravity effect
    vec2 grav = GRAVITY * age * age * 0.5;

    return mainPos + drift + grav;
}

Step 4: Buffer-Stored Particle State (Stateful System)

What: Use one row of pixels in a Buffer texture to store all particles, with each pixel = one particle's (pos.x, pos.y, vel.x, vel.y).

Why: When inter-frame persistent state is needed (physics collisions, force field interactions, N-body simulations), particle data must be written to a Buffer and read back the next frame. In ShaderToy, each pixel is a storage cell.

Design points:

• fragCoord.y > 0.5: Only the first row of pixels stores particles; remaining pixels are discarded
• fragCoord.x corresponds to particle ID; each pixel's RGBA stores (pos.x, pos.y, vel.x, vel.y)
• iFrame < 5: First few frames are initialization, randomly distributing particle positions
• Force accumulation: boundary repulsion + inter-particle attraction/repulsion + friction
• Clamp velocity and acceleration after integration to prevent numerical explosion

// === Buffer A: Particle physics update ===
#define NUM_PARTICLES 40    // adjustable: particle count
#define MAX_VEL 0.5         // adjustable: maximum velocity
#define MAX_ACC 3.0         // adjustable: maximum acceleration
#define RESIST 0.2          // adjustable: drag coefficient
#define DT 0.03             // adjustable: time step

// Read the i-th particle's data
vec4 loadParticle(float i) {
    return texelFetch(iChannel0, ivec2(i, 0), 0);
}

void mainImage(out vec4 fragColor, in vec2 fragCoord) {
    if (fragCoord.y > 0.5 || fragCoord.x > float(NUM_PARTICLES)) discard;

    float id = floor(fragCoord.x);
    vec2 res = iResolution.xy / iResolution.y;

    // Initialization
    if (iFrame < 5) {
        vec2 rng = hash12(id);
        fragColor = vec4(0.1 + 0.8 * rng * res, 0.0, 0.0);
        return;
    }

    // Read current state
    vec4 particle = loadParticle(id); // xy = pos, zw = vel
    vec2 pos = particle.xy;
    vec2 vel = particle.zw;

    // === Force accumulation ===
    vec2 force = vec2(0.0);

    // Boundary repulsion force
    force += 0.8 * (1.0 / abs(pos) - 1.0 / abs(res - pos));

    // Inter-particle interaction (attraction/repulsion)
    for (float i = 0.0; i < float(NUM_PARTICLES); i++) {
        if (i == id) continue;
        vec4 other = loadParticle(i);
        vec2 w = pos - other.xy;
        float d = length(w);
        if (d > 0.0)
            force -= w * (6.3 + log(d * d * 0.02)) / exp(d * d * 2.4) / d;
    }

    // Friction force
    force -= vel * RESIST / DT;

    // === Integration ===
    vec2 acc = force;
    float a = length(acc);
    acc *= a > MAX_ACC ? MAX_ACC / a : 1.0; // limit acceleration

    vel += acc * DT;
    float v = length(vel);
    vel *= v > MAX_VEL ? MAX_VEL / v : 1.0; // limit velocity

    pos += vel * DT;

    fragColor = vec4(pos, vel);
}

Step 5: Particle Rendering — Point Light / Metaball Style

What: Iterate over all particles in the Image pass and render each as a soft glowing point.

Why: 1/dot(p,p) produces natural inverse-square distance falloff; when multiple particles overlap, the result resembles metaball fusion. This is the most classic particle rendering method.

Rendering principle:

• dot(p, p) is dist²; using it as the denominator produces inverse-square distance falloff
• BRIGHTNESS controls point size — larger values produce bigger glow points
• totalWeight accumulates the metaball contribution of all particles
• Color interpolates between COLOR_START and COLOR_END based on particle velocity
• mix(col, pcol, mb / totalWeight) implements contribution-weighted color blending, with nearby particles having higher color weight
• Final normalize + clamp prevents overexposure

#define BRIGHTNESS 0.002        // adjustable: particle brightness
#define COLOR_START vec3(0.0, 0.64, 0.2)  // adjustable: start color
#define COLOR_END vec3(0.06, 0.35, 0.85)  // adjustable: end color

vec3 renderParticles(vec2 uv) {
    vec3 col = vec3(0.0);
    float totalWeight = 0.0;

    for (int i = 0; i < NUM_PARTICLES; i++) {
        vec4 particle = loadParticle(float(i));
        vec2 p = uv - particle.xy;

        // Metaball-style falloff: radius / distance²
        float mb = BRIGHTNESS / dot(p, p);
        totalWeight += mb;

        // Interpolate color based on particle attributes
        float ratio = length(particle.zw) / MAX_VEL;
        vec3 pcol = mix(COLOR_START, COLOR_END, ratio);
        col = mix(col, pcol, mb / totalWeight);
    }

    totalWeight /= float(NUM_PARTICLES);
    col = normalize(col) * clamp(totalWeight, 0.0, 0.4);
    return col;
}

Step 6: Frame Feedback Motion Blur

What: Blend the current frame with the previous frame's Buffer to produce motion trails.

Why: Single-frame particles are just discrete dots; through temporal accumulation (feedback blending), continuous trails/afterimage effects are produced. The blend coefficient controls trail length.

Design points:

• TRAIL_DECAY closer to 1 produces longer trails (0.99 = very long trail, 0.9 = short trail)
• Requires an extra Buffer pass: Buffer B handles trail accumulation, Image pass reads from Buffer B
• prev * TRAIL_DECAY + current: Decay old frame + overlay new frame
• This method can also simulate high particle density with few particles + long trails

#define TRAIL_DECAY 0.95  // adjustable: trail decay rate, closer to 1 = longer trail

// In the rendering Buffer's mainImage:
void mainImage(out vec4 fragColor, in vec2 fragCoord) {
    vec2 uv = fragCoord / iResolution.xy;

    // Read previous frame
    vec3 prev = texture(iChannel0, uv).rgb * TRAIL_DECAY;

    // Draw current frame particles
    vec3 current = renderParticles(fragCoord / iResolution.y);

    // Overlay
    fragColor = vec4(prev + current, 1.0);
}

Step 7: HSV Coloring and Star Glare Effects

What: Color particles using HSV color space and add cross/star diffraction spike lines.

Why: HSV makes it easy to rotate hue for rainbow effects; star glare (diffraction spikes) simulates real lens optical effects, giving light points more visual quality.

HSV coloring principle:

• hsv.x = Hue, 0-1 maps to one revolution of the color wheel
• hsv.y = Saturation, 0 = gray, 1 = pure color
• hsv.z = Value, 0 = black, 1 = brightest
• Cosine waves approximate the RGB channel hue response curves

Star glare principle:

• Star glare is caused by diffraction from lens aperture blades in real photography
• Implemented by stretching the distance field in specific directions: one horizontal, one vertical, one at each 45-degree diagonal
• stretch parameter controls the stretch ratio; larger values produce thinner, longer lines
• 0.707 is the approximation of cos(45°) = sin(45°), used to rotate to diagonal directions

// HSV -> RGB conversion
vec3 hsv2rgb(vec3 c) {
    vec4 K = vec4(1.0, 2.0 / 3.0, 1.0 / 3.0, 3.0);
    vec3 p = abs(fract(c.xxx + K.xyz) * 6.0 - K.www);
    return c.z * mix(K.xxx, clamp(p - K.xxx, 0.0, 1.0), c.y);
}

// Star glare effect: produces elongated light rays in horizontal/vertical/diagonal directions
float starGlare(vec2 relPos, float intensity) {
    // Horizontal/vertical branches
    vec2 stretch = vec2(9.0, 0.32); // adjustable: stretch ratio
    float dh = length(relPos * stretch);
    float dv = length(relPos * stretch.yx);

    // Diagonal branches
    vec2 diagPos = 0.707 * vec2(dot(relPos, vec2(1, 1)), dot(relPos, vec2(1, -1)));
    float dd1 = length(diagPos * vec2(13.0, 0.61));
    float dd2 = length(diagPos * vec2(0.61, 13.0));

    float glare = 0.25 / (dh * 3.0 + 0.01)
                + 0.25 / (dv * 3.0 + 0.01)
                + 0.19 / (dd1 * 3.0 + 0.01)
                + 0.19 / (dd2 * 3.0 + 0.01);
    return glare * intensity;
}

Common Variant Details

Variant 1: Metaball Polar Coordinate Particles

Difference from base version: Particles are uniformly distributed in polar coordinates and expand outward, using 1/dot(p,p) metaball fusion instead of point lights, producing organic "blob-like" effects.

Design points:

• Particle positions change from Cartesian to polar coordinates: angles uniformly distributed around the circle, distance cycles with fract over time
• fract(time * speed + hash) produces particles expanding from center outward then respawning
• Metaball rendering: 0.84 / dot(p, p) values naturally accumulate where particles overlap, forming fused organic shapes
• Color interpolates between start and end colors based on distance d
• mb / totalSum ensures color blending is weighted by contribution

// Particle position changed to polar coordinate expansion
float d = fract(time * 0.51 + 48934.4238 * sin(float(i) * 692.7398));
float angle = TAU * float(i) / float(NUM_PARTICLES);
vec2 particlePos = d * vec2(cos(angle), sin(angle)) * 4.0;

// Metaball rendering replaces point light
vec2 p = uv - particlePos;
float mb = 0.84 / dot(p, p);  // adjustable: 0.84 = metaball radius
col = mix(col, mix(startColor, endColor, d), mb / totalSum);

Variant 2: Buffer Storage + Boids Flocking Behavior

Difference from base version: Changes from stateless to stateful, with each particle stored in a Buffer pixel, enabling N-body attraction/repulsion interaction and Boids emergent behavior.

Design points:

• Each particle iterates over all other particles, computing the net attraction/repulsion force
• The force formula (6.3 + log(d² * 0.02)) / exp(d² * 2.4) produces:
  • Short-range repulsion (exponential decay dominates)
  • Medium-range attraction (logarithmic term dominates)
  • Long-range no effect (exponential decay approaches zero)
• Friction vel * 0.2 / dt prevents infinite acceleration
• Overall effect: particles self-organize into group motion patterns, exhibiting fish-school/bird-flock emergent behavior

// Buffer A: force accumulation
vec2 sumForce = vec2(0.0);
for (float j = 0.0; j < NUM_PARTICLES; j++) {
    if (j == id) continue;
    vec4 other = texelFetch(iChannel0, ivec2(j, 0), 0);
    vec2 w = pos - other.xy;
    float d = length(w);
    // Combined attraction+repulsion: short-range repulsion, long-range attraction
    sumForce -= w * (6.3 + log(d * d * 0.02)) / exp(d * d * 2.4) / d;
}
sumForce -= vel * 0.2 / dt; // friction

Variant 3: Verlet Integration Cloth Simulation

Difference from base version: Particles are connected through spring constraints (grid topology), using Verlet integration instead of Euler method — no need to explicitly store velocity.

Design points:

• Verlet integration: newPos = 2 * current - previous + acc * dt²
  • Velocity is implicit in current - previous
  • No separate velocity storage needed; RGBA can store (current.xy, previous.xy)
  • More stable than Euler in constraint solving (won't blow up from high-frequency oscillation)
• Spring constraints: each pair of adjacent particles has a "rest length"
  • Compute the difference between current distance and rest length
  • Move particles toward the rest length by a small step (0.1 is the relaxation coefficient)
  • Multiple constraint-solving iterations converge to a stable state
• Grid topology: particle IDs arranged in rows and columns, each connected to its up/down/left/right neighbors

// Verlet integration: velocity is implicit in (current position - previous position)
// particle.xy = current position, particle.zw = previous position
vec2 newPos = 2.0 * particle.xy - particle.zw + vec2(0.0, -0.6) * dt * dt;
particle.zw = particle.xy;
particle.xy = newPos;

// Spring constraint solving
vec4 neighbor = texelFetch(iChannel0, neighborId, 0);
vec2 delta = neighbor.xy - particle.xy;
float dist = length(delta);
float restLength = 0.1; // adjustable: rest length
particle.xy += 0.1 * (dist - restLength) * (delta / dist);

Variant 4: 3D Particles + Ray Rendering

Difference from base version: Particles are stored in 3D space; rendering uses rays cast from the camera, computing the closest distance from each ray to each particle for coloring.

Design points:

• Camera at (0, 0, 2.5), ray direction determined by screen UV
• Point-to-line distance formula: |cross(P-O, D)|, where O is ray origin, D is ray direction, P is particle position
• dot(cross(...), cross(...)) computes the squared distance (avoiding sqrt)
• × 1000.0 is the distance scaling factor controlling visual particle size
• Difference from 2D rendering: 2D uses length of uv - pos, 3D uses closest distance from ray to point

// 3D rendering: closest distance from ray to particle
vec3 ro = vec3(0.0, 0.0, 2.5);
vec3 rd = normalize(vec3(uv, -0.5));
for (int i = 0; i < numParticles; i++) {
    vec3 pos = texture(iChannel0, vec2(i, 100.0) * w).rgb;
    // Squared distance from point to line
    float d = dot(cross(pos - ro, rd), cross(pos - ro, rd));
    d *= 1000.0;
    float glow = 0.14 / (pow(d, 1.1) + 0.03);
    col += glow * particleColor;
}

Variant 5: Raindrop Particles (3D Scene Integration)

Difference from base version: Particles move in 3D world space (gravity + wind + jitter), rendered as screen-space water drops with normal mapping to simulate refraction. Respawned randomly after lifecycle ends.

Design points:

• speedScale includes sin(π/2 * pow(age/lifetime, 2)) to accelerate the fall
• Wind force is projected onto the camera's right/up directions via dot product
• Jitter randVec2 * jitterSpeed simulates air turbulence
• Death and respawn: particle.z accumulates age; when it exceeds particle.a (lifespan), position and lifespan are reset
• Rendering can overlay raindrop SDF + refraction normal mapping to simulate realistic raindrop optics

// 3D force accumulation
float speedScale = 0.0015 * (0.1 + 1.9 * sin(PI * 0.5 * pow(age / lifetime, 2.0)));
particle.x += (windShieldOffset.x + windIntensity * dot(rayRight, windDir)) * fallSpeed * speedScale * dt;
particle.y += (windShieldOffset.y + windIntensity * dot(rayUp, windDir)) * fallSpeed * speedScale * dt;
// Jitter
particle.xy += 0.001 * (randVec2(particle.xy + iTime) - 0.5) * jitterSpeed * dt;
// Death and respawn
if (particle.z > particle.a) {
    particle.xy = vec2(rand(seedX), rand(seedY)) * iResolution.xy;
    particle.a = lifetimeMin + rand(pid) * (lifetimeMax - lifetimeMin);
    particle.z = 0.0;
}

Variant 6: Vortex/Storm Particle System

Difference from base version: Particles move along spiral trajectories, forming storm/sandstorm/blizzard effects. Stateless single pass.

Design points:

• Differential rotation: inner circles rotate faster than outer circles (angularSpeed = k / (offset + radius)), producing natural vortices
• Particle color must be significantly brighter than the background (2-3x), otherwise invisible against similar-colored backgrounds
• Brightness budget: with 150 particles, numerator=0.002, epsilon=0.008 (peak=0.25) is safe
• Vortex center dark zone implemented with smoothstep(innerR, outerR, dist) mask

Variant 7: Meteor/Trailing Line Rendering

Difference from base version: Particles are rendered as elongated glow lines rather than circular light points.

Design points:

• Must have a clearly visible static starfield background: Call starField() function; stars rendered as sharp Gaussian points using exp(-dist²*k), with peak brightness >= 0.3
• Meteor trail must be bright enough: core multiplier >= 0.15; after dividing by sample count, each step still needs >= 0.005 visible contribution
• Do not use 1/(distPerp² + tiny_epsilon) for lines — use exp(-distPerp / width) for safe glow
• Meteor head headGlow = 0.005 / (dist² + 0.0008) ensures bright visibility
• trailLen range 0.15-0.35 ensures sufficient trail length

Variant 8: Fountain/Upward Jet Particle System

Difference from base version: Particles jet upward from a single point, follow parabolic arcs, then fall back. Stateless single pass.

Design points:

• Must include three layers: (1) Main water column particles (upward jet + parabola) (2) Splash particles (flying sideways upon hitting water) (3) Water surface/pool visuals
• Particles must be sharp, visible individual points: Use small epsilon (<= 0.002) with small numerator; must not only produce diffuse glow
• Parabolic motion: pos = base + vel0 * t + 0.5 * gravity * t²
• Ground clipping: if (pos.y < waterLevel) continue;
• Brightness budget: 60 main particles + 40 splash particles, each with epsilon in the 0.001-0.002 range

Variant 10: Spiral Array/Magic Particle System

Difference from base version: Particles arranged in spiral or geometric arrays, producing magic circles, magic dust, and similar effects. Particles feature iridescent color variation. Stateless single pass.

Design points:

• Must have discrete visible particles: Each particle must be an individually visible small light point, not just blurry glow blobs. Use small epsilon (0.0004-0.0006) for sufficient sharpness
• Spiral trajectory: angle = baseAngle + norm * spiralTurns + time * rotSpeed, radius increases with norm
• Magic circles use independent ring particle layers with uniformly distributed angles + time rotation, using elliptical projection to simulate 3D perspective
• Iridescent effect: hue = fract(particleId / total + time * speed + norm * shift), hue varies continuously with ID and time, covering the full color wheel
• Starlight shimmer: shimmer = 0.7 + 0.3 * sin(time * freq + particleId * phase) controls each particle's brightness pulsation
• Two-layer structure: (1) Spiral ascending particle stream (2) Horizontally rotating magic circle light point rings

Brightness Budget Quick Reference

Single pass system: N × (numerator / epsilon) < 5.0

Particle Count	Recommended numerator	Recommended epsilon	Single Particle Peak	Total Peak (no fade)
40	0.015	0.03	0.5	20 → Reinhard OK
80	0.008	0.015	0.53	42 → Reinhard OK
150	0.002	0.008	0.25	37 → Reinhard OK
200	0.001	0.005	0.2	40 → Reinhard OK

Multi-pass ping-pong system: N × (numerator / epsilon) × 1/(1-decay) < 10.0

decay	Amplification Factor	20 Particle Peak Limit	50 Particle Peak Limit	100 Particle Peak Limit
0.88	8.3x	0.06	0.024	0.012
0.92	12.5x	0.04	0.016	0.008
0.95	20x	0.025	0.01	0.005

Performance Optimization In-Depth Analysis

1. Particle Count and Loop Overhead

• Bottleneck: Every pixel iterates over all particles (O(W×H×N)); particle count is the biggest performance lever.
• Optimization: Reducing particle count from 200 to 80 may have little visual difference but doubles performance. Early exit optimization can also help:

float dist = length(uv - ppos);
if (dist > 0.1) continue; // adjustable: skip particles beyond influence range

• Note: The early exit threshold must be tuned based on particle brightness / influence radius; too small causes abrupt particle edge cutoff

2. Frame Feedback as Substitute for High Particle Count

• Technique: Few particles + frame feedback trails (prev * 0.95 + current) visually equals many more particles. Drawing 50 particles per frame + accumulation = visual density far exceeding 50.
• This approach has the additional benefit of producing natural motion blur
• Requires an additional Buffer pass for accumulated frames

3. N-body Interaction Complexity

• Bottleneck: Each particle interacts with all others = O(N²). Becomes very slow when N > 100.
• Optimization A: Only interact with K nearest neighbors (using Voronoi tracking acceleration structure, see "Combining with Voronoi Spatial Acceleration Structure" below).
• Optimization B: Divide space into grid cells, only check particles in adjacent cells. Implementing the grid on GPU requires additional Buffer passes to maintain grid data.

4. Sub-frame Stepping

• Problem: High-speed particles move multiple pixels per frame, leaving discontinuous trajectories.
• Optimization: Perform multiple small steps per frame for each particle, accumulating rendering along the way:

const int stepsPerFrame = 7; // adjustable
for (int j = 0; j < stepsPerFrame; j++) {
    // Render particle contribution at this position
    pos += vel * 0.002 * 0.2;
}
col /= float(stepsPerFrame);

• More sub-frames produce more continuous trajectories but linearly increase computational cost
• Suitable for firework explosions, high-speed bullet curtains, etc.

5. Precision and Numerical Stability

• Velocity and acceleration need clamping to prevent numerical explosion:

float v = length(vel);
vel *= v > MAX_VEL ? MAX_VEL / v : 1.0;

• Verlet integration is more stable than Euler in constraint solving, especially for cloth and spring networks
• For long-running simulations, be aware of floating-point precision errors accumulating over time

Combination Suggestions with Complete Code

Combining with Raymarching Scenes

Particle systems are often embedded in Raymarching scenes (e.g., rain, sparks, dust). Method: During the Raymarching step loop, simultaneously sample particle density/positions and overlay onto scene color. Or render particles to a separate Buffer and blend during final compositing.

Combining with Noise / Flow Fields

Use Simplex/Perlin noise to generate a velocity field; particles move along the noise gradient:

// Use noise to drive particle velocity
vel += hash33(vel + time) * 7.0; // random perturbation
vel = mix(vel, -pos * pow(length(pos), 0.75), 0.5 + 0.5 * sin(time)); // center attraction

This combination is suitable for "neural synapse", "smoke flow", and other organic effects.

Combining with Post-Processing

• Bloom: Apply Gaussian blur to particle rendering output and overlay, enhancing the glow.
• Chromatic Aberration: Offset-sample R/G/B channels separately, adding a lens effect.
• Tone Mapping: Apply Reinhard mapping col = col / (1.0 + col) to HDR particle brightness.

Combining with SDF Shape Rendering

Render particles as specific SDF shapes (fish, water drops, sparks) instead of abstract light points. Method: Rotate local coordinates based on particle velocity direction, then compute SDF distance in that coordinate system:

float sdFish(vec2 p, float angle) {
    float c = cos(angle), s = sin(angle);
    p *= 20.0 * mat2(c, s, -s, c);
    return max(min(length(p), length(p - vec2(0.56, 0.0))) - 0.3, -min(length(p - vec2(0.8, 0.0)) - 0.45, length(p + vec2(0.14, 0.0)) - 0.12)) * 0.05;
}

Combining with Voronoi Spatial Acceleration Structure

For large-scale particles (thousands), use Voronoi tracking acceleration structure instead of brute-force traversal. Each pixel maintains the IDs of the 4 nearest particles, updated through neighborhood propagation. This reduces rendering and physics queries from O(N) to O(1) (fixed neighborhood query per pixel). Suitable for fluid simulation and large-scale swarm behavior.