Examples of the PlaneGeometry package

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) )
#> [1] 2.621359 1.158114
# check whether the Exeter point is also on (FC')
lineFCprime$includes(ExeterPoint)
#> [1] TRUE

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
(inradius^2 + semiperimeter^2) / (4*inradius)
#> [1] 11.15942
ApolloniusCircle$radius
#> [1] 11.15942

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
)

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)
#> [1] 425

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)
)
#>       a       c 
#>   "you" "great"

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)
#> [1] 158

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:

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
)

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
xlim <- center0[1] + c(-radius0 - 0.1, radius0 + 0.1)
ylim <- center0[2] + c(-radius0 - 0.1, radius0 + 0.1)
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){
    for(circ in gasket[[i]]){
      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(),
  "ApollonianGasket.gif",
  width = 512, height = 512,
  delay = 0.1
)

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)

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
)

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)
  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[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)

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