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)[[2]] Bprime <- intersectionCircleLine(circumcircle, medianB)[[2]] Cprime <- intersectionCircleLine(circumcircle, medianC)[[1]] 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[[1]]; B <- intersections[[2]] points(rbind(A,B), pch = 19, col = c("red", "blue")) theta1 <- Arg((A-O1)[1] + 1i*(A-O1)[2]) theta2 <- Arg((B-O1)[1] + 1i*(B-O1)[2]) path1 <- Arc$new(O1, circ1$radius, theta1, theta2, FALSE)$path()

theta1 <- Arg((A-O2)[1] + 1i*(A-O2)[2])
theta2 <- Arg((B-O2)[1] + 1i*(B-O2)[2])
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[[1]] <- Line$new(inters$I1, inters$I2, FALSE, FALSE) for(i in 2:n){ theta <- atan2(Q0[2], Q0[1]) 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[1]+T[2]*tan(theta), 0) PT <- sqrt(c(crossprod(T-P))) xTprime <- P[1]+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[1], imaginary = xy[2]) 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[[1]] <- circles0 allCircles[[2]] <- 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[1], center[2], 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[1] + 1i * M[2] } 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[2] > B[2]){
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[1] - A[1]) / 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[2], OP[1])
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[[1]]
C2 <- mcircles[[2]]
C3 <- mcircles[[3]]
triangles3d(cbind(rbind(A, B, C), 0), col = "yellow", alpha = 0.2)
holes <- list(
side.circle.circle(A, B, C1, C2, color = colors[1]),
side.circle.circle(B, C, C2, C3, color = colors[1]),
side.circle.circle(C, A, C3, C1, color = colors[1]),
innerSoddyCircle(C1, C2, C3, color = colors[1]),
side.side.circle(A, B, C, C1, color = colors[1]),
side.side.circle(B, A, C, C2, color = colors[1]),
side.side.circle(C, A, B, C3, color = colors[1])
)
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[[1]], n1, phi1, shift, depth2, colors2)
drawCircularGasket(mcircles[[2]], n2, phi2, shift, depth2, colors2)
drawCircularGasket(mcircles[[3]], 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.