library(pracma)
phi <- seq(-pi/2, pi/2, by=0.1)
theta <- seq(0, 2*pi, by=0.1)
mgrd <- meshgrid(phi, theta)
phi <- mgrd$X
theta <-  mgrd$Y
scale = qnorm(1-0.01/2)
x <- cos(phi) * cos(theta) * scale
dim(x) <- NULL
y <- cos(phi) * sin(theta) * scale
dim(y) <- NULL
z <- sin(phi) * scale
dim(z) <- NULL

sphere = as.matrix(rbind(x, y, z))
N <- 100
x <- rnorm(N)
y <- rnorm(N)
z <- rnorm(N)


C <- t(matrix(c(4, 0, 0, 0, 2, 0, 0, 0, 1), nrow = 3))

theta <- pi/4
phi <- pi/6

rot1 <- t(matrix(c(sin(theta), cos(theta), 0,
                   cos(theta), -sin(theta), 0,
                   0, 0, 1), nrow=3))

rot2 <- t(matrix(c(1, 0, 0,
                   0, sin(phi), cos(phi),
                   0, cos(phi), -sin(phi)),nrow=3))

rot = rot1 %*% rot2
C <- rot %*% C

points <- setNames(data.frame(t(C %*% as.matrix(rbind(x, y, z)))), c('x', 'y', 'z'))
ell <- setNames(data.frame(t(C %*% sphere)), c('x', 'y', 'z'))
print(points[1:5,])

##            x          y           z
## 1  0.6342832  3.3748724 -0.03145166
## 2 -1.7891032 -2.3421015  0.52137608
## 3  1.2805058  0.6694286  1.17129645
## 4  5.6448206  4.4006883 -0.15911536
## 5  2.2407650  1.5251863 -0.62399069

library(plotly)

## Loading required package: ggplot2
## 
## Attaching package: 'plotly'
## 
## The following object is masked from 'package:ggplot2':
## 
##     last_plot
## 
## The following object is masked from 'package:graphics':
## 
##     layout

p <- plot_ly(points, x = x, y = y, z = z, type = "scatter3d", mode = "markers") %>% layout(scene = list(aspectmode = 'data'))

p 

p <- add_trace(p, data = ell, x=x, y=y, z=z, type='mesh3d', alphahull = 0, opacity=0.4, hoverinfo='none') 
p

arrow_profile <- function(len, width=0.2, headlength=2, headwidth=0.5, steps=50) {
  #returns list of X's and Y's
  #Y=f(X)
  # the graph of Y is an arrow profile
  X <- c(-len, seq(-len, len-headlength, length.out = steps %/% 2),
         len-headlength,
         seq(len-headlength+1E-5, len, length.out = steps %/% 2))
  Y <- c(0, rep(width, length(X)-1))
  for (i in 1:length(X)) {
    if (X[i] > len-headlength) {
      Y[i] = (len - X[i])*headwidth/(headlength)
    }
  }
  list(X=X, Y=Y)
}

rotate <- function(profile, thetasteps=50) {
  # rotate profile Y = f(X) around X axes
  X <- profile$X
  R <- profile$Y
  Theta <- seq(0, 2*pi+2*pi/thetasteps, length.out = thetasteps)
  mgrd <- meshgrid(X, Theta)
  X <- mgrd$X
  R <- meshgrid(R, Theta)$X
  Theta <- mgrd$Y
  Y <- R*cos(Theta)
  Z <- R*sin(Theta)
  list(x=X, y=Y, z=Z)
}

apply_matrix <- function(C, points) {
  x <- points$x
  y <- points$y
  z <- points$z
  n <- nrow(x)
  m <- ncol(x)
  outx <- matrix(rep(0, n*m), nrow = n)
  outy <- matrix(rep(0, n*m), nrow = n)
  outz <- matrix(rep(0, n*m), nrow = n)
  
  for (i in 1:nrow(x)) {
    for (j in 1:ncol(x)) {
      out <- C %*% c(x[i, j], y[i, j], z[i, j])
      outx[i,j] <- out[1,]
      outy[i,j] <- out[2,]
      outz[i,j] <- out[3,]
    }
  }
  list(x=outx,y=outy,z=outz)
}

p <- add_trace(p, data = apply_matrix(rot,rotate(arrow_profile(14))), x=x, y=y, z=z, type='surface', hoverinfo='none', showscale=F, colorscale=list(list(0, "#009E73"), list(1, "#009E73")))
p

scale = 2
project_along_PC1 <- list(camera=list(eye=list(x=rot[1,1]*scale, y=rot[2,1]*scale, z=rot[3,1]*scale)))
                       
                       
p %>% layout(scene=project_along_PC1)

p

p %>% layout(
  scene=list(
    camera=list(
      eye=list(
        x=rot[1,3]*scale,
        y=rot[2,3]*scale,
        z=rot[3,3]*scale),
      up=list(
        x=rot[1,2],
        y=rot[2,2],
        z=rot[3,2])
      )
    )
  )