library(plot3D)

arrows3D(x0 = runif(10), y0 = runif(10), z0 = runif(10),
x1 = runif(10), y1 = runif(10), z1 = runif(10),
colvar = 1:10, code = 1:3, main = "arrows3D", colkey = FALSE)

x <- c(0,0,0.8)
y <- c(0,0.8,0)
z <- c(0,0,0)

polygon3D(x, y, z, border = "black", lwd = 1, alpha = 0.8, theta=230,
           bty='n',main = "polygon3D",xlim=range(-.2,1),ylim=range(-.2,1))

th<-30
O<-c(0,0,0)
Y<-c(0,3,2)
Yhat<-c(0,3,0)
fm1<-rbind(O,Y,Yhat)
Y1<-c(0.5,2.5,0)
fm2<-rbind(O,Y1,Yhat)
fm3<-rbind(Y,Y1,Yhat)
fm4<-rbind(O,Y1,Y)


fm1eps<-rbind(fm1,fm2[2,])
fm1eps[1,2]<-fm1[1,2]-0.35

fm1eps[3,2]<-fm1[2,3]-0.2

fm1eps[2,3]<-fm1[2,3]+0.4
fm1eps[3,1]<-fm1[3,1]+0.55
fm1eps[4,]<-fm2[3,]-c(-0.1,0.5,0.1)

polygon3D(fm1[1,1],fm1[1,2],fm1[1,3], border = "black", lwd = 1, alpha = 0.8, theta=60,
           bty='b2',main = "",xlim=range(-0.25,0.65),ylim=range(-.45,4),zlim=range(0,5),
	     xlab="",ylab="",zlab=expression(y))
points3D(fm1[1,1],fm1[1,2],fm1[1,3],add=T,pch=19,colkey=FALSE)
text3D(fm1eps[1,1],fm1eps[1,2],fm1eps[1,3],c('O'),add=T)
dev.print(device=pdf,file='OrthPred1.pdf',width=8,height=8)

polygon3D(fm1[,1],fm1[,2],fm1[,3], border = "black", lwd = 1, alpha = 0.8, theta=60,
           bty='b2',main = "",xlim=range(-0.25,0.65),ylim=range(-.45,4),zlim=range(0,5),
	     xlab="",ylab="",zlab=expression(y),axes=FALSE)
segments3D(0,0,0,0,1,0,lty=1,add=T)
text3D(fm1eps[c(1,2,4),1],fm1eps[c(1,2,4),2],fm1eps[c(1,2,4),3],c('O',expression(Y),expression(hat(Y))),add=T)
points3D(fm1[,1],fm1[,2],fm1[,3],add=T,pch=19,colkey=FALSE)
dev.print(device=pdf,file='OrthPred2.pdf',width=8,height=8)

polygon3D(fm4[,1],fm4[,2],fm4[,3], border = "black", col='green',lwd = 1, alpha = 0.2,theta=60,
           bty='b2',main = "",xlim=range(-0.25,0.65),ylim=range(-.45,4),zlim=range(0,5),
	     xlab="",ylab="",zlab=expression(y),axes=FALSE)
segments3D(Y[1],Y[2],Y[3],Yhat[1],Yhat[2],Yhat[3],lty=2,add=T)
segments3D(0,0,0,Yhat[1],Yhat[2],Yhat[3],lty=2,add=T)
segments3D(Y1[1],Y1[2],Y1[3],Yhat[1],Yhat[2],Yhat[3],lty=2,add=T)
arrows3D(0,0,0,1.6*Y1[1],1.6*Y1[2],1.6*Y1[3],lty=2,add=T)
text3D(1.6*Y1[1],1.6*Y1[2],1.6*Y1[3],expression(X^{(1)}),add=T)
points3D(fm4[,1],fm4[,2],fm4[,3],add=T,pch=19,colkey=FALSE)
text3D(fm1eps[c(1,2,3),1],fm1eps[c(1,2,3),2],fm1eps[c(1,2,3),3],c('O',expression(Y),expression(hat(Y)^{(1)})),add=T)
dev.print(device=pdf,file='OrthPred3.pdf',width=8,height=8)


polygon3D(fm1[,1],fm1[,2],fm1[,3], border = "black", lwd = 1, alpha = 0.8, theta=60,
           bty='b2',main = "",xlim=range(-0.25,0.65),ylim=range(-.45,4),zlim=range(0,5),
	     xlab="",ylab="",zlab=expression(y),axes=FALSE)
text3D(fm1eps[c(1,2,4),1],fm1eps[c(1,2,4),2],fm1eps[c(1,2,4),3],c('O',expression(Y),expression(hat(Y))),add=T)
points3D(fm1[,1],fm1[,2],fm1[,3],add=T,pch=19,colkey=FALSE)
points3D(fm4[,1],fm4[,2],fm4[,3],add=T,pch=19,colkey=FALSE)
dev.print(device=pdf,file='OrthPred4.pdf',width=8,height=8)

polygon3D(fm3[,1],fm3[,2],fm3[,3], border = "black", col='blue',lwd = 1, alpha = 0.5, add=T)
text3D(fm1eps[c(1,2,3),1],fm1eps[c(1,2,3),2],fm1eps[c(1,2,3),3],c('O',expression(Y),expression(hat(Y)^{(1)})),add=T)
arrows3D(0,0,0,1.6*Y1[1],1.6*Y1[2],1.6*Y1[3],lty=2,add=T)
text3D(1.6*Y1[1],1.6*Y1[2],1.6*Y1[3],expression(X^{(1)}),add=T)
dev.print(device=pdf,file='OrthPred5.pdf',width=8,height=8)

polygon3D(fm2[,1],fm2[,2],fm2[,3], border = "black", col='red',lwd = 1, alpha = 0.2, add=T)
text3D(fm1eps[,1],fm1eps[,2],fm1eps[,3],c('O',expression(Y),expression(hat(Y)^{(1)}),expression(hat(Y))),add=T)
dev.print(device=pdf,file='OrthPred6.pdf',width=8,height=8)


##########################################################################
X       <- seq(-0.25, 1, length.out = 50)
Y       <- seq(-0.25, 4, length.out = 50)
M       <- mesh(X, Y)
phi     <- M$x
theta   <- M$y

# x, y and z grids
r <- sin(4*phi)^3 + cos(2*phi)^3 + sin(6*theta)^2 + cos(6*theta)^4
x <- phi
y <- theta
z <- 0*r * sin(phi) * sin(theta)-0.1
surf3D(x, y, z, col='gray', shade = 0.5, colkey=FALSE,box = FALSE,add=T)