# Bookstein Orange Book page 27
require(sm)
require(scatterplot3d)
n<-50
x<- (-n):n
x<- x/(n/2)
y<- x
 r<- matrix(0, nrow=2*n, ncol=2*n)
 z<- r
 mx<- r
 my<- r
 for (i in 1:(2*n)) {for (j in 1:(2*n))
   {r[i,j]<- (x[i]^2 + x[j]^2)^0.5
    z[i,j]<- -(r[i,j]^2)*log(r[i,j]^2)
    mx[i,j]<- x[i]
    my[i,j]<- y[j]}}
  image(z, col = heat.colors(2*n))
  contour(z, ad = TRUE)
  title("Bookstein orange book, Fig 2.2.1", font = 4)
pause()
# col = terrain.colors
  i<-(n/2):(2*n-n/2)
  image(z[i,i], col = terrain.colors(2*n))
  contour(z[i,i], add = TRUE)
  title("Bookstein orange book, Fig 2.2.1", font = 4)

pause()
# color <- rep("green", length(z))
 color <- "red"
 mz<- z
 x<- as.vector(mx[i,i])
 y<- as.vector(my[i,i])
 z<- as.vector(z[i,i])

 s3d<- scatterplot3d(x, y, z, main="Bookstein Fig 2.2.1", type="l", grid=TRUE, color=0)
 for (j in i) {
 s3d$points3d(mx[i,j], my[i,j], mz[i,j], type="l", col="red", lwd=1)
              }
 for (j in i) {
 s3d$points3d(mx[j,i], my[j,i], mz[j,i], type="l", col="green", lwd=1)
              }
              
pause()
# require(Rcmdr)
# s3d<- scatter3d(x, y, z, main="Bookstein Fig 2.2.1", type="l", grid=TRUE, lwd=1)
div<-20
lim<-50
del<-lim/div
nm<- 1+2*div
x<- seq(-lim,lim,del)
y<- seq(-lim,lim,del)
z<- matrix(0,nm,nm)
nxyz<- nm*nm
xyz<- matrix(0,nxyz,3)
for (i in 1:nm) {
    xi<-x[i]
    for (j in 1:nm) {
      yj<- y[j]
      z[i,j]<- 4*(yj^2 - xi^2)/(1+xi^2+yj^2)
      xyz[j+(i-1)*nm,]<- c(xi, yj, z[i,j])
                    }
                }
                
 xp<- as.vector(xyz[,1])
 yp<- as.vector(xyz[,2])
 zp<- as.vector(xyz[,3])

 s3d<- scatterplot3d(xp, yp, zp, main="Bookstein Fig 2.2.2", color=0, 
                        angle = 30, grid=TRUE)
 for(j in 1:nm) {
 inde<- nm*(j-1)+ (1:nm)
 s3d$points3d(xp[inde], yp[inde], zp[inde], type="l", col="red", lwd=1)
                }
 for(j in 1:nm) {
 inde<- (j-1) + seq(1,nxyz,nm)
 s3d$points3d(xp[inde], yp[inde], zp[inde], type="l", col="green", lwd=1)
                }
              
pause()

 xp<- as.vector(xyz[,1])
 yp<- as.vector(xyz[,2]+xyz[,3])

 plot(xp, yp, main="Bookstein Fig 2.2.2", typ='n')
 for(j in 1:nm) {
 inde<- nm*(j-1)+ (1:nm)
 lines(xp[inde], yp[inde], col="red", lwd=1)
                }
 for(j in 1:nm) {
 inde<- (j-1) + seq(1,nxyz,nm)
 lines(xp[inde], yp[inde], col="green", lwd=1)
                }