Apollonian circles playground

apollonian_circles.jscad

"use strict"
const jscad = require('@jscad/modeling')
const math = require('mathjs')

const { extrudeLinear } = jscad.extrusions
const { hullChain } = jscad.hulls
const { circle } = jscad.primitives
const { vectorText } = jscad.text
const { scale, translate } = jscad.transforms

const th=1
const sc=0.03
const seg=360
const f=100


function soddyCenters(A, B, r1, r2) {
    const k1 =  1 / r1;
    const k2 =  1 / r2;

    const Q = k1*k2 -(k2 + k1);
    const sqrtQ = Math.sqrt(Q);

    const k3a = k1 + k2 - 1 + 2*sqrtQ;
    const k3b = k1 + k2 - 1 - 2*sqrtQ;

    const z1 = math.complex(A[0], A[1]);
    const z2 = math.complex(B[0], B[1]);

    function computeCenter(k3, sign) {
        const T1 = math.multiply(z1, k1);
        const T2 = math.multiply(z2, k2);

        const S = math.add(T1, T2);

        const c12 = math.multiply(z1, z2);

        const under = math.multiply(c12, k1*k2)

        const mag = Math.sqrt(under.re*under.re + under.im*under.im);
        const phase = Math.atan2(under.im, under.re);
        const sqrtUnder = math.complex(
            Math.sqrt(mag) * Math.cos(phase/2),
            Math.sqrt(mag) * Math.sin(phase/2)
        );

        const numerator = math.add(S, math.multiply(sqrtUnder, 2*sign));
        return [numerator.re / k3, numerator.im / k3, k3] 
    }

    const Cplus  = computeCenter(k3a, +1);
    const Cminus = computeCenter(k3b, -1);

    return { Cplus, Cminus };
}

function outerCenter(A, B, r1, r2) {

    function side(A, B, P) {
        return (B[0] - A[0]) * (P[1] - A[1]) - (B[1] - A[1]) * (P[0] - A[0]);
    }

    const { Cplus, Cminus } = soddyCenters(A, B, r1, r2);

    const sOuter = side(A, B, [0,0]);
    const sPlus  = side(A, B, Cplus);
    const sMinus = side(A, B, Cminus);

    if (Math.sign(sPlus) !== Math.sign(sOuter)) {
        return Cplus;
    } else {
        return Cminus;
    }
}

// const gap = chooseGapCircle([2/3, 1/2 ], [2/3, 0 ], 1/6, 1/3);

function C(c,notext=false) {
  return [
    notext?[]:buildFlatText([f*(c[0]-sc*25/c[2]),f*(c[1]-0.3/c[2])], ""+c[2], th, f*sc/c[2]),
    circle({segments: seg, center: [f*c[0],f*c[1]], radius: f/Math.abs(c[2])})
  ]
}

const issquare = (n) => { var i=isqrt(n); return n === i*i }
const isqrt = (n) => { return Math.floor(Math.sqrt(n)) }

function nxt(B1,B2,B3) {
  var a,b,c
  [a,b,c]=[B1[2],B2[2],B3[2]]
  console.assert(issquare(a*b+a*c+b*c))
  return a+b+c+2*isqrt(a*b+a*c+b*c)
}

function Center(B1,B2,B3,k4) {
  var a=r1+r2
  var b=r1+r3
  var c=r2+r3
  var k1=B1[2]
  var k2=B2[2]
  var k3=B3[2]
  var r1=1/k1
  var r2=1/k2
  var r3=1/k3
  var z1 = math.complex(B1[0],B1[1])
  var z2 = math.complex(B2[0],B2[1])
  var z3 = math.complex(B3[0],B3[1])
  // z4=(k1*z1+k2*z2+k3*z3+2*sqrt(k1*k2*z1*z2+k2*k3*z2*z3+k3*k1*z3*z1))/k4
  var z4 =
    math.multiply(
      math.add(
        math.add(
          math.multiply(k1,z1),
          math.add(
            math.multiply(k2,z2),
            math.multiply(k3,z3)
          )
        ),
        math.multiply(2,
          math.sqrt(
            math.add(
              math.add(
                math.multiply(
                  k1*k2,
                  math.multiply(z1,z2)
                ),
                math.multiply(k2*k3,
                  math.multiply(z2,z3)
                )
              ),
              math.multiply(k3*k1,
                math.multiply(z3,z1)
              )
            )
          )
        )
      ),
      1/k4
    )

  return [z4.re, z4.im]
}

// k3*z3 = k1*z1 + k2*z2 + k0*z0 ± 2 sqrt(k1*k2*z1*z2 + k2*k0*z2*z0 + k0*k1*z0*z1)
// z0=(0,0)
// k3*z3 = k1*z1 + k2*z2 ± 2 sqrt(k1*k2*z1*z2)

function i(B1,B2,B3) {
  var n=nxt(B1,B2,B3)
  var c=Center(B1,B2,B3,n).concat([n])
  console.log(c)
  return C(c)
}

function o(B1,B2) {
  const out = outerCenter(B1, B2, 1/B1[2], 1/B2[2]);
  console.log(out)
  return C(out)
}

function r(B1,B2,B3) {
  return [
    o(B1,B3),
    o(B2,B3),
    i(B1,B2,B3)
  ]
}

function T(B1,B2,B3,l) {
  if (l===0)  return []
  
  const A13 = outerCenter(B1, B3, 1/B1[2], 1/B3[2])
  const A23 = outerCenter(B2, B3, 1/B2[2], 1/B3[2])
  const n=nxt(B1,B2,B3)
  const A123=Center(B1,B2,B3,n).concat([n])

  if (l===1)  return [C(A13), C(A23), C(A123)]

  return [C(A13), C(A23), C(A123),
          T(B1,B3,A13,l-1),T(B2,B3,A23,l-1)]
}

function F(B1,B2,B3,l) {
  if (l===0)  return []
  
  const n=nxt(B1,B2,B3)
  const A123=Center(B1,B2,B3,n).concat([n])

  if (l===1)  return [C(A123)]

  return [C(A123),T(B1,B3,A13,l-1),T(B2,B3,A23,l-1)]
}

function main() {
  var A=[[0,0,-1],[0,1/2,2],[0,-1/2,2],[2/3,0,3], [2/3,0.5,6]]
    console.log(A[0])
    console.log(A[1])
    console.log(A[2])
//    console.log(A[3])
//  const out = outerCenter([2/3, 1/2 ], [2/3, 0 ], 1/6, 1/3);
  return [
    C(A[0], true), buildFlatText([-5,f], "-1", th, f*sc),
    C(A[1]),
    C(A[2]),
    C(A[3]),


    T(A[1],A[2],A[3],6),
//    o(A[1],A[3],[1,1]),
//    i(A[1],A[2],A[3]),
//    o(A[1],A[3],[0,0]),
//    i(A[1],A[3],A[4]),
//    o(A[1],A[4],[0,0]),
//    o(A[3],A[4],[0,0]),
//    o(A[2],A[3],[0,0]),
  ]
}

const buildFlatText = (pos, message, characterLineWidth, s) => {
  if (message === undefined || message.length === 0) return []

  const lineSegments = []
  vectorText({ x: 0, y: 0, input: message }).forEach((segmentPoints) => {
    const corners = segmentPoints.map((point) => 
      translate(point, circle({ radius: characterLineWidth / 2})))
    lineSegments.push(hullChain(corners))
  })
  return translate(pos, scale([s,s,s],extrudeLinear({ height: 0.00001 }, lineSegments)))
}

module.exports = { main }

Initial version.

• compute bend=1/radius of inner (to 2,2,3 bend circles) circle as 15
• compute center of that circle so that it "kisses" the three other circles
• lots of mathjs complex number computations

jscad.app share link:
https://jscad.app/#https://gist.githubusercontent.com/Hermann-SW/374019946acdb2d9399f7619247ae4c3/raw/fa33731936b415cebc6a25968d7ff75f1aa36b82/apollonian_circles.jscad

recursion is not completely correct, and obviously some integers are not integer ;-)
https://jscad.app/#https://gist.githubusercontent.com/Hermann-SW/374019946acdb2d9399f7619247ae4c3/raw/ab7a65d6ecd17bf27ebbd6e5d382cd992393fa8c/apollonian_circles.jscad
The displayed integers correspond to unit fraction radii: