JSFiddle - React, Tailwind, and code Playground

HTML

<canvas id="c"></canvas>

JavaScript

b = document.body;
c = document.getElementsByTagName('canvas')[0];
a = c.getContext('2d');

with(Math)S=sin,C=cos,R=random,Q=sqrt,P=PI,M=max,N=min;
for(Z in a)a[Z[0]+(Z[6]||Z[2])]=a[Z];
console.log(a);

var blobs = {};
var nb = 0;
var np = 20;
var inc = 2 * P / np;
setInterval(function() {
    W=c.width=1000;
    H=c.height=500;
    if (R() < .5 && nb < 200) {
        var blob = blobs[nb++] = [];
        blob.cx = cx = R() * W;
        blob.cy = cy = R() * H;
        var r = blob.r = 10+R()*30;
        for (var j = 0; j < np; j++) {
            var p = blob[j] = {};
            p.xp = p.x = cx + r * C(j * inc);
            p.yp = p.y = cy + r * S(j * inc);
            p.dx = p.dy = p.nx = p.ny = 0;
        }
        blob[np] = blob[0];
        blob.color = 'hsl('+R()*256+',100%,50%)';
        blob.a = 2*P*r*r;
    }
    // How to correctly handle blob size change
    // You have to change r so that the collision handling works,
    // and a so the volume preservation works
    
    // Collide blobs
    for (var i = 0;i<nb;i++) {
        var blob = blobs[i];
        for (var j=i+1;j<nb;j++) {
            var coll = blobs[j];
            var bound = (blob.r+coll.r)*1.5;
            var vx = coll.cx - blob.cx;
            if (vx > -bound && vx < bound) {
                var vy = coll.cy - blob.cy;
                if (vy > -bound && vy < bound) {
                    var len = Q(vx * vx + vy * vy);
                    var dx = vx / len;
                    var dy = vy / len;
                    var rsum = blob.r+coll.r;
                    var l1 = blob.r/rsum*len;
                    var l2 = coll.r/rsum*len;
                    for (k = 0; k < np; k++) {
                        var p1 = blob[k];
                        var p2 = coll[k];

                        var dp = M(0, (p1.x - blob.cx) * dx + (p1.y - blob.cy) * dy - l1);
                        p1.x -= dp * dx;
                        p1.y -= dp * dy;
                        dp = N(0, (p2.x - coll.cx) * dx +...