outer_product <- function (a, b) {
  c(a[2]*b[3] - a[3]*b[2], a[3]*b[1] - a[1]*b[3], a[1]*b[2] - a[2]*b[1])
}

d <- function(t) {
  s <- t/sqrt(3)
  h <- sqrt(1 - ((1 + sqrt(3)/2)^2 + (1/2)^2)*s^2)

  A <- c(s, 0, 0)
  B <- c(-1/2*s, sqrt(3)/2*s, 0)
  A_ <- c(-sqrt(3)/2*s, 1/2*s, h)
  B_ <- c(0, -s, h)

  AA_ <- A_ - A
  BB_ <- B_ - B
  BA <- A - B

  AA_xBB_ <- outer_product(AA_, BB_)
  ((A - B) %*% AA_xBB_) / sqrt(AA_xBB_ %*% AA_xBB_)
}

ts <- (1:89)/100
ds <- c()
for (i in 1:89) {
  ds[i] <- d(ts[i])
}
plot(ts, ds, type="l")