Apollonian circle packing

From JSXGraph Wiki
Revision as of 21:13, 27 April 2009 by A WASSERMANN (talk | contribs)
The printable version is no longer supported and may have rendering errors. Please update your browser bookmarks and please use the default browser print function instead.

The underlying JavaScript code

var brd = JXG.JSXGraph.initBoard('box', {originX: 300, originY: 300, grid:true, unitX: 100, unitY: 100});
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)];
}   

// The outer circle    
p0 = brd.createElement('point',[0,0],{name:'',visible:false});
b0 = -0.5;
c0 = brd.createElement('circle',[p0,Math.abs(1.0/b0)],{strokeWidth:1});
c0.curvature = function(){ return b0;}; // constant

// The left initial inner circle    
a = brd.createElement('line',[[0,0],[2,0]],{straightFirst:false,straightLast:false,visible:false});
p1 = brd.createElement('glider',[1.3,0,a],{name:''});
c1 = brd.createElement('circle',[p1,function(){return 2-p1.X();}],{strokeWidth:1});
c1.curvature = function(){ return 1/(2-p1.X());};

// The right initial inner circle    
p2 = brd.createElement('point',[function(){return p1.X()-2;},0],{name:'',visible:false});
c2 = brd.createElement('circle',[p2,function(){return p1.X();}],{strokeWidth:1});
c2.curvature = function(){ return 1/(p1.X());};

// Now, we construct the fourth initial circle
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];
}

p3 = brd.createElement('point',[thirdCircleX,thirdCircleY],{name:'',visible:false});
c3 = brd.createElement('circle',[p3,thirdCircleRadius],{strokeWidth:1});
c3.curvature = function(){ return 1.0/this.radius;};


// Recursion
// In an Apollonian configuration replace the fourth circle by its
// second possible variant.
otherCirc = function(circs,level) {
    var p,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];
        return 1/b[4];
    }
    p = brd.createElement('point',[fx,fy],{name:'',visible:false});
    c = brd.createElement('circle',[p,fr],{strokeWidth:1,fillColor:HSV2RGB(360/(level),0.8,0.8),fillOpacity:0.5});
    c.curvature = function(){ return 1/this.radius;};
    
    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;
}

//-------------------------------------------------------
brd.suspendUpdate();
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);
brd.unsuspendUpdate();

References