Nothing
Examples of the PlaneGeometry package
In PlaneGeometry: Plane Geometry
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
(inradius^2 + semiperimeter^2) / (4*inradius)
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
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
)
{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)
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)
{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.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
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 (inradius^2 + semiperimeter^2) / (4*inradius) 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 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 )
{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) 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)
{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.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.