Share JSXGraph: example "Apollonian circle packing"

JSXGraph
Share JSXGraph: example "Apollonian circle packing"
This website is a beta version. The official release will be in **2023**.

Apollonian circle packing

__References__ - [https://www.ams.org/featurecolumn/archive/kissing.html](https://www.ams.org/featurecolumn/archive/kissing.html) - Jeffrey C. Lagarias, Colin L. Mallows, Allan R. Wilks: Beyond the Descartes circle theorem - [https://en.wikipedia.org/wiki/Apollonian_gasket](https://en.wikipedia.org/wiki/Apollonian_gasket) - Weisstein, Eric W. ["Apollonian Gasket."](https://mathworld.wolfram.com/ApollonianGasket.html) From MathWorld - A Wolfram Web Resource - Weisstein, Eric W. ["Soddy Circles."](https://mathworld.wolfram.com/SoddyCircles.html) From MathWorld - A Wolfram Web Resource
// Define the id of your board in BOARDID

const board = JXG.JSXGraph.initBoard(BOARDID, {
    boundingbox: [-2, 2, 2, -2],
    keepaspectratio: true
});

var b0, c0, c1, c2, c3,
    a, p1,
    solveQ2, thirdCircleX, thirdCircleY, thirdCircleRadius,
    otherCirc;

solveQ2 = function(x1, x2, x3, off) {
    var a, b, c, d;
    a = 0.5;
    b = -(x1 + x2 + x3);
    c = x1 * x1 + x2 * x2 + x3 * x3 - 0.5 * (x1 + x2 + x3) * (x1 + x2 + x3) - off;
    d = b * b - 4 * a * c;
    if (Math.abs(d) < 0.00000001) d = 0.0;
    return [(-b + Math.sqrt(d)) / (2.0 * a), (-b - Math.sqrt(d)) / (2.0 * a)];
}

a = board.create('segment', [
    [0, 0],
    [2, 0]
], {
    visible: false
});

p1 = board.create('glider', [1.3, 0, a], {
    name: 'Drag me'
});
b0 = -0.5;

c0 = board.create('circle', [
    [0, 0], Math.abs(1.0 / b0)
], {
    strokeWidth: 1
});
c1 = board.create('circle', [p1, function() {
    return 2 - p1.X();
}], {
    strokeWidth: 1
});

c2 = board.create('circle', [
    [function() {
        return p1.X() - 2;
    }, 0],
    function() {
        return p1.X();
    }
], {
    strokeWidth: 1
});

// Constant curvature
c0.curvature = function() {
    return b0;
};
c1.curvature = function() {
    return 1 / (2 - p1.X());
};
c2.curvature = function() {
    return 1 / (p1.X());
};

thirdCircleX = function() {
    var b0, b1, b2, x0, x1, x2, b3, bx3;
    b0 = c0.curvature();
    b1 = c1.curvature();
    b2 = c2.curvature();
    x0 = c0.midpoint.X();
    x1 = c1.midpoint.X();
    x2 = c2.midpoint.X();

    b3 = solveQ2(b0, b1, b2, 0);
    bx3 = solveQ2(b0 * x0, b1 * x1, b2 * x2, 2);
    return bx3[0] / b3[0];
}
thirdCircleY = function() {
    var b0, b1, b2, y0, y1, y2, b3, by3;
    b0 = c0.curvature();
    b1 = c1.curvature();
    b2 = c2.curvature();
    y0 = c0.midpoint.Y();
    y1 = c1.midpoint.Y();
    y2 = c2.midpoint.Y();

    b3 = solveQ2(b0, b1, b2, 0);
    by3 = solveQ2(b0 * y0, b1 * y1, b2 * y2, 2);
    return by3[0] / b3[0];
}
thirdCircleRadius = function() {
    var b0, b1, b2, b3, bx3, by3;
    b0 = c0.curvature();
    b1 = c1.curvature();
    b2 = c2.curvature();
    b3 = solveQ2(b0, b1, b2, 0);
    return 1.0 / b3[0];
}

c3 = board.create('circle', [
    [thirdCircleX, thirdCircleY], thirdCircleRadius
], {
    strokeWidth: 1
});
c3.curvature = function() {
    return 1.0 / this.radius;
};

otherCirc = function(circs, level) {
    var c, fx, fy, fr;
    if (level <= 0) return;

    fx = function() {
        var b, x, i;
        b = [];
        x = [];
        for (i = 0; i < 4; i++) {
            b[i] = circs[i].curvature();
            x[i] = circs[i].midpoint.X();
        }

        b[4] = 2 * (b[0] + b[1] + b[2]) - b[3];
        x[4] = (2 * (b[0] * x[0] + b[1] * x[1] + b[2] * x[2]) - b[3] * x[3]) / b[4];
        return x[4];
    }
    fy = function() {
        var b, y, i;
        b = [];
        y = [];
        for (i = 0; i < 4; i++) {
            b[i] = circs[i].curvature();
            y[i] = circs[i].midpoint.Y();
        }

        b[4] = 2 * (b[0] + b[1] + b[2]) - b[3];
        y[4] = (2 * (b[0] * y[0] + b[1] * y[1] + b[2] * y[2]) - b[3] * y[3]) / b[4];
        return y[4];
    }
    fr = function() {
        var b, i;
        b = [];
        for (i = 0; i < 4; i++) {
            b[i] = circs[i].curvature();
        }
        b[4] = 2 * (b[0] + b[1] + b[2]) - b[3];
        if (isNaN(b[4])) {
            return 1000.0;
        } else {
            return 1 / b[4];
        }
    }
    c = board.create('circle', [
        [fx, fy], fr
    ], {
        strokeWidth: 1,
        name: '',
        fillColor: JXG.hsv2rgb(50 * level, 0.8, 0.8),
        highlightFillColor: JXG.hsv2rgb(50 * level, 0.5, 0.8),
        fillOpacity: 0.5,
        highlightFillOpacity: 0.5
    });
    c.curvature = function() {
        return 1 / this.radius;
    };

    // Recursion
    otherCirc([circs[0], circs[1], c, circs[2]], level - 1);
    otherCirc([circs[0], circs[2], c, circs[1]], level - 1);
    otherCirc([circs[1], circs[2], c, circs[0]], level - 1);
    return c;
}

//-------------------------------------------------------

var level = 4;
otherCirc([c0, c1, c2, c3], level);
otherCirc([c3, c1, c2, c0], level);
otherCirc([c0, c2, c3, c1], level);
otherCirc([c0, c1, c3, c2], level);