---
title: "Examples of the PlaneGeometry package"
output: 
  rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{Examples of the PlaneGeometry package}
  %\VignetteEncoding{UTF-8}
  %\VignetteEngine{knitr::rmarkdown}
editor_options: 
  chunk_output_type: console
---

```{r, include = FALSE}
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>"
)
```

```{r setup}
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`.

```{r}
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.

```{r 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)
```


# 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?

```{r}
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()
```

```{r 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)
```


# 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`.

```{r}
# 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:

```{r 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)
```

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:

```{r}
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?

```{r 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)
```


# 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.

```{r}
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()
}
```


```{r 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)
)
```

Here is a version which allows to fill the hyperbolic triangles:

```{r}
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()
}
```


```{r 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)
)
```


# Director circle of an ellipse 

Let's draw the [director circle](https://en.wikipedia.org/wiki/Director_circle) 
of an ellipse. We start by constructing the minimum bounding box of the ellipse.

```{r}
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.

```{r, 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)
```

```{r 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)
```


# 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.

```{r}
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)
```

```{r 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)
```

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

```{r, 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
)
```


![](https://github.com/stla/PlaneGeometry/raw/master/inst/imgs/SteinerChainElliptical.gif){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.

```{r}
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)
})
```

```{r 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)
```

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.

```{r 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)
```

Run the code below to see an animated trajectory:

```{r, 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)
```


# 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:

```{r 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]]
)
```

There are 425 generalized circles:

```{r}
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:

```{r 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]
}
```

For example:

```{r}
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:

```{r}
gcircles <- uniqueWith(
  gcircles,
  function(g1, g2){
    class(g1)[1L] == class(g2)[1L] && g1$isEqual(g2)
  }
)
```

Now it remains 161 generalized circles:

```{r}
length(gcircles)
```

Let's write a helper function to draw these generalized circles:

```{r 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], ...)
  }
}
```

And now let's draw them:

```{r 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 circles

This construction is taken from the book *Indra's Pearls: The Vision of Felix Klein*.

```{r 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]]])
}
```

Here is the picture:

```{r 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)
```

Here is the same picture but with better quality:

![](https://github.com/stla/PlaneGeometry/raw/master/inst/imgs/SchottkyCircles_tikz.png){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](https://www3.risc.jku.at/publications/download/risc_5011/DiplomaThesisPonweiser.pdf) written by Thomas Ponweiser. This is the paper 
which motivated me to implement the generalized power of a Möbius transformation.


```{r 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)
    }
  }
}
```

To get the animation, run:

```{r, 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
)
```

![](https://github.com/stla/PlaneGeometry/raw/master/inst/imgs/modularTessellation_compressed.gif){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.

```{r 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
}
```

Let's apply our function:

```{r 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)
```

We can do an animation now:

```{r 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
)
```

![](https://github.com/stla/PlaneGeometry/raw/master/inst/imgs/ApollonianGasket_compressed.gif){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

```{r}
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.

```{r 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)
  }
}
```

```{r, eval=FALSE}
library(gifski)
save_gif(
  fanim(),
  "ApollonianMobius.gif",
  width = 512, height = 512,
  delay = 0.1
)
par(opar)
```

![](https://github.com/stla/PlaneGeometry/raw/master/inst/imgs/ApollonianMobius.gif){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:

```{r 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)
```

Let's plot the fractal in 3D with the help of the **rgl** package.

```{r 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), ...)
}
```

Now here is how to plot the fractal and make an animation:

```{r 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
)
```

![](https://github.com/stla/PlaneGeometry/raw/master/inst/imgs/Apollony_compressed.gif){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.

```{r 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)
```

And we get the 3D picture:

```{r 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)
```

![](https://github.com/stla/PlaneGeometry/raw/master/inst/imgs/MalfattiGasket.png){width=50%}

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


# Circular Malfatti gasket

Yet another Malfatti based Apollonian gasket. This one uses the outer 
Soddy circle, which has a negative radius!

```{r 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)
}
```

```{r 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")
)
```


![](https://github.com/stla/PlaneGeometry/raw/master/inst/imgs/MalfattiCircularGasket.png){width=80%}


# Hyperbola

Finally, in order to illustrate the `Hyperbola` objects, we will check the 
so-called "triangle tangent-asymptotes" theorem. Observe the figure below. 
The point `P` has been taken arbitrarily on the hyperbola, and the blue line 
is the tangent at `P`.

```{r 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)
```

The theorem claims that the area of the yellow triangle does not depend on the 
location of `P`, and this area is equal to the product of the two semi-axes of 
the hyperbola.