Steiner chain with Julia

SteinerChain.jl

using Meshes
using MeshViz
using GLMakie
using LinearAlgebra
using Makie
using Printf

# image of circle (center, radius) by the inversion 
#   with center c and power k
function iotaCircle(c, k, center, radius)
  r = sqrt(abs(k)) 
  z1 = sign(k) * (center-c)/r 
  D1 = (radius/r)^2 - dot(z1, z1)
  z2 = -z1/D1
  R2 = sqrt(dot(z2, z2) + 1/D1)
  return (center = r*z2 + c, radius = r*R2)
end

# n: vector, the numbers of spheres at each step
# -1 < phi < 1, phi != 0
function steiner(n, phi, shift = 0, Center = [0, 0], radius = 2)
  depth = length(n)
  invphi = 1/phi 
  I = [radius*invphi, 0] + Center
  k = radius*radius * (1 - invphi*invphi)
  m = n[1]
  if depth > 1
    popfirst!(n)
  end
  sine = sin(pi/m)
  Coef = 1 / (1 + sine)
  O1x = 2*radius*invphi
  CRadius = Coef*radius
  CSide = CRadius*sine
  for i in 1:m
    beta = (i + shift)*2*pi/m 
    pti = [CRadius*cos(beta), CRadius*sin(beta)] + Center
    cc = iotaCircle(I, k, pti, CSide)
    center = cc.center - [O1x, 0]
    r = cc.radius
    if depth == 1
      push!(
        spheres, 
        discretize(Meshes.Sphere(Meshes.Point(center[1], center[2], 0), r))
      )
    end
    if depth > 1
      steiner(n, phi, -shift, center, r)
    end
  end
end

# get the spheres
function getSpheres(shift)
  global spheres = Vector{SimpleMesh}(undef, 0)
  steiner([4, 3], 0.2, shift)
  spheres
end

# draw the Steiner chain
function drawSC(nplots)
  shifts = LinRange(0, 1, nplots+1)[1:nplots]
  for i in 1:nplots
    spheres = getSpheres(shifts[i])
    fig, ax, _ = viz( # invisible bounding box
      Meshes.Box(Meshes.Point(-2, -2, -0.5), Meshes.Point(2, 2, 0.5)); 
      alpha = 0
    )
    ax.show_axis = false
    viz!(reduce(merge, spheres); color = :purple)
    scale!(ax.scene, 1.65, 1.65, 1.65)
    Makie.rotate!(ax.scene, Makie.Vec3f(-1, 0, 0), pi/9)
    png = @sprintf "zzpic%03d.png" i
    Makie.save(png, fig)
  end
end