Constant KERNEL_FORCE_FIELD_SAMPLE 

pub const KERNEL_FORCE_FIELD_SAMPLE: &str = r#"
// Force field sampling kernel
// Reads particle positions, evaluates force fields, accumulates forces into velocity.

layout(local_size_x = 256) in;

struct Particle {
    vec4 pos_age;
    vec4 vel_life;
};

layout(std430, binding = 0) buffer Particles {
    Particle particles[];
};

// Force field types: 0=attractor, 1=repulsor, 2=vortex, 3=directional, 4=noise, 5=drag
struct ForceField {
    vec4 pos_strength;   // xyz = position, w = strength
    vec4 dir_radius;     // xyz = direction, w = radius
    vec4 params;         // x = falloff, y = frequency, z = type, w = unused
};

layout(std430, binding = 2) readonly buffer ForceFields {
    ForceField fields[];
};

uniform uint u_particle_count;
uniform uint u_field_count;
uniform float u_dt;
uniform float u_time;

float hash31(vec3 p) {
    p = fract(p * vec3(443.8975, 441.4230, 437.1950));
    p += dot(p, p.yzx + 19.19);
    return fract((p.x + p.y) * p.z);
}

vec3 eval_field(ForceField f, vec3 pos, float time) {
    vec3 field_pos = f.pos_strength.xyz;
    float strength = f.pos_strength.w;
    float radius = f.dir_radius.w;
    vec3 direction = f.dir_radius.xyz;
    float falloff = f.params.x;
    float freq = f.params.y;
    int ftype = int(f.params.z);

    vec3 delta = field_pos - pos;
    float dist = length(delta);
    float atten = 1.0;
    if (radius > 0.0) {
        atten = 1.0 - clamp(dist / radius, 0.0, 1.0);
        atten = pow(atten, falloff);
    }

    vec3 force = vec3(0.0);

    if (ftype == 0) {
        // Attractor: pull toward center
        if (dist > 0.001) {
            force = normalize(delta) * strength * atten;
        }
    } else if (ftype == 1) {
        // Repulsor: push away from center
        if (dist > 0.001) {
            force = -normalize(delta) * strength * atten;
        }
    } else if (ftype == 2) {
        // Vortex: swirl around axis (direction)
        vec3 axis = normalize(direction);
        vec3 radial = delta - dot(delta, axis) * axis;
        if (length(radial) > 0.001) {
            vec3 tangent = cross(axis, normalize(radial));
            force = tangent * strength * atten;
        }
    } else if (ftype == 3) {
        // Directional: constant force in a direction within radius
        force = normalize(direction) * strength * atten;
    } else if (ftype == 4) {
        // Noise: pseudo-random force based on position and time
        float n1 = hash31(pos * freq + vec3(time));
        float n2 = hash31(pos * freq + vec3(0.0, time, 0.0));
        float n3 = hash31(pos * freq + vec3(0.0, 0.0, time));
        force = (vec3(n1, n2, n3) * 2.0 - 1.0) * strength * atten;
    } else if (ftype == 5) {
        // Drag: opposes velocity (we approximate by opposing position change)
        force = -normalize(delta) * strength * atten * dist;
    }

    return force;
}

void main() {
    uint idx = gl_GlobalInvocationID.x;
    if (idx >= u_particle_count) return;

    Particle p = particles[idx];
    vec3 pos = p.pos_age.xyz;
    vec3 vel = p.vel_life.xyz;

    // Skip dead particles
    if (p.pos_age.w >= p.vel_life.w) return;

    // Accumulate forces from all fields
    vec3 total_force = vec3(0.0);
    for (uint i = 0u; i < u_field_count; i++) {
        total_force += eval_field(fields[i], pos, u_time);
    }

    // Apply accumulated force
    vel += total_force * u_dt;
    particles[idx].vel_life.xyz = vel;
}
"#;

Force field sampling kernel: evaluate multiple force fields at particle positions.