knitr::opts_chunk$set( collapse = TRUE, comment = "#>" )  library(PlaneGeometry)  # Exeter point The Exeter point is defined as follows on Wikipedia. Let$ABC$be any given triangle. Let the medians through the vertices$A$,$B$,$C$meet the circumcircle of triangle$ABC$at$A'$,$B'$and$C'$respectively. Let$DEF$be the triangle formed by the tangents at$A$,$B$, and$C$to the circumcircle of triangle$ABC$. (Let$D$be the vertex opposite to the side formed by the tangent at the vertex$A$, let$E$be the vertex opposite to the side formed by the tangent at the vertex$B$, and let$F$be the vertex opposite to the side formed by the tangent at the vertex$C$.) The lines through$DA'$,$EB'$and$FC'$are concurrent. The point of concurrence is the Exeter point of triangle$ABC$. Let's construct it with the PlaneGeometry package. We do not need to construct the triangle$DEF$: it is the tangential triangle of$ABC$, and is provided by the tangentialTriangle method of the R6 class Triangle. A <- c(0,2); B <- c(5,4); C <- c(5,-1) t <- Triangle$new(A, B, C)
circumcircle <- t$circumcircle() centroid <- t$centroid()
medianA <- Line$new(A, centroid) medianB <- Line$new(B, centroid)
medianC <- Line$new(C, centroid) Aprime <- intersectionCircleLine(circumcircle, medianA)[] Bprime <- intersectionCircleLine(circumcircle, medianB)[] Cprime <- intersectionCircleLine(circumcircle, medianC)[] DEF <- t$tangentialTriangle()
lineDAprime <- Line$new(DEF$A, Aprime)
lineEBprime <- Line$new(DEF$B, Bprime)
lineFCprime <- Line$new(DEF$C, Cprime)
( ExeterPoint <- intersectionLineLine(lineDAprime, lineEBprime) )
# check whether the Exeter point is also on (FC')
lineFCprime$includes(ExeterPoint)  Let's draw a figure now. opar <- par(mar = c(0,0,0,0)) plot(NULL, asp = 1, xlim = c(-2,9), ylim = c(-6,7), xlab = NA, ylab = NA, axes = FALSE) draw(t, lwd = 2, col = "black") draw(circumcircle, lwd = 2, border = "cyan") draw(Triangle$new(Aprime,Bprime,Cprime), lwd = 2, col = "green")
draw(DEF, lwd = 2, col = "blue")
draw(Line$new(ExeterPoint, DEF$A, FALSE, FALSE), lwd = 2, col = "red")
draw(Line$new(ExeterPoint, DEF$B, FALSE, FALSE), lwd = 2, col = "red")
draw(Line$new(ExeterPoint, DEF$C, FALSE, FALSE), lwd = 2, col = "red")
points(rbind(ExeterPoint), pch = 19, col = "red")
par(opar)


# Circles tangent to three circles

Let $\mathcal{C}_1$, $\mathcal{C}_2$ and $\mathcal{C}_3$ be three circles with respective radii $r_1$, $r_2$ and $r_3$ such that $r_3 < r_1$ and $r_3 < r_2$. How to construct some circles simultaneously tangent to these three circles?

C1 <- Circle$new(c(0,0), 2) C2 <- Circle$new(c(5,5), 3)
C3 <- Circle$new(c(6,-2), 1) # inversion swapping C1 and C3 with positive power iota1 <- inversionSwappingTwoCircles(C1, C3, positive = TRUE) # inversion swapping C2 and C3 with positive power iota2 <- inversionSwappingTwoCircles(C2, C3, positive = TRUE) # take an arbitrary point on C3 M <- C3$pointFromAngle(0)
# invert it with iota1 and iota2
M1 <- iota1$invert(M); M2 <- iota2$invert(M)
# take the circle C passing through M, M1, M2
C <- Triangle$new(M,M1,M2)$circumcircle()
# take the line passing through the two inversion poles
cl <- Line$new(iota1$pole, iota2$pole) # take the radical axis of C and C3 L <- C$radicalAxis(C3)
# let H bet the intersection of these two lines
H <- intersectionLineLine(L, cl)
# take the circle Cp with diameter [HO3]
O3 <- C3$center Cp <- CircleAB(H, O3) # get the two intersection points T0 and T1 of C3 with Cp T0_and_T1 <- intersectionCircleCircle(C3, Cp) T0 <- T0_and_T1[[1L]]; T1 <- T0_and_T1[[2L]] # invert T0 with respect to the two inversions T0p <- iota1$invert(T0); T0pp <- iota2$invert(T0) # the circle passing through T0 and its two images is a solution Csolution0 <- Triangle$new(T0, T0p, T0pp)$circumcircle() # invert T1 with respect to the two inversions T1p <- iota1$invert(T1); T1pp <- iota2$invert(T1) # the circle passing through T1 and its two images is another solution Csolution1 <- Triangle$new(T1, T1p, T1pp)$circumcircle()  opar <- par(mar = c(0,0,0,0)) plot(NULL, asp = 1, xlim = c(-4,9), ylim = c(-4,9), xlab = NA, ylab = NA, axes = FALSE) draw(C1, col = "yellow", border = "red") draw(C2, col = "yellow", border = "red") draw(C3, col = "yellow", border = "red") draw(Csolution0, lwd = 2, border = "blue") draw(Csolution1, lwd = 2, border = "blue") par(opar)  # Apollonius circle of a triangle There are several circles called "Apollonius circle". We take the one defined as follows, with respect to a reference triangle: the circle which touches all three excircles of the reference triangle and encompasses them. It can be constructed as the inversive image of the nine-point circle with respect to the circle orthogonal to the excircles of the reference triangle. This inversion can be obtained in PlaneGeometry with the function inversionFixingThreeCircles. # reference triangle t <- Triangle$new(c(0,0), c(5,3), c(3,-1))
# nine-point circle
npc <- t$orthicTriangle()$circumcircle()
# excircles
excircles <- t$excircles() # inversion with respect to the circle orthogonal to the excircles iota <- inversionFixingThreeCircles(excircles$A, excircles$B, excircles$C)
# Apollonius circle
ApolloniusCircle <- iota$invertCircle(npc)  Let's do a figure: opar <- par(mar = c(0,0,0,0)) plot(NULL, asp = 1, xlim = c(-10,14), ylim = c(-5, 18), xlab = NA, ylab = NA, axes = FALSE) draw(t, lwd = 2) draw(excircles$A, lwd = 2, border = "blue")
draw(excircles$B, lwd = 2, border = "blue") draw(excircles$C, lwd = 2, border = "blue")
draw(ApolloniusCircle, lwd = 2, border = "red")
par(opar)


The radius of the Apollonius circle is $\frac{r^2+s^2}{4r}$ where $r$ is the inradius of the triangle and $s$ its semiperimeter. Let's check this point:

inradius <- t$inradius() semiperimeter <- sum(t$edges()) / 2
ApolloniusCircle$radius  # Filling the lapping area of two circles Let two circles intersecting at two points. How to fill the lapping area of the two circles? O1 <- c(2,5); circ1 <- Circle$new(O1, 2)
O2 <- c(4,4); circ2 <- Circle$new(O2, 3) opar <- par(mar = c(0,0,0,0)) plot(NULL, asp = 1, xlim = c(0,8), ylim = c(0,8), xlab = NA, ylab = NA) draw(circ1, border = "purple", lwd = 2) draw(circ2, border = "forestgreen", lwd = 2) intersections <- intersectionCircleCircle(circ1, circ2) A <- intersections[]; B <- intersections[] points(rbind(A,B), pch = 19, col = c("red", "blue")) theta1 <- Arg((A-O1) + 1i*(A-O1)) theta2 <- Arg((B-O1) + 1i*(B-O1)) path1 <- Arc$new(O1, circ1$radius, theta1, theta2, FALSE)$path()

theta1 <- Arg((A-O2) + 1i*(A-O2))
theta2 <- Arg((B-O2) + 1i*(B-O2))
path2 <- Arc$new(O2, circ2$radius, theta2, theta1, FALSE)$path() polypath(rbind(path1,path2), col = "yellow") par(opar)  # Hyperbolic tessellation In the help page of the Circle R6 class (?Circle), we show how to draw a hyperbolic triangle with the help of the method orthogonalThroughTwoPointsOnCircle(). Here we will use this method to draw a hyperbolic tessellation. tessellation <- function(depth, Thetas0, colors){ stopifnot( depth >= 3, is.numeric(Thetas0), length(Thetas0) == 3L, is.character(colors), length(colors) >= depth ) circ <- Circle$new(c(0,0), 3)

arcs <- lapply(seq_along(Thetas0), function(i){
ip1 <- ifelse(i == length(Thetas0), 1L, i+1L)
circ$orthogonalThroughTwoPointsOnCircle(Thetas0[i], Thetas0[ip1], arc = TRUE) }) inversions <- lapply(arcs, function(arc){ Inversion$new(arc$center, arc$radius^2)
})

Ms <- vector("list", depth)

Ms[[1L]] <- lapply(Thetas0, function(theta) c(cos(theta), sin(theta)))

Ms[[2L]] <- vector("list", 3L)
for(i in 1L:3L){
im1 <- ifelse(i == 1L, 3L, i-1L)
M <- inversions[[i]]$invert(Ms[[1L]][[im1]]) attr(M, "iota") <- i Ms[[2L]][[i]] <- M } for(d in 3L:depth){ n1 <- length(Ms[[d-1L]]) n2 <- 2L*n1 Ms[[d]] <- vector("list", n2) k <- 0L while(k < n2){ for(j in 1L:n1){ M <- Ms[[d-1L]][[j]] for(i in 1L:3L){ if(i != attr(M, "iota")){ k <- k + 1L newM <- inversions[[i]]$invert(M)
attr(newM, "iota") <- i
Ms[[d]][[k]] <- newM
}
}
}
}
}

# plot ####
opar <- par(mar = c(0,0,0,0), bg = "black")
plot(NULL, asp = 1, xlim = c(-4,4), ylim = c(-4,4),
xlab = NA, ylab = NA, axes = FALSE)
draw(circ, border = "white")
invisible(lapply(arcs, draw, col = colors[1L], lwd = 2))

Thetas <- lapply(
rapply(Ms, function(M){
Arg(M[1L] + 1i*M[2L])
}, how="replace"),
unlist)

for(d in 2L:depth){
thetas <- sort(unlist(Thetas[1L:d]))
for(i in 1L:length(thetas)){
ip1 <- ifelse(i == length(thetas), 1L, i+1L)
arc <- circ$orthogonalThroughTwoPointsOnCircle(thetas[i], thetas[ip1], arc = TRUE) draw(arc, lwd = 2, col = colors[d]) } } par(opar) invisible() }  tessellation( depth = 5L, Thetas0 = c(0, 2, 3.8), colors = viridisLite::viridis(5) )  Here is a version which allows to fill the hyperbolic triangles: tessellation2 <- function(depth, Thetas0, colors){ stopifnot( depth >= 3, is.numeric(Thetas0), length(Thetas0) == 3L, is.character(colors), length(colors)-1L >= depth ) circ <- Circle$new(c(0,0), 3)

arcs <- lapply(seq_along(Thetas0), function(i){
ip1 <- ifelse(i == length(Thetas0), 1L, i+1L)
circ$orthogonalThroughTwoPointsOnCircle(Thetas0[i], Thetas0[ip1], arc = TRUE) }) inversions <- lapply(arcs, function(arc){ Inversion$new(arc$center, arc$radius^2)
})

Ms <- vector("list", depth)

Ms[[1L]] <- lapply(Thetas0, function(theta) c(cos(theta), sin(theta)))

Ms[[2L]] <- vector("list", 3L)
for(i in 1L:3L){
im1 <- ifelse(i == 1L, 3L, i-1L)
M <- inversions[[i]]$invert(Ms[[1L]][[im1]]) attr(M, "iota") <- i Ms[[2L]][[i]] <- M } for(d in 3L:depth){ n1 <- length(Ms[[d-1L]]) n2 <- 2L*n1 Ms[[d]] <- vector("list", n2) k <- 0L while(k < n2){ for(j in 1L:n1){ M <- Ms[[d-1L]][[j]] for(i in 1L:3L){ if(i != attr(M, "iota")){ k <- k + 1L newM <- inversions[[i]]$invert(M)
attr(newM, "iota") <- i
Ms[[d]][[k]] <- newM
}
}
}
}
}

# plot ####
opar <- par(mar = c(0,0,0,0), bg = "black")
plot(NULL, asp = 1, xlim = c(-4,4), ylim = c(-4,4),
xlab = NA, ylab = NA, axes = FALSE)

path <- do.call(rbind, lapply(rev(arcs), function(arc) arc$path())) polypath(path, col = colors[1L]) invisible(lapply(arcs, function(arc){ path1 <- arc$path()
B <- arc$startingPoint() A <- arc$endingPoint()
alpha1 <- Arg(A[1L] + 1i*A[2L])
alpha2 <- Arg(B[1L] + 1i*B[2L])
path2 <- Arc$new(c(0,0), 3, alpha1, alpha2, FALSE)$path()
polypath(rbind(path1,path2), col = colors[2L])
}))

Thetas <- lapply(
rapply(Ms, function(M){
Arg(M[1L] + 1i*M[2L])
}, how="replace"),
unlist)

for(d in 2L:depth){
thetas <- sort(unlist(Thetas[1L:d]))
for(i in 1L:length(thetas)){
ip1 <- ifelse(i == length(thetas), 1L, i+1L)
arc <- circ$orthogonalThroughTwoPointsOnCircle(thetas[i], thetas[ip1], arc = TRUE) path1 <- arc$path()
B <- arc$startingPoint() A <- arc$endingPoint()
alpha1 <- Arg(A[1L] + 1i*A[2L])
alpha2 <- Arg(B[1L] + 1i*B[2L])
path2 <- Arc$new(c(0,0), 3, alpha1, alpha2, FALSE)$path()
polypath(rbind(path1,path2), col = colors[d+1L])
}
}

draw(circ, border = "white")

par(opar)

invisible()
}

tessellation2(
depth = 5L,
Thetas0 = c(0, 2, 3.8),
colors = viridisLite::viridis(6)
)


# Director circle of an ellipse

Let's draw the director circle of an ellipse. We start by constructing the minimum bounding box of the ellipse.

ell <- Ellipse$new(c(1,1), 5, 2, 30) majorAxis <- ell$diameter(0)
minorAxis <- ell$diameter(pi/2) v1 <- (majorAxis$B - majorAxis$A) / 2 v2 <- (minorAxis$B - minorAxis$A) / 2 # sides of the minimum bounding box side1 <- majorAxis$translate(v2)
side2 <- majorAxis$translate(-v2) side3 <- minorAxis$translate(v1)
side4 <- minorAxis$translate(-v1) # take a vertex of the bounding box A <- side1$A
# director circle
circ <- CircleOA(ell$center, A)  Now let's take a tangent$T_1$to the ellipse, construct the half-line directed by$T_1$with origin the point of tangency, determine the intersection point of this half-line with the director circle, and draw the perpendicular$T_2$of$T_1$passing by this intersection point. Then$T_2$is another tangent to the ellipse. T1 <- ell$tangent(0.3)
halfT1 <- T1$clone(deep = TRUE) halfT1$extendA <- FALSE
I <- intersectionCircleLine(circ, halfT1, strict = TRUE)
T2 <- T1$perpendicular(I)  opar <- par(mar=c(0,0,0,0)) plot(NULL, asp = 1, xlim = c(-3,6), ylim = c(-5,7), xlab = NA, ylab = NA) # draw the ellipse draw(ell, col = "blue") # draw the bounding box draw(side1, lwd = 2, col = "green") draw(side2, lwd = 2, col = "green") draw(side3, lwd = 2, col = "green") draw(side4, lwd = 2, col = "green") # draw the director circle draw(circ, lwd = 2, border = "red") # draw the two tangents draw(T1); draw(T2) # restore the graphical parameters par(opar)  # Playing with Steiner chains The PlaneGeometry package has a function SteinerChain which generates a Steiner chain of circles. ## Elliptical Steiner chain By applying an affine transformation to a Steiner chain, we can get an elliptical Steiner chain. c0 <- Circle$new(c(3,0), 3) # exterior circle
circles <- SteinerChain(c0, 3, -0.2, 0.5)
# take an ellipse
ell <- Ellipse$new(c(-4,0), 4, 2.5, 140) # take the affine transformation which maps the exterior circle to this ellipse f <- AffineMappingEllipse2Ellipse(c0, ell) # take the images of the Steiner circles by this transformation ellipses <- lapply(circles, f$transformEllipse)

opar <- par(mar = c(0,0,0,0))
plot(NULL, asp = 1, xlim = c(-8,6), ylim = c(-4,4),
xlab = NA, ylab = NA, axes = FALSE)
# draw the Steiner chain
invisible(lapply(circles, draw, lwd = 2, col = "blue"))
draw(c0, lwd = 2)
# draw the elliptical Steiner chain
invisible(lapply(ellipses, draw, lwd = 2, col = "red", border = "forestgreen"))
draw(ell, lwd = 2, border = "forestgreen")
par(opar)


Here is how I got the animation below, by varying the shift parameter of the Steiner chain.

library(gifski)

c0 <- Circle$new(c(3,0), 3) ell <- Ellipse$new(c(-4,0), 4, 2.5, 140)
f <- AffineMappingEllipse2Ellipse(c0, ell)

fplot <- function(shift){
circles <- SteinerChain(c0, 3, -0.2, shift)
ellipses <- lapply(circles, f$transformEllipse) opar <- par(mar = c(0,0,0,0)) plot(NULL, asp = 1, xlim = c(-8,0), ylim = c(-4,4), xlab = NA, ylab = NA, axes = FALSE) invisible(lapply(ellipses, draw, lwd = 2, col = "blue", border = "black")) draw(ell, lwd = 2) par(opar) invisible() } shift_ <- seq(0, 3, length.out = 100)[-1L] save_gif( for(shift in shift_){ fplot(shift) }, "SteinerChainElliptical.gif", 512, 512, 1/12, res = 144 ) {width=60%} ## Nested Steiner chains We can choose the exterior circle of the Steiner chain. Therefore, given a circle of a Steiner chain, we can nest another Steiner chain in this circle. c0 <- Circle$new(c(3,0), 3) # exterior circle
circles <- SteinerChain(c0, 3, -0.2, 0.5)
# Steiner chain for each circle, except the small one (it is too small)
chains <- lapply(circles[1:3], function(c0){
SteinerChain(c0, 3, -0.2, 0.5)
})

opar <- par(mar = c(0,0,0,0))
plot(NULL, asp = 1, xlim = c(0,6), ylim = c(-4,4),
xlab = NA, ylab = NA, axes = FALSE)
# draw the big Steiner chain
invisible(lapply(circles, draw, lwd = 2, border = "blue"))
draw(c0, lwd = 2)
# draw the nested Steiner chain
invisible(lapply(chains, function(circles){
lapply(circles, draw, lwd = 2, border = "red")
}))
par(opar)


Of course you can also nest elliptical Steiner chains, and animate the picture!

# Elliptical billiard

The following code plots the trajectory of a ball on an elliptical billiard.

reflect <- function(incidentDir, normalVec){
incidentDir - 2*c(crossprod(normalVec, incidentDir)) * normalVec
}

# n: number of segments; P0: initial point; v0: initial direction
trajectory <- function(n, P0, v0){
out <- vector("list", n)
L <- Line$new(P0, P0+v0) inters <- intersectionEllipseLine(ell, L) Q0 <- inters$I2
out[] <- Line$new(inters$I1, inters$I2, FALSE, FALSE) for(i in 2:n){ theta <- atan2(Q0, Q0) t <- ell$theta2t(theta, degrees = FALSE)
nrmlVec <- ell$normal(t) v <- reflect(Q0-P0, nrmlVec) inters <- intersectionEllipseLine(ell, Line$new(Q0, Q0+v))
out[[i]] <- Line$new(inters$I1, inters$I2, FALSE, FALSE) P0 <- Q0 Q0 <- if(isTRUE(all.equal(Q0, inters$I1))) inters$I2 else inters$I1
}
out
}

ell <- Ellipse$new(c(0,0), 6, 3, 0) P0 <- ell$pointFromAngle(60)
v0 <- c(cos(pi+0.8), sin(pi+0.8))
traj <- trajectory(150, P0, v0)

opar <- par(mar = c(0,0,0,0))
plot(NULL, asp = 1, xlim = c(-7,7), ylim = c(-4,4),
xlab = NA, ylab = NA, axes = FALSE)
draw(ell, border = "red", col = "springgreen", lwd = 3)
invisible(lapply(traj, draw))
par(opar)


Run the code below to see an animated trajectory:

opar <- par(mar = c(0,0,0,0))
plot(NULL, asp = 1, xlim = c(-7,7), ylim = c(-4,4),
xlab = NA, ylab = NA, axes = FALSE)
draw(ell, border = "red", col = "springgreen", lwd = 3)
for(i in 1:length(traj)){
draw(traj[[i]])
Sys.sleep(0.3)
}
par(opar)


# An illustration of inversions

A generalized circle is either a circle or a line. The following code generates a family of generalized circles by repeatedly applying inversions:

# generation 0
angles <- c(0, pi/2, pi, 3*pi/2)
bigCircle <- Circle$new(center = c(0, 0), radius = 2) attr(bigCircle, "gen") <- 0L attr(bigCircle, "base") <- length(angles) + 1L gen0 <- c( lapply(seq_along(angles), function(i){ beta <- angles[i] circle <- Circle$new(center = c(cos(beta), sin(beta)), radius = 1)
attr(circle, "gen") <- 0L
attr(circle, "base") <- i
circle
}),
list(
bigCircle
)
)
n0 <- length(gen0)

# generations 1, 2, 3
generations <- vector("list", length = 4L)
generations[[1L]] <- gen0
for(g in 2L:4L){
gen <- generations[[g-1L]]
n <- length(gen)
n1 <- n*(n0 - 1L)
gen_new <- vector("list", length = n1)
k <- 0L
while(k < n1){
for(j in 1L:n){
gcircle_j <- gen[[j]]
base <- attr(gcircle_j, "base")
for(i in 1L:n0){
if(i != base){
k <- k + 1L
circ <- gen0[[i]]
iota <- Inversion$new(pole = circ$center, power = circ$radius^2) gcircle <- iota$invertGcircle(gcircle_j)
attr(gcircle, "gen") <- g - 1L
attr(gcircle, "base") <- i
gen_new[[k]] <- gcircle
}
}
}
}
generations[[g]] <- gen_new
}

gcircles <- c(
generations[[1L]], generations[[2L]], generations[[3L]], generations[[4L]]
)


There are 425 generalized circles:

length(gcircles)


But some of them are duplicated. In order to remove the duplicates, I will use the following function uniqueWith, which takes as arguments a list or a vector and a function representing a binary relation between the elements of this collection:

uniqueWith <- function(v, f){
size <- length(v)
for(i in seq_len(size-1L)){
j <- i + 1L
while(j <= size){
if(f(v[[i]], v[[j]])){
v <- v[-j]
size <- size - 1L
}else{
j <- j + 1L
}
}
}
v[1L:size]
}


For example:

uniqueWith(
c(a = "you", b = "are", c = "great"),
function(x, y) nchar(x) == nchar(y)
)


So we can remove the duplicated generalized circles as follows:

gcircles <- uniqueWith(
gcircles,
function(g1, g2){
class(g1)[1L] == class(g2)[1L] && g1$isEqual(g2) } )  Now it remains 161 generalized circles: length(gcircles)  Let's write a helper function to draw these generalized circles: drawGcircle <- function(gcircle, colors = rainbow(4L), ...){ gen <- attr(gcircle, "gen") if(is(gcircle, "Circle")){ draw(gcircle, border = colors[1L + gen], ...) }else{ draw(gcircle, col = colors[1L + gen], ...) } }  And now let's draw them: opar <- par(mar = c(0,0,0,0), bg = "black") plot(0, 0, type = "n", xlim = c(-2.3, 2.3), ylim = c(-2.3, 2.3), asp = 1, axes = FALSE, xlab = NA, ylab = NA) invisible(lapply(gcircles, drawGcircle, lwd = 2)) par(opar)  # Schottky circles This construction is taken from the book Indra's Pearls: The Vision of Felix Klein. library(freegroup) a <- alpha(1) A <- inverse(a) b <- alpha(2) B <- inverse(b) # words of size n n <- 6L G <- do.call(expand.grid, rep(list(c("a", "A", "b", "B")), n)) G <- split(as.matrix(G), 1:nrow(G)) G <- lapply(G, function(w){ sum(do.call(c.free, lapply(w, function(x) switch(x, a=a, A=A, b=b, B=B)))) }) G <- uniqueWith(G, free_equal) sizes <- vapply(G, total, numeric(1L)) Gn <- G[sizes == n] # starting circles #### Ca <- Line$new(c(-1,0), c(1,0))
Rc <- sqrt(2)/4
yI <- -3*sqrt(2)/4
CA <- Circle$new(c(0,yI), Rc) theta <- -0.5 T <- c(Rc*cos(theta), yI+Rc*sin(theta)) P <- c(T+T*tan(theta), 0) PT <- sqrt(c(crossprod(T-P))) xTprime <- P+PT xPprime <- -yI/tan(theta) PprimeTprime <- abs(xTprime-xPprime) Rcprime <- abs(yI*PprimeTprime/xPprime) Cb <- Circle$new(c(xTprime, -Rcprime), Rcprime)
CB <- Circle$new(c(-xTprime, -Rcprime), Rcprime) GCIRCLES <- list(a = Ca, A = CA, b = Cb, B = CB) # Mobius transformations #### Mob_a <- Mobius$new(rbind(c(sqrt(2), 1i), c(-1i, sqrt(2))))
Mob_A <- Mob_a$inverse() toCplx <- function(xy) complex(real = xy, imaginary = xy) Mob_b <- Mobius$new(rbind(
c(toCplx(Cb$center), c(crossprod(Cb$center))-Cb$radius^2), c(1, -toCplx(CB$center))
))
Mob_B <- Mob_b$inverse() MOBS <- list(a = Mob_a, A = Mob_A, b = Mob_b, B = Mob_B) # Conversion word of size n to circle word2seq <- function(g){ seq <- c() gr <- reduce(g)[[1L]] for(j in 1L:ncol(gr)){ monomial <- gr[, j] t <- c("a", "b")[monomial[1L]] i <- monomial[2L] if(i < 0L){ i <- -i t <- toupper(t) } seq <- c(seq, rep(t, i)) } seq } word2circle <- function(g){ seq <- word2seq(g) mobs <- MOBS[seq] mobius <- Reduce(function(M1, M2) M1$compose(M2), mobs[-n])
mobius$transformGcircle(GCIRCLES[[seq[n]]]) }  Here is the picture: opar <- par(mar = c(0,0,0,0), bg = "black") plot(NULL, asp = 1, xlim = c(-3,3), ylim = c(-3,3), axes = FALSE, xlab = NA, ylab = NA) draw(Ca); draw(CA); draw(Cb); draw(CB) C1 <- Mob_A$transformCircle(CA)
C2 <- Mob_A$transformCircle(CB) C3 <- Mob_A$transformCircle(Cb)
draw(C1, lwd = 2, border = "red")
draw(C2, lwd = 2, border = "red")
draw(C3, lwd = 2, border = "red")
C1 <- Mob_a$transformLine(Ca) C2 <- Mob_a$transformCircle(Cb)
C3 <- Mob_a$transformCircle(CB) draw(C1, lwd = 2, border = "green") draw(C2, lwd = 2, border = "green") draw(C3, lwd = 2, border = "green") C1 <- Mob_b$transformLine(Ca)
C2 <- Mob_b$transformCircle(CA) C3 <- Mob_b$transformCircle(Cb)
draw(C1, lwd = 2, border = "blue")
draw(C2, lwd = 2, border = "blue")
draw(C3, lwd = 2, border = "blue")
C1 <- Mob_B$transformLine(Ca) C2 <- Mob_B$transformCircle(CA)
C3 <- Mob_B$transformCircle(CB) draw(C1, lwd = 2, border = "yellow") draw(C2, lwd = 2, border = "yellow") draw(C3, lwd = 2, border = "yellow") for(g in Gn){ circ <- word2circle(g) draw(circ, lwd = 2, border = "orange") } par(opar)  Here is the same picture but with better quality: {width=60%} I realized this picture with the tikzDevice package. # Modular tesselation I did this animation after I came across the paper Complex Variables Visualized written by Thomas Ponweiser. This is the paper which motivated me to implement the generalized power of a Möbius transformation. library(elliptic) # for the unimodular matrices # Möbius transformations T <- Mobius$new(rbind(c(0,-1), c(1,0)))
U <- Mobius$new(rbind(c(1,1), c(0,1))) R <- U$compose(T)
# R**t, generalized power
Rt <- function(t){
R$gpower(t) } # starting circles I <- Circle$new(c(0, 1.5), 0.5)
TI <- T$transformCircle(I) # modified Cayley transformation Phi <- Mobius$new(rbind(c(1i, 1), c(1, 1i)))

draw_pair <- function(M, u, compose = FALSE){
if(compose) M <- M$compose(T) A <- M$compose(Rt(u))$compose(Phi) C <- A$transformCircle(I)
draw(C, col = "magenta")
C <- A$transformCircle(TI) draw(C, col = "magenta") if(!compose){ draw_pair(M, u, compose=TRUE) } } n <- 8L transfos <- unimodular(n) fplot <- function(u){ opar <- par(mar = c(0,0,0,0), bg = "black") plot(NULL, asp = 1, xlim = c(-1.1, 1.1), ylim = c(-1.1, 1.1), axes = FALSE, xlab = NA, ylab = NA) for(i in 1L:dim(transfos)[3L]){ transfo <- transfos[, , i] M <- Mobius$new(transfo)
draw_pair(M, u)
M <- M$inverse() draw_pair(M, u) diag(transfo) <- -diag(transfo) M <- Mobius$new(transfo)
draw_pair(M, u)
M <- M$inverse() draw_pair(M, u) d <- diag(transfo) if(d[1L] != d[2L]){ diag(transfo) <- rev(diag(transfo)) M <- Mobius$new(transfo)
draw_pair(M, u)
M <- M$inverse() draw_pair(M, u) } } }  To get the animation, run: library(gifski) u_ <- seq(0, 3, length.out = 181)[-1] save_gif( for(u in u_){ fplot(u) }, width = 512, height = 512, delay = 0.1 ) {width=60%} # Apollonian gasket It is not hard to draw an Apollonian gasket with PlaneGeometry. We do a function, in order to use it later to do an animation. # function to construct the "children" #### ApollonianChildren <- function(inversions, circles1){ m <- length(inversions) n <- length(circles1) circles2 <- list() for(i in 1:n){ circ <- circles1[[i]] k <- attr(circ, "inversion") for(j in 1:m){ if(j != k){ circle <- inversions[[j]]$invertCircle(circ)
attr(circle, "inversion") <- j
circles2 <- append(circles2, circle)
}
}
}
circles2
}

ApollonianGasket <- function(c0, n, phi, shift, depth){
circles0 <- SteinerChain(c0, n, phi, shift)
# construct the inversions ####
inversions <- vector("list", n + 1)
for(i in 1:n){
inversions[[i]] <- inversionFixingThreeCircles(
c0, circles0[[i]], circles0[[(i %% n) + 1]]
)
}
inversions[[n+1]] <- inversionSwappingTwoCircles(c0, circles0[[n+1]])
# first generation of children
circles1 <- list()
for(i in 1:n){
ip1 <- (i %% n) + 1
for(j in 1:(n+1)){
if(j != i && j != ip1){
circle <- inversions[[i]]$invertCircle(circles0[[j]]) attr(circle, "inversion") <- i circles1 <- append(circles1, circle) } } } # construct children #### allCircles <- vector("list", depth) allCircles[] <- circles0 allCircles[] <- circles1 for(i in 3:depth){ allCircles[[i]] <- ApollonianChildren(inversions, allCircles[[i-1]]) } allCircles }  Let's apply our function: library(viridisLite) # for the colors c0 <- Circle$new(c(0,0), 3) # the exterior circle
depth <- 5
colors <- plasma(depth)
ApollonianCircles <- ApollonianGasket(c0, n = 3, phi = 0.3, shift = 0, depth)
# plot ####
center0 <- c0$center radius0 <- c0$radius
opar <- par(mar = c(0, 0, 0, 0))
plot(NULL, type = "n", xlim = xlim, ylim = ylim,
xlab = NA, ylab = NA, axes = FALSE, asp = 1)
draw(c0, border = "black", lwd = 2)
for(i in 1:depth){
for(circ in ApollonianCircles[[i]]){
draw(circ, col = colors[i])
}
}
par(opar)


We can do an animation now:

fplot <- function(shift){
gasket <- ApollonianGasket(c0, n = 3, phi = 0.3, shift = shift, depth)
par(mar = c(0, 0, 0, 0))
plot(NULL, type = "n", xlim = xlim, ylim = ylim,
xlab = NA, ylab = NA, axes = FALSE, asp = 1)
draw(c0, border = "black", lwd = 2)
for(i in 1:depth){
draw(circ, col = colors[i])
}
}
}

fanim <- function(){
shifts <- seq(0, 3, length.out = 101)[-101]
for(shift in shifts){
fplot(shift)
}
}

library(gifski)
save_gif(
fanim(),
width = 512, height = 512,
delay = 0.1
) {width=50%}

We can also animate the Apollonian gasket with the help of a Möbius transformation. Consider the following complex matrix:

$$M = \begin{pmatrix} i & \gamma \ \bar\gamma & -i \end{pmatrix}$$ with $|\gamma| < 1$.

The Möbius transformation associated to $M$ maps the unit disk to the unit disk and it is of order $2$. Its powers map the unit disk to the unit dis as well. By the way, after some calculus, one can give the expression of $M^t$. We find

Mt <- function(gamma, t){
h <- sqrt(1 - Mod(gamma)^2)
d2 <- h^t * (cos(t*pi/2) + 1i*sin(t*pi/2))
d1 <- Conj(d2)
A11 <- Re(d1) - 1i*Im(d1)/h
A12 <- Im(d2) * gamma / h
rbind(
c(A11, A12),
c(Conj(A11), Conj(A12))
)
}


Now let's do the animation.

c0 <- Circle$new(c(0,0), 1) # the exterior circle depth <- 5 ApollonianCircles <- ApollonianGasket(c0, n = 3, phi = 0.1, shift = 0.5, depth) xlim <- c(-1.1, 1.1) ylim <- c(-1.1, 1.1) opar <- par(mar = c(0, 0, 0, 0)) fplot <- function(gamma, t){ plot(NULL, type = "n", xlim = xlim, ylim = ylim, xlab = NA, ylab = NA, axes = FALSE, asp = 1) draw(c0, border = "black", lwd = 2) Mob <- Mt(gamma, t) for(i in 1:depth){ for(circ in ApollonianCircles[[i]]){ draw(Mob$transformCircle(circ), col = colors[i])
}
}
}
fanim <- function(){
gamma <- 0.5 + 0.4i
t_ <- seq(0, 2, length.out = 91)[-91]
for(t in t_){
fplot(gamma, t)
}
}

library(gifski)
save_gif(
fanim(),
"ApollonianMobius.gif",
width = 512, height = 512,
delay = 0.1
)
par(opar) {width=60%}

# Another Apollonian fractal

Let's do another Apollonian fractal which uses the inner Soddy circle. As you can see, the code is short:

apollony <- function(c1, c2, c3, n){
soddycircle <- soddyCircle(c1, c2, c3)
if(n == 1){
soddycircle
}else{
c(
apollony(c1, c2, soddycircle, n-1),
apollony(c1, soddycircle, c3, n-1),
apollony(soddycircle, c2, c3, n-1)
)
}
}

fractal <- function(n){
c1 = Circle$new(c(1, -1/sqrt(3)), 1) c2 = Circle$new(c(-1, -1/sqrt(3)), 1)
c3 = Circle$new(c(0, sqrt(3) - 1/sqrt(3)), 1) do.call(c, lapply(1:n, function(i) apollony(c1, c2, c3, i))) } circs <- fractal(4)  Let's plot the fractal in 3D with the help of the rgl package. library(rgl) # the spheres in rgl, obtained with the spheres3d function, are not smooth; # the way we use below provides pretty spheres unitSphere <- subdivision3d(icosahedron3d(), depth = 4) unitSphere$vb[4, ] <-
apply(unitSphere$vb[1:3, ], 2, function(x) sqrt(sum(x * x))) unitSphere$normals <- unitSphere$vb drawSphere <- function(circle, ...){ center <- circle$center
radius <- circle$radius sphere <- scale3d(unitSphere, radius, radius, radius) shade3d(translate3d(sphere, center, center, 0), ...) }  Now here is how to plot the fractal and make an animation: # plot #### open3d(windowRect = c(50, 50, 562, 562)) bg3d(color = "#363940") view3d(35, 60, zoom = 0.95) for(circ in circs){ drawSphere(circ, color = "darkred") } # animation #### movie3d( spin3d(axis = c(0, 0, 1), rpm = 15), duration = 4, fps = 15, movie = "Apollony", dir = ".", convert = "magick convert -dispose previous -loop 0 -delay 1x%d %s*.png %s.%s", startTime = 1/60 ) {width=40%} # Malfatti gasket 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 + 1i * M } 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 > B){
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 - A) / distance(A, B))
zX1 <- exp(-1i * alpha) * (zoA - zB)
zX2 <- exp(-1i * alpha) * (zoB - zB)
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, OP)
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 1:3){ 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[]
C2 <- mcircles[]
C3 <- mcircles[]
triangles3d(cbind(rbind(A, B, C), 0), col = "yellow", alpha = 0.2)
holes <- list(
side.circle.circle(A, B, C1, C2, color = colors),
side.circle.circle(B, C, C2, C3, color = colors),
side.circle.circle(C, A, C3, C1, color = colors),
innerSoddyCircle(C1, C2, C3, color = colors),
side.side.circle(A, B, C, C1, color = colors),
side.side.circle(B, A, C, C2, color = colors),
side.side.circle(C, A, B, C3, color = colors)
)
for(d in 1:depth){
n_holes <- length(holes)
Holes <- list()
for(i in 1:n_holes){
Holes <- append(Holes, Newholes(holes[[i]], colors[d + 1]))
}
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 <- 3
colors <- viridis(depth + 1)
n1 <- 3
n2 <- 4
n3 <- 5
depth2 <- 3
phi1 <- 0.2
phi2 <- 0.3
phi3 <- 0.4
shift <- 0
colors2 <- plasma(depth2)


And we get the 3D picture:

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[], n1, phi1, shift, depth2, colors2)
drawCircularGasket(mcircles[], n2, phi2, shift, depth2, colors2)
drawCircularGasket(mcircles[], n3, phi3, shift, depth2, colors2) {width=50%}

As an exercise, you can add some animation to this picture, by animating the three circular gaskets.

## Try the PlaneGeometry package in your browser

Any scripts or data that you put into this service are public.

PlaneGeometry documentation built on May 31, 2023, 7:37 p.m.