.0973));\\n p3 += dot(p3, p3.yzx+33.33);\\n return fract((p3.xx+p3.yz)*p3.zy);\\n}\\n\\nvec3 udir(vec2 rng)\\n{\\n float phi = 2.*PI*rng.x;\\n float ctheta = 2.*rng.y-1.;\\n float stheta = sqrt(1.0-ctheta*ctheta);\\n return vec3(cos(phi)*stheta, sin(phi)*stheta, ctheta);\\n}\\n\\nstruct Particle \\n{\\n uint mass;\\n vec3 pos;\\n vec3 vel;\\n vec3 force;\\n float density;\\n};\\n\\n//5 bits for shared exponent, 9 bits for each component\\nuint packvec3(vec3 v)\\n{\\n //get the exponent\\n float maxv = max(abs(v.x), max(abs(v.y), abs(v.z)));\\n int exp = clamp(int(ceil(log2(maxv))), -15, 15);\\n float scale = exp2(-float(exp));\\n uvec3 sv = uvec3(round(clamp(v*scale, -1.0, 1.0) * 255.0) + 255.0);\\n uint packed = uint(exp + 15) | (sv.x << 5) | (sv.y << 14) | (sv.z << 23);\\n return packed;\\n}\\n\\nvec3 unpackvec3(uint packed)\\n{\\n int exp = int(packed & 0x1Fu) - 15;\\n vec3 sv = vec3((packed >> 5) & 0x1FFu, (packed >> 14) & 0x1FFu, (packed >> 23) & 0x1FFu);\\n vec3 v = (sv - 255.0) / 255.0;\\n v *= exp2(float(exp));\\n return v;\\n}\\n\\nuint packmasspos(uint mass, vec3 p0)\\n{\\n uvec3 pos0 = uvec3(clamp(p0, 0.0, 1.0) * 255.0);\\n return mass | (pos0.x << 8) | (pos0.y << 16) | (pos0.z << 24);\\n}\\n\\nuint packMassPos(uint mass, vec3 pos)\\n{\\n uvec3 pos0 = uvec3(clamp(pos, 0.0, 1.0) * 255.0); // Assuming pos range [0, 1] in a cell\\n uint data1 = mass | (pos0.x << 8) | (pos0.y << 16) | (pos0.z << 24);\\n return data1;\\n}\\n\\nvoid unpackMassPos(uint packed, out uint mass, out vec3 pos)\\n{\\n mass = packed & 0xFFu;\\n uvec3 pos0 = uvec3((packed >> 8) & 0xFFu, (packed >> 16) & 0xFFu, (packed >> 24) & 0xFFu);\\n pos = vec3(pos0) / 255.0;\\n}\\n\\nvec4 packParticles(Particle p0, Particle p1, vec3 pos)\\n{\\n p0.pos -= pos;\\n p1.pos -= pos;\\n\\n uvec4 data = uvec4(\\n packMassPos(p0.mass, p0.pos),\\n packMassPos(p1.mass, p1.pos),\\n packvec3(p0.vel),\\n packvec3(p1.vel)\\n );\\n\\n return uintBitsToFloat(data);\\n}\\n\\nvoid unpackParticles(vec4 packed, vec3 pos, out Particle p0, out Particle p1)\\n{\\n uvec4 data = floatBitsToUint(packed);\\n \\n unpackMassPos(data.x, p0.mass, p0.pos);\\n unpackMassPos(data.y, p1.mass, p1.pos);\\n\\n p0.pos += pos;\\n p1.pos += pos;\\n\\n p0.vel = unpackvec3(data.z);\\n p1.vel = unpackvec3(data.w);\\n}\\n\\nfloat sdBox( in vec2 p, in vec2 b )\\n{\\n vec2 d = abs(p)-b;\\n return length(max(d,0.0)) + min(max(d.x,d.y),0.0);\\n}\\n\\nint ClosestCluster(Particle p0, Particle p1, Particle incoming)\\n{\\n //first try to choose the particle with significantly smaller mass\\n if(float(p0.mass) < 0.01*float(p1.mass) || float(p1.mass) < 0.01*float(p0.mass))\\n {\\n return p0.mass < p1.mass ? 0 : 1;\\n }\\n\\n //otherwise choose the closest one\\n float d0 = length(p0.pos - incoming.pos);\\n float d1 = length(p1.pos - incoming.pos);\\n return d0 < d1 ? 0 : 1;\\n}\\n\\nvoid BlendParticle(inout Particle p, in Particle incoming)\\n{\\n uint newMass = p.mass + incoming.mass;\\n vec2 weight = vec2(p.mass, incoming.mass) / float(newMass);\\n p.pos = p.pos*weight.x + incoming.pos*weight.y;\\n p.vel = p.vel*weight.x + incoming.vel*weight.y;\\n p.mass = newMass;\\n}\\n\\nvoid Clusterize(inout Particle p0, inout Particle p1, in Particle incoming, vec3 pos, float timeStep)\\n{\\n if(incoming.mass == 0u) return;\\n //move incoming particle forwards\\n incoming.pos += incoming.vel*timeStep;\\n //check if the incoming particle is in the cell\\n if(!all(equal(pos, floor(incoming.pos)))) return;\\n \\n int closest = ClosestCluster(p0, p1, incoming);\\n if(closest == 0)\\n {\\n BlendParticle(p0, incoming);\\n }\\n else\\n {\\n BlendParticle(p1, incoming);\\n }\\n}\\n\\nvoid SplitParticle(inout Particle p1, inout Particle p2)\\n{\\n vec3 pos = p1.pos;\\n uint newMass = p1.mass;\\n p1.mass = newMass/2u;\\n p2.mass = newMass - p1.mass;\\n

@compute.toys