Now we will do something a bit more complicated. We will take a
triangle, fill each of its three Malfatti circles with an Apollonian
gasket, and fill the rest of the triangle with tangent circles. In fact,
tangent spheres: we will draw the result in 3D, this will be more
pretty.
toCplx <- function(M) {
M[1L] + 1i * M[2L]
}
fromCplx <- function(z) {
c(Re(z), Im(z))
}
distance <- function(A, B) {
sqrt(c(crossprod(B - A)))
}
innerSoddyRadius <- function(r1, r2, r3) {
1 / (1/r1 + 1/r2 + 1/r3 + 2 * sqrt(1/r1/r2 + 1/r2/r3 + 1/r3/r1))
}
innerSoddyCircle <- function(c1, c2, c3, ...) {
radius <- innerSoddyRadius(c1$radius, c3$radius, c3$radius)
center <- Triangle$new(c1$center, c2$center, c3$center)$equalDetourPoint()
c123 <- Circle$new(center, radius)
drawSphere(c123, ...)
list(
list(type = "ccc", c1 = c123, c2 = c1, c3 = c2),
list(type = "ccc", c1 = c123, c2 = c2, c3 = c3),
list(type = "ccc", c1 = c123, c2 = c1, c3 = c3)
)
}
side.circle.circle <- function(A, B, cA, cB, ...) {
if(A[2L] > B[2L]){
return(side.circle.circle(B, A, cB, cA, ...))
}
rA <- cA$radius
rB <- cB$radius
zoA <- toCplx(cA$center)
zoB <- toCplx(cB$center)
zB <- toCplx(A)
alpha <- acos((B[1L] - A[1L]) / distance(A, B))
zX1 <- exp(-1i * alpha) * (zoA - zB)
zX2 <- exp(-1i * alpha) * (zoB - zB)
soddyR <- innerSoddyRadius(rA, rB, Inf)
if(Re(zX1) < Re(zX2)) {
Y <- (2 * rA * sqrt(rB) / (sqrt(rA) + sqrt(rB)) + Re(zX1)) +
sign(Im(zX1)) * 1i * soddyR
} else {
Y <- (2 * rB * sqrt(rA) / (sqrt(rA) + sqrt(rB)) + Re(zX2)) +
sign(Im(zX1)) * 1i * soddyR
}
M <- fromCplx(Y * exp(1i * alpha) + zB)
cAB <- Circle$new(M, soddyR)
drawSphere(cAB, ...)
list(
list(type = "ccc", c1 = cAB, c2 = cA, c3 = cB),
list(type = "ccl", cA = cA, cB = cAB, A = A, B = B),
list(type = "ccl", cA = cAB, cB = cB, A = A, B = B)
)
}
side.side.circle <- function(A, B, C, circle, ...) {
zA <- toCplx(A)
zO <- toCplx(circle$center)
vec <- zA - zO
P <- fromCplx(zO + circle$radius * vec / Mod(vec))
OP <- P - circle$center
onTangent <- P + c(-OP[2L], OP[1L])
L1 <- Line$new(P, onTangent)
P1 <- intersectionLineLine(L1, Line$new(A, C))
P2 <- intersectionLineLine(L1, Line$new(A, B))
incircle <- Triangle$new(A, P1, P2)$incircle()
drawSphere(incircle, ...)
list(
list(type = "cll", A = A, B = B, C = C, circle = incircle),
list(type = "ccl", cA = circle, cB = incircle, A = A, B = B),
list(type = "ccl", cA = circle, cB = incircle, A = A, B = C)
)
}
Newholes <- function(holes, color) {
newholes <- list()
for(i in 1L:3L) {
hole <- holes[[i]]
holeType <- hole[["type"]]
if(holeType == "ccc") {
x <- with(hole, innerSoddyCircle(c1, c2, c3, color = color))
} else if(holeType == "ccl") {
x <- with(hole, side.circle.circle(A, B, cA, cB, color = color))
} else if (holeType == "cll") {
x <- with(hole, side.side.circle(A, B, C, circle, color = color))
}
newholes <- c(newholes, list(x))
}
newholes
}
MalfattiCircles <- function(A, B, C) {
Triangle$new(A, B, C)$MalfattiCircles()
}
drawTriangularGasket <- function(mcircles, A, B, C, colors, depth) {
C1 <- mcircles[[1L]]
C2 <- mcircles[[2L]]
C3 <- mcircles[[3L]]
triangles3d(cbind(rbind(A, B, C), 0), col = "yellow", alpha = 0.2)
holes <- list(
side.circle.circle(A, B, C1, C2, color = colors[1L]),
side.circle.circle(B, C, C2, C3, color = colors[1L]),
side.circle.circle(C, A, C3, C1, color = colors[1L]),
innerSoddyCircle(C1, C2, C3, color = colors[1L]),
side.side.circle(A, B, C, C1, color = colors[1L]),
side.side.circle(B, A, C, C2, color = colors[1L]),
side.side.circle(C, A, B, C3, color = colors[1L])
)
for(d in 1L:depth) {
n_holes <- length(holes)
Holes <- list()
for(i in 1L:n_holes) {
Holes <- append(Holes, Newholes(holes[[i]], colors[d + 1L]))
}
holes <- do.call(list, Holes)
}
}
drawCircularGasket <- function(c0, n, phi, shift, depth, colors) {
ApollonianCircles <- ApollonianGasket(c0, n, phi, shift, depth)
for(i in 1:depth) {
for(circ in ApollonianCircles[[i]]){
drawSphere(circ, color = colors[i])
}
}
}
library(viridisLite)
A <- c(-5, -4)
B <- c(5, -2)
C <- c(0, 6)
mcircles <- MalfattiCircles(A, B, C)
depth <- 3L
colors <- viridis(depth + 1L)
n1 <- 3L
n2 <- 4L
n3 <- 5L
depth2 <- 3L
phi1 <- 0.2
phi2 <- 0.3
phi3 <- 0.4
shift <- 0
colors2 <- plasma(depth2)
library(rgl)
open3d(windowRect = c(50, 50, 562, 562), zoom = 0.9)
bg3d(rgb(54, 57, 64, maxColorValue = 255))
drawTriangularGasket(mcircles, A, B, C, colors, depth)
drawCircularGasket(mcircles[[1L]], n1, phi1, shift, depth2, colors2)
drawCircularGasket(mcircles[[2L]], n2, phi2, shift, depth2, colors2)
drawCircularGasket(mcircles[[3L]], n3, phi3, shift, depth2, colors2)
As an exercise, you can add some animation to this picture, by
animating the three circular gaskets.
