ussianGrad(vec3 dx, float d)\\n{\\n float norm = 1.0/(cub(d)*sqrt(cub(TWO_PI)));\\n float r = length(dx);\\n return - (norm/sqr(d)) * exp(-0.5*sqr(r/d)) * dx;\\n}\\n\\n#define GD(x, R) Gaussian(length(x),R)\\n#define GGRAD(x, R) GaussianGrad(x, R)\\n#define GS(x) exp(-dot(x,x))\\n#define DIR(phi) vec2(cos(TWO_PI*phi),sin(TWO_PI*phi))\\n\\n#define loop(i,x) for(int i = 0; i < x; i++)\\n#define range(i,a,b) for(int i = a; i <= b; i++)\\n\\n#define pixel(a, p, s) texture(a, (p+0.5)/vec2(s))\\n\\nvec4 voxel(int ch, vec3 p3d)\\n{\\n return pixel(ch, dim2from3(p3d), size2d);\\n}\\n\\n//trilinearSample interpolation = linear interp between layers\\nvec4 trilinearSample(int ch, vec3 p3d)\\n{\\n return mix(voxel(ch, vec3(p3d.xy, floor(p3d.z))),voxel(ch, vec3(p3d.xy, ceil(p3d.z))), fract(p3d.z));\\n}\\n\\nfloat sdBox( vec3 p, vec3 b )\\n{\\n vec3 d = abs(p) - b;\\n return min(max(d.x,max(d.y,d.z)),0.0) + length(max(d,0.0));\\n}\\n\\nvec2 hash21(float p)\\n{\\n\\tvec3 p3 = fract(vec3(p) * vec3(.1031, .1030, .0973));\\n\\tp3 += dot(p3, p3.yzx + 33.33);\\n return fract((p3.xx+p3.yz)*p3.zy);\\n}\\n\\nvec2 hash23(vec3 p3)\\n{\\n\\tp3 = fract(p3 * vec3(.1031, .1030, .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 close

@compute.toys