Nothing
inst/doc/examples.R
In PlaneGeometry: Plane Geometry
## ---- include = FALSE---------------------------------------------------------
knitr::opts_chunk$set(
collapse = TRUE,
comment = "#>"
)
## ----setup--------------------------------------------------------------------
library(PlaneGeometry)
## -----------------------------------------------------------------------------
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)
## ----Exeter, fig.width=4, fig.height=4----------------------------------------
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)
## -----------------------------------------------------------------------------
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()
## ----tangentCircles, fig.width=4, fig.height=4--------------------------------
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)
## -----------------------------------------------------------------------------
# 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)
## ----apollonius, fig.width=4, fig.height=4------------------------------------
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)
## -----------------------------------------------------------------------------
inradius <- t$inradius()
semiperimeter <- sum(t$edges()) / 2
(inradius^2 + semiperimeter^2) / (4*inradius)
ApolloniusCircle$radius
## ----lappingArea, fig.width=4, fig.height=4-----------------------------------
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)
## -----------------------------------------------------------------------------
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, fig.width=4, fig.height=4, eval=require("viridisLite", quietly=TRUE)----
tessellation(
depth = 5L,
Thetas0 = c(0, 2, 3.8),
colors = viridisLite::viridis(5L)
)
## -----------------------------------------------------------------------------
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, fig.width=4, fig.height=4, eval=require("viridisLite", quietly=TRUE)----
tessellation2(
depth = 5L,
Thetas0 = c(0, 2, 3.8),
colors = viridisLite::viridis(6L)
)
## -----------------------------------------------------------------------------
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)
## ---- message=FALSE-----------------------------------------------------------
T1 <- ell$tangent(0.3)
halfT1 <- T1$clone(deep = TRUE)
halfT1$extendA <- FALSE
I <- intersectionCircleLine(circ, halfT1, strict = TRUE)
T2 <- T1$perpendicular(I)
## ----directorCircle, fig.width=4, fig.height=4--------------------------------
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)
## -----------------------------------------------------------------------------
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)
## ----ellipticalSteiner, fig.width=4, fig.height=4-----------------------------
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)
## ---- eval=FALSE--------------------------------------------------------------
# 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
# )
## -----------------------------------------------------------------------------
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)
})
## ----nestedSteiner, fig.width=4, fig.height=4---------------------------------
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)
## ----billiard, fig.width=4, fig.height=4--------------------------------------
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)
## ---- eval=FALSE--------------------------------------------------------------
# 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)
## ----gcircles-----------------------------------------------------------------
# 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]]
)
## -----------------------------------------------------------------------------
length(gcircles)
## ----uniqueWith---------------------------------------------------------------
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]
}
## -----------------------------------------------------------------------------
uniqueWith(
c(a = "you", b = "are", c = "great"),
function(x, y) nchar(x) == nchar(y)
)
## -----------------------------------------------------------------------------
gcircles <- uniqueWith(
gcircles,
function(g1, g2){
class(g1)[1L] == class(g2)[1L] && g1$isEqual(g2)
}
)
## -----------------------------------------------------------------------------
length(gcircles)
## ----drawGcircle--------------------------------------------------------------
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], ...)
}
}
## ----draw_gcircles, fig.width=4, fig.height=4---------------------------------
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, message=FALSE, eval=require("freegroup", quietly=TRUE)---------
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]]])
}
## ----draw_schottky, fig.width=4, fig.height=4, eval=require("freegroup", quietly=TRUE)----
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)
## ----modularTesselation, message=FALSE, eval=require("elliptic", quietly=TRUE)----
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)
}
}
}
## ---- eval=FALSE--------------------------------------------------------------
# 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
# )
## ----gasket-------------------------------------------------------------------
# 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
}
## ----drawgasket, fig.height=4, fig.width=4, eval=require("viridisLite", quietly=TRUE)----
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)
## ----animgasket, eval=FALSE---------------------------------------------------
# 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
# )
## -----------------------------------------------------------------------------
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))
)
}
## ----ApollonianMobius---------------------------------------------------------
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)
}
}
## ---- eval=FALSE--------------------------------------------------------------
# library(gifski)
# save_gif(
# fanim(),
# "ApollonianMobius.gif",
# width = 512, height = 512,
# delay = 0.1
# )
# par(opar)
## ----apollony, eval=FALSE-----------------------------------------------------
# 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)
## ----rglsphere, eval=require("rgl", quietly=TRUE)-----------------------------
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 = 4L)
unitSphere$vb[4L, ] <-
apply(unitSphere$vb[1L:3L, ], 2L, function(x) sqrt(sum(x * x)))
unitSphere$normals <- unitSphere$vb
drawSphere <- function(circle, ...) {
center <- circle$center
radius <- abs(circle$radius)
sphere <- scale3d(unitSphere, radius, radius, radius)
shade3d(translate3d(sphere, center[1L], center[2L], 0), ...)
}
## ----rglapollony, eval=FALSE--------------------------------------------------
# # 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
# )
## ----malfattigasket, eval=require("viridisLite", quietly=TRUE)----------------
toCplx <- function(M) {
M[1L] + 1i * M[2L]
}
fromCplx <- function(z) {
c(Re(z), Im(z))
}
distance <- function(A, B) {
sqrt(c(crossprod(B - A)))
}
innerSoddyRadius <- function(r1, r2, r3) {
1 / (1/r1 + 1/r2 + 1/r3 + 2 * sqrt(1/r1/r2 + 1/r2/r3 + 1/r3/r1))
}
innerSoddyCircle <- function(c1, c2, c3, ...) {
radius <- innerSoddyRadius(c1$radius, c3$radius, c3$radius)
center <- Triangle$new(c1$center, c2$center, c3$center)$equalDetourPoint()
c123 <- Circle$new(center, radius)
drawSphere(c123, ...)
list(
list(type = "ccc", c1 = c123, c2 = c1, c3 = c2),
list(type = "ccc", c1 = c123, c2 = c2, c3 = c3),
list(type = "ccc", c1 = c123, c2 = c1, c3 = c3)
)
}
side.circle.circle <- function(A, B, cA, cB, ...) {
if(A[2L] > B[2L]){
return(side.circle.circle(B, A, cB, cA, ...))
}
rA <- cA$radius
rB <- cB$radius
zoA <- toCplx(cA$center)
zoB <- toCplx(cB$center)
zB <- toCplx(A)
alpha <- acos((B[1L] - A[1L]) / distance(A, B))
zX1 <- exp(-1i * alpha) * (zoA - zB)
zX2 <- exp(-1i * alpha) * (zoB - zB)
soddyR <- innerSoddyRadius(rA, rB, Inf)
if(Re(zX1) < Re(zX2)) {
Y <- (2 * rA * sqrt(rB) / (sqrt(rA) + sqrt(rB)) + Re(zX1)) +
sign(Im(zX1)) * 1i * soddyR
} else {
Y <- (2 * rB * sqrt(rA) / (sqrt(rA) + sqrt(rB)) + Re(zX2)) +
sign(Im(zX1)) * 1i * soddyR
}
M <- fromCplx(Y * exp(1i * alpha) + zB)
cAB <- Circle$new(M, soddyR)
drawSphere(cAB, ...)
list(
list(type = "ccc", c1 = cAB, c2 = cA, c3 = cB),
list(type = "ccl", cA = cA, cB = cAB, A = A, B = B),
list(type = "ccl", cA = cAB, cB = cB, A = A, B = B)
)
}
side.side.circle <- function(A, B, C, circle, ...) {
zA <- toCplx(A)
zO <- toCplx(circle$center)
vec <- zA - zO
P <- fromCplx(zO + circle$radius * vec / Mod(vec))
OP <- P - circle$center
onTangent <- P + c(-OP[2L], OP[1L])
L1 <- Line$new(P, onTangent)
P1 <- intersectionLineLine(L1, Line$new(A, C))
P2 <- intersectionLineLine(L1, Line$new(A, B))
incircle <- Triangle$new(A, P1, P2)$incircle()
drawSphere(incircle, ...)
list(
list(type = "cll", A = A, B = B, C = C, circle = incircle),
list(type = "ccl", cA = circle, cB = incircle, A = A, B = B),
list(type = "ccl", cA = circle, cB = incircle, A = A, B = C)
)
}
Newholes <- function(holes, color) {
newholes <- list()
for(i in 1L:3L) {
hole <- holes[[i]]
holeType <- hole[["type"]]
if(holeType == "ccc") {
x <- with(hole, innerSoddyCircle(c1, c2, c3, color = color))
} else if(holeType == "ccl") {
x <- with(hole, side.circle.circle(A, B, cA, cB, color = color))
} else if (holeType == "cll") {
x <- with(hole, side.side.circle(A, B, C, circle, color = color))
}
newholes <- c(newholes, list(x))
}
newholes
}
MalfattiCircles <- function(A, B, C) {
Triangle$new(A, B, C)$MalfattiCircles()
}
drawTriangularGasket <- function(mcircles, A, B, C, colors, depth) {
C1 <- mcircles[[1L]]
C2 <- mcircles[[2L]]
C3 <- mcircles[[3L]]
triangles3d(cbind(rbind(A, B, C), 0), col = "yellow", alpha = 0.2)
holes <- list(
side.circle.circle(A, B, C1, C2, color = colors[1L]),
side.circle.circle(B, C, C2, C3, color = colors[1L]),
side.circle.circle(C, A, C3, C1, color = colors[1L]),
innerSoddyCircle(C1, C2, C3, color = colors[1L]),
side.side.circle(A, B, C, C1, color = colors[1L]),
side.side.circle(B, A, C, C2, color = colors[1L]),
side.side.circle(C, A, B, C3, color = colors[1L])
)
for(d in 1L:depth) {
n_holes <- length(holes)
Holes <- list()
for(i in 1L:n_holes) {
Holes <- append(Holes, Newholes(holes[[i]], colors[d + 1L]))
}
holes <- do.call(list, Holes)
}
}
drawCircularGasket <- function(c0, n, phi, shift, depth, colors) {
ApollonianCircles <- ApollonianGasket(c0, n, phi, shift, depth)
for(i in 1:depth) {
for(circ in ApollonianCircles[[i]]){
drawSphere(circ, color = colors[i])
}
}
}
library(viridisLite)
A <- c(-5, -4)
B <- c(5, -2)
C <- c(0, 6)
mcircles <- MalfattiCircles(A, B, C)
depth <- 3L
colors <- viridis(depth + 1L)
n1 <- 3L
n2 <- 4L
n3 <- 5L
depth2 <- 3L
phi1 <- 0.2
phi2 <- 0.3
phi3 <- 0.4
shift <- 0
colors2 <- plasma(depth2)
## ----rgl, eval=FALSE----------------------------------------------------------
# library(rgl)
# open3d(windowRect = c(50, 50, 562, 562), zoom = 0.9)
# bg3d(rgb(54, 57, 64, maxColorValue = 255))
# drawTriangularGasket(mcircles, A, B, C, colors, depth)
# drawCircularGasket(mcircles[[1L]], n1, phi1, shift, depth2, colors2)
# drawCircularGasket(mcircles[[2L]], n2, phi2, shift, depth2, colors2)
# drawCircularGasket(mcircles[[3L]], n3, phi3, shift, depth2, colors2)
## ----malfattiCircular, warning=FALSE, eval=require("rgl", quietly=TRUE)-------
library(rgl)
iteration <- function(circlesWithIndicator, inversions) {
out <- list()
for(j in seq_along(circlesWithIndicator)) {
circle <- circlesWithIndicator[[j]][["circle"]]
indic <- circlesWithIndicator[[j]][["indic"]]
for(i in 1L:4L) {
if(i != indic) {
circleWithIndicator <- list(
"circle" = inversions[[i]]$invertCircle(circle),
"indic" = i
)
out <- append(out, list(circleWithIndicator))
}
}
}
out
}
gasket <- function(circlesWithIndicator, inversions, depth, colors) {
if(depth > 0){
circlesWithIndicator <- iteration(circlesWithIndicator, inversions)
for(i in seq_along(circlesWithIndicator)) {
drawSphere(circlesWithIndicator[[i]]$circle, color = colors[1L])
}
colors <- colors[-1L]
gasket(circlesWithIndicator, inversions, depth-1L, colors)
}
}
drawGasket <- function(triangle, depth, colors) {
Mcircles <- triangle$MalfattiCircles()
Mtriangle <- Triangle$new(
Mcircles[[1L]]$center, Mcircles[[2L]]$center, Mcircles[[3L]]$center
)
soddyO <- Mtriangle$outerSoddyCircle()
Mcircles <- append(Mcircles, list(soddyO))
for(i in 1L:4L) {
lines3d(
cbind(Mcircles[[i]]$asEllipse()$path(), 0),
color = "black", lwd = 2
)
}
inversions <- vector("list", 4L)
circlesWithIndicator <- vector("list", 4L)
inversions[[1L]] <- inversionFixingThreeCircles(
soddyO, Mcircles[[2L]], Mcircles[[3L]]
)
inversions[[2L]] <- inversionFixingThreeCircles(
soddyO, Mcircles[[1L]], Mcircles[[3L]]
)
inversions[[3L]] <- inversionFixingThreeCircles(
soddyO, Mcircles[[1L]], Mcircles[[2L]]
)
inversions[[4L]] <- inversionFixingThreeCircles(
Mcircles[[1L]], Mcircles[[2L]], Mcircles[[3L]]
)
for(i in 1L:4L) {
circlesWithIndicator[[i]] <-
list("circle" = inversions[[i]]$invertCircle(Mcircles[[i]]), "indic" = i)
drawSphere(circlesWithIndicator[[i]]$circle, color = colors[1L])
}
colors <- colors[-1L]
gasket(circlesWithIndicator, inversions, depth, colors)
}
CircularMalfattiGasket <- function(C, depth, colors) {
A <- c(0,0); B <- c(1,0)
t <- Triangle$new(A, B, C)
Mcircles <- t$MalfattiCircles()
Mtriangle <- Triangle$new(
Mcircles[[1L]]$center, Mcircles[[2L]]$center, Mcircles[[3L]]$center
)
soddyO <- Mtriangle$outerSoddyCircle()
center <- soddyO$center; radius = -soddyO$radius
A1 <- (A-center)/radius; B1 <- (B-center)/radius; C1 <- (C-center)/radius;
t1 <- Triangle$new(A1, B1, C1);
drawGasket(t1, depth, colors)
}
## ----plotMalfattiCircularGasket, eval=FALSE-----------------------------------
# open3d(windowRect = 50 + c(0, 0, 900, 300))
# mfrow3d(1, 3)
# view3d(0, 0, zoom = 0.7)
# CircularMalfattiGasket(
# C = c(0, sqrt(3/2)), depth = 2L,
# colors = c("yellow", "orangered", "darkmagenta")
# )
# next3d()
# view3d(0, 0, zoom = 0.7)
# CircularMalfattiGasket(
# C = c(1, sqrt(3/2)), depth = 2L,
# colors = c("yellow", "orangered", "darkmagenta")
# )
# next3d()
# view3d(0, 0, zoom = 0.7)
# CircularMalfattiGasket(
# C = c(2, sqrt(3/2)), depth = 2L,
# colors = c("yellow", "orangered", "darkmagenta")
# )
## ----hyperbola, fig.width=7, fig.height=5-------------------------------------
# take a hyperbola
L1 <- LineFromInterceptAndSlope(0, 2) # asymptote 1
L2 <- LineFromInterceptAndSlope(-2, -0.15) # asymptote 2
M <- c(2, 3) # a point on the hyperbola
hyperbola <- Hyperbola$new(L1, L2, M)
# take a point on the hyperbola and the tangent at this point
OAB <- hyperbola$OAB()
O <- OAB$O; A <- OAB$A; B <- OAB$B
t <- 0.1
P <- O + cosh(t)*A + sinh(t)*B
tgt <- Line$new(P, P + sinh(t)*A + cosh(t)*B)
# the triangle of interest
C <- intersectionLineLine(L1, tgt)
D <- intersectionLineLine(L2, tgt)
trgl <- Triangle$new(O, C, D)
# plot
opar <- par(mar = c(4, 4, 1, 1))
hyperbola$plot(lwd = 2)
draw(L1, col = "red")
draw(L2, col = "red")
text(t(O), "O", pos = 3)
points(t(P), pch = 19, col = "blue")
text(t(P), "P", pos = 4)
draw(tgt, col = "blue", lwd = 2)
text(t(C), "C", pos = 2)
text(t(D), "D", pos = 4)
trgl$plot(add = TRUE, col = "yellow")
par(opar)
# theorem checking: area of the triangle does not depend on
# the choice of P; more precisely, it is equal to ab
trgl$area()
with(hyperbola$abce(), a * b)
Try the PlaneGeometry package in your browser
Any scripts or data that you put into this service are public.
PlaneGeometry documentation built on Aug. 10, 2023, 1:09 a.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.
## ---- include = FALSE---------------------------------------------------------
knitr::opts_chunk$set(
collapse = TRUE,
comment = "#>"
)
## ----setup--------------------------------------------------------------------
library(PlaneGeometry)
## -----------------------------------------------------------------------------
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)
## ----Exeter, fig.width=4, fig.height=4----------------------------------------
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)
## -----------------------------------------------------------------------------
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()
## ----tangentCircles, fig.width=4, fig.height=4--------------------------------
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)
## -----------------------------------------------------------------------------
# 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)
## ----apollonius, fig.width=4, fig.height=4------------------------------------
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)
## -----------------------------------------------------------------------------
inradius <- t$inradius()
semiperimeter <- sum(t$edges()) / 2
(inradius^2 + semiperimeter^2) / (4*inradius)
ApolloniusCircle$radius
## ----lappingArea, fig.width=4, fig.height=4-----------------------------------
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)
## -----------------------------------------------------------------------------
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, fig.width=4, fig.height=4, eval=require("viridisLite", quietly=TRUE)----
tessellation(
depth = 5L,
Thetas0 = c(0, 2, 3.8),
colors = viridisLite::viridis(5L)
)
## -----------------------------------------------------------------------------
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, fig.width=4, fig.height=4, eval=require("viridisLite", quietly=TRUE)----
tessellation2(
depth = 5L,
Thetas0 = c(0, 2, 3.8),
colors = viridisLite::viridis(6L)
)
## -----------------------------------------------------------------------------
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)
## ---- message=FALSE-----------------------------------------------------------
T1 <- ell$tangent(0.3)
halfT1 <- T1$clone(deep = TRUE)
halfT1$extendA <- FALSE
I <- intersectionCircleLine(circ, halfT1, strict = TRUE)
T2 <- T1$perpendicular(I)
## ----directorCircle, fig.width=4, fig.height=4--------------------------------
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)
## -----------------------------------------------------------------------------
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)
## ----ellipticalSteiner, fig.width=4, fig.height=4-----------------------------
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)
## ---- eval=FALSE--------------------------------------------------------------
# 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
# )
## -----------------------------------------------------------------------------
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)
})
## ----nestedSteiner, fig.width=4, fig.height=4---------------------------------
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)
## ----billiard, fig.width=4, fig.height=4--------------------------------------
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)
## ---- eval=FALSE--------------------------------------------------------------
# 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)
## ----gcircles-----------------------------------------------------------------
# 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]]
)
## -----------------------------------------------------------------------------
length(gcircles)
## ----uniqueWith---------------------------------------------------------------
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]
}
## -----------------------------------------------------------------------------
uniqueWith(
c(a = "you", b = "are", c = "great"),
function(x, y) nchar(x) == nchar(y)
)
## -----------------------------------------------------------------------------
gcircles <- uniqueWith(
gcircles,
function(g1, g2){
class(g1)[1L] == class(g2)[1L] && g1$isEqual(g2)
}
)
## -----------------------------------------------------------------------------
length(gcircles)
## ----drawGcircle--------------------------------------------------------------
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], ...)
}
}
## ----draw_gcircles, fig.width=4, fig.height=4---------------------------------
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, message=FALSE, eval=require("freegroup", quietly=TRUE)---------
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]]])
}
## ----draw_schottky, fig.width=4, fig.height=4, eval=require("freegroup", quietly=TRUE)----
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)
## ----modularTesselation, message=FALSE, eval=require("elliptic", quietly=TRUE)----
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)
}
}
}
## ---- eval=FALSE--------------------------------------------------------------
# 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
# )
## ----gasket-------------------------------------------------------------------
# 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
}
## ----drawgasket, fig.height=4, fig.width=4, eval=require("viridisLite", quietly=TRUE)----
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)
## ----animgasket, eval=FALSE---------------------------------------------------
# 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
# )
## -----------------------------------------------------------------------------
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))
)
}
## ----ApollonianMobius---------------------------------------------------------
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)
}
}
## ---- eval=FALSE--------------------------------------------------------------
# library(gifski)
# save_gif(
# fanim(),
# "ApollonianMobius.gif",
# width = 512, height = 512,
# delay = 0.1
# )
# par(opar)
## ----apollony, eval=FALSE-----------------------------------------------------
# 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)
## ----rglsphere, eval=require("rgl", quietly=TRUE)-----------------------------
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 = 4L)
unitSphere$vb[4L, ] <-
apply(unitSphere$vb[1L:3L, ], 2L, function(x) sqrt(sum(x * x)))
unitSphere$normals <- unitSphere$vb
drawSphere <- function(circle, ...) {
center <- circle$center
radius <- abs(circle$radius)
sphere <- scale3d(unitSphere, radius, radius, radius)
shade3d(translate3d(sphere, center[1L], center[2L], 0), ...)
}
## ----rglapollony, eval=FALSE--------------------------------------------------
# # 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
# )
## ----malfattigasket, eval=require("viridisLite", quietly=TRUE)----------------
toCplx <- function(M) {
M[1L] + 1i * M[2L]
}
fromCplx <- function(z) {
c(Re(z), Im(z))
}
distance <- function(A, B) {
sqrt(c(crossprod(B - A)))
}
innerSoddyRadius <- function(r1, r2, r3) {
1 / (1/r1 + 1/r2 + 1/r3 + 2 * sqrt(1/r1/r2 + 1/r2/r3 + 1/r3/r1))
}
innerSoddyCircle <- function(c1, c2, c3, ...) {
radius <- innerSoddyRadius(c1$radius, c3$radius, c3$radius)
center <- Triangle$new(c1$center, c2$center, c3$center)$equalDetourPoint()
c123 <- Circle$new(center, radius)
drawSphere(c123, ...)
list(
list(type = "ccc", c1 = c123, c2 = c1, c3 = c2),
list(type = "ccc", c1 = c123, c2 = c2, c3 = c3),
list(type = "ccc", c1 = c123, c2 = c1, c3 = c3)
)
}
side.circle.circle <- function(A, B, cA, cB, ...) {
if(A[2L] > B[2L]){
return(side.circle.circle(B, A, cB, cA, ...))
}
rA <- cA$radius
rB <- cB$radius
zoA <- toCplx(cA$center)
zoB <- toCplx(cB$center)
zB <- toCplx(A)
alpha <- acos((B[1L] - A[1L]) / distance(A, B))
zX1 <- exp(-1i * alpha) * (zoA - zB)
zX2 <- exp(-1i * alpha) * (zoB - zB)
soddyR <- innerSoddyRadius(rA, rB, Inf)
if(Re(zX1) < Re(zX2)) {
Y <- (2 * rA * sqrt(rB) / (sqrt(rA) + sqrt(rB)) + Re(zX1)) +
sign(Im(zX1)) * 1i * soddyR
} else {
Y <- (2 * rB * sqrt(rA) / (sqrt(rA) + sqrt(rB)) + Re(zX2)) +
sign(Im(zX1)) * 1i * soddyR
}
M <- fromCplx(Y * exp(1i * alpha) + zB)
cAB <- Circle$new(M, soddyR)
drawSphere(cAB, ...)
list(
list(type = "ccc", c1 = cAB, c2 = cA, c3 = cB),
list(type = "ccl", cA = cA, cB = cAB, A = A, B = B),
list(type = "ccl", cA = cAB, cB = cB, A = A, B = B)
)
}
side.side.circle <- function(A, B, C, circle, ...) {
zA <- toCplx(A)
zO <- toCplx(circle$center)
vec <- zA - zO
P <- fromCplx(zO + circle$radius * vec / Mod(vec))
OP <- P - circle$center
onTangent <- P + c(-OP[2L], OP[1L])
L1 <- Line$new(P, onTangent)
P1 <- intersectionLineLine(L1, Line$new(A, C))
P2 <- intersectionLineLine(L1, Line$new(A, B))
incircle <- Triangle$new(A, P1, P2)$incircle()
drawSphere(incircle, ...)
list(
list(type = "cll", A = A, B = B, C = C, circle = incircle),
list(type = "ccl", cA = circle, cB = incircle, A = A, B = B),
list(type = "ccl", cA = circle, cB = incircle, A = A, B = C)
)
}
Newholes <- function(holes, color) {
newholes <- list()
for(i in 1L:3L) {
hole <- holes[[i]]
holeType <- hole[["type"]]
if(holeType == "ccc") {
x <- with(hole, innerSoddyCircle(c1, c2, c3, color = color))
} else if(holeType == "ccl") {
x <- with(hole, side.circle.circle(A, B, cA, cB, color = color))
} else if (holeType == "cll") {
x <- with(hole, side.side.circle(A, B, C, circle, color = color))
}
newholes <- c(newholes, list(x))
}
newholes
}
MalfattiCircles <- function(A, B, C) {
Triangle$new(A, B, C)$MalfattiCircles()
}
drawTriangularGasket <- function(mcircles, A, B, C, colors, depth) {
C1 <- mcircles[[1L]]
C2 <- mcircles[[2L]]
C3 <- mcircles[[3L]]
triangles3d(cbind(rbind(A, B, C), 0), col = "yellow", alpha = 0.2)
holes <- list(
side.circle.circle(A, B, C1, C2, color = colors[1L]),
side.circle.circle(B, C, C2, C3, color = colors[1L]),
side.circle.circle(C, A, C3, C1, color = colors[1L]),
innerSoddyCircle(C1, C2, C3, color = colors[1L]),
side.side.circle(A, B, C, C1, color = colors[1L]),
side.side.circle(B, A, C, C2, color = colors[1L]),
side.side.circle(C, A, B, C3, color = colors[1L])
)
for(d in 1L:depth) {
n_holes <- length(holes)
Holes <- list()
for(i in 1L:n_holes) {
Holes <- append(Holes, Newholes(holes[[i]], colors[d + 1L]))
}
holes <- do.call(list, Holes)
}
}
drawCircularGasket <- function(c0, n, phi, shift, depth, colors) {
ApollonianCircles <- ApollonianGasket(c0, n, phi, shift, depth)
for(i in 1:depth) {
for(circ in ApollonianCircles[[i]]){
drawSphere(circ, color = colors[i])
}
}
}
library(viridisLite)
A <- c(-5, -4)
B <- c(5, -2)
C <- c(0, 6)
mcircles <- MalfattiCircles(A, B, C)
depth <- 3L
colors <- viridis(depth + 1L)
n1 <- 3L
n2 <- 4L
n3 <- 5L
depth2 <- 3L
phi1 <- 0.2
phi2 <- 0.3
phi3 <- 0.4
shift <- 0
colors2 <- plasma(depth2)
## ----rgl, eval=FALSE----------------------------------------------------------
# library(rgl)
# open3d(windowRect = c(50, 50, 562, 562), zoom = 0.9)
# bg3d(rgb(54, 57, 64, maxColorValue = 255))
# drawTriangularGasket(mcircles, A, B, C, colors, depth)
# drawCircularGasket(mcircles[[1L]], n1, phi1, shift, depth2, colors2)
# drawCircularGasket(mcircles[[2L]], n2, phi2, shift, depth2, colors2)
# drawCircularGasket(mcircles[[3L]], n3, phi3, shift, depth2, colors2)
## ----malfattiCircular, warning=FALSE, eval=require("rgl", quietly=TRUE)-------
library(rgl)
iteration <- function(circlesWithIndicator, inversions) {
out <- list()
for(j in seq_along(circlesWithIndicator)) {
circle <- circlesWithIndicator[[j]][["circle"]]
indic <- circlesWithIndicator[[j]][["indic"]]
for(i in 1L:4L) {
if(i != indic) {
circleWithIndicator <- list(
"circle" = inversions[[i]]$invertCircle(circle),
"indic" = i
)
out <- append(out, list(circleWithIndicator))
}
}
}
out
}
gasket <- function(circlesWithIndicator, inversions, depth, colors) {
if(depth > 0){
circlesWithIndicator <- iteration(circlesWithIndicator, inversions)
for(i in seq_along(circlesWithIndicator)) {
drawSphere(circlesWithIndicator[[i]]$circle, color = colors[1L])
}
colors <- colors[-1L]
gasket(circlesWithIndicator, inversions, depth-1L, colors)
}
}
drawGasket <- function(triangle, depth, colors) {
Mcircles <- triangle$MalfattiCircles()
Mtriangle <- Triangle$new(
Mcircles[[1L]]$center, Mcircles[[2L]]$center, Mcircles[[3L]]$center
)
soddyO <- Mtriangle$outerSoddyCircle()
Mcircles <- append(Mcircles, list(soddyO))
for(i in 1L:4L) {
lines3d(
cbind(Mcircles[[i]]$asEllipse()$path(), 0),
color = "black", lwd = 2
)
}
inversions <- vector("list", 4L)
circlesWithIndicator <- vector("list", 4L)
inversions[[1L]] <- inversionFixingThreeCircles(
soddyO, Mcircles[[2L]], Mcircles[[3L]]
)
inversions[[2L]] <- inversionFixingThreeCircles(
soddyO, Mcircles[[1L]], Mcircles[[3L]]
)
inversions[[3L]] <- inversionFixingThreeCircles(
soddyO, Mcircles[[1L]], Mcircles[[2L]]
)
inversions[[4L]] <- inversionFixingThreeCircles(
Mcircles[[1L]], Mcircles[[2L]], Mcircles[[3L]]
)
for(i in 1L:4L) {
circlesWithIndicator[[i]] <-
list("circle" = inversions[[i]]$invertCircle(Mcircles[[i]]), "indic" = i)
drawSphere(circlesWithIndicator[[i]]$circle, color = colors[1L])
}
colors <- colors[-1L]
gasket(circlesWithIndicator, inversions, depth, colors)
}
CircularMalfattiGasket <- function(C, depth, colors) {
A <- c(0,0); B <- c(1,0)
t <- Triangle$new(A, B, C)
Mcircles <- t$MalfattiCircles()
Mtriangle <- Triangle$new(
Mcircles[[1L]]$center, Mcircles[[2L]]$center, Mcircles[[3L]]$center
)
soddyO <- Mtriangle$outerSoddyCircle()
center <- soddyO$center; radius = -soddyO$radius
A1 <- (A-center)/radius; B1 <- (B-center)/radius; C1 <- (C-center)/radius;
t1 <- Triangle$new(A1, B1, C1);
drawGasket(t1, depth, colors)
}
## ----plotMalfattiCircularGasket, eval=FALSE-----------------------------------
# open3d(windowRect = 50 + c(0, 0, 900, 300))
# mfrow3d(1, 3)
# view3d(0, 0, zoom = 0.7)
# CircularMalfattiGasket(
# C = c(0, sqrt(3/2)), depth = 2L,
# colors = c("yellow", "orangered", "darkmagenta")
# )
# next3d()
# view3d(0, 0, zoom = 0.7)
# CircularMalfattiGasket(
# C = c(1, sqrt(3/2)), depth = 2L,
# colors = c("yellow", "orangered", "darkmagenta")
# )
# next3d()
# view3d(0, 0, zoom = 0.7)
# CircularMalfattiGasket(
# C = c(2, sqrt(3/2)), depth = 2L,
# colors = c("yellow", "orangered", "darkmagenta")
# )
## ----hyperbola, fig.width=7, fig.height=5-------------------------------------
# take a hyperbola
L1 <- LineFromInterceptAndSlope(0, 2) # asymptote 1
L2 <- LineFromInterceptAndSlope(-2, -0.15) # asymptote 2
M <- c(2, 3) # a point on the hyperbola
hyperbola <- Hyperbola$new(L1, L2, M)
# take a point on the hyperbola and the tangent at this point
OAB <- hyperbola$OAB()
O <- OAB$O; A <- OAB$A; B <- OAB$B
t <- 0.1
P <- O + cosh(t)*A + sinh(t)*B
tgt <- Line$new(P, P + sinh(t)*A + cosh(t)*B)
# the triangle of interest
C <- intersectionLineLine(L1, tgt)
D <- intersectionLineLine(L2, tgt)
trgl <- Triangle$new(O, C, D)
# plot
opar <- par(mar = c(4, 4, 1, 1))
hyperbola$plot(lwd = 2)
draw(L1, col = "red")
draw(L2, col = "red")
text(t(O), "O", pos = 3)
points(t(P), pch = 19, col = "blue")
text(t(P), "P", pos = 4)
draw(tgt, col = "blue", lwd = 2)
text(t(C), "C", pos = 2)
text(t(D), "D", pos = 4)
trgl$plot(add = TRUE, col = "yellow")
par(opar)
# theorem checking: area of the triangle does not depend on
# the choice of P; more precisely, it is equal to ab
trgl$area()
with(hyperbola$abce(), a * b)
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.