#     Rotated sinc function.

     x <- seq(-10, 10, length=50)
     y <- x
     f <- function(x,y)
     {
         r <- sqrt(x^2+y^2)
         10 * sin(r)/r
     }
     z <- outer(x, y, f)
     z[is.na(z)] <- 1
     zlims <- range(z, na.rm = TRUE)

     par(bg = "white")
     persp(x, y, z, theta = 30, phi = 30, expand = 0.5, col = "lightblue",
           xlab = "X", ylab = "Y", zlab = "Z")


     persp(x, y, z, theta = 30, phi = 30, expand = 0.5, col = "lightblue",
           ltheta = 120, shade = 0.75, ticktype = "detailed",
           xlab = "X", ylab = "Y", zlab = "Z")

     par(new=TRUE)
     f <- function(x,y)
     {
         r <- sqrt(x^2+y^2)
         sin(r)
     }
     z <- outer(x, y, f)
     z[is.na(z)] <- 1
     persp(x, y, z, theta = 30, phi = 30, expand = 0.5, col = "green", zlim=zlims)

     #     sin(x+y)
     f <- function(x,y)
     {
	sin(x+y)
     }
     z <- outer(x, y, f)
     z[is.na(z)] <- 1
     par(bg = "white")
     persp(x, y, z, theta = 30, phi = 30, expand = 0.5, col = "lightblue",
           xlab = "X", ylab = "Y", zlab = "Z")
     persp(x, y, z, theta = 30, phi = 30, expand = 0.5, col = "lightblue",
           ltheta = 120, shade = 0.75, ticktype = "detailed",
           xlab = "X", ylab = "Y", zlab = "Z")

     par(new=TRUE)
     f <- function(x,y)
     {
	cos(x+y)
     }
     z <- outer(x, y, f)
     z[is.na(z)] <- 1
     persp(x, y, z, theta = 30, phi = 30, expand = 0.5, col = "green",
           xlab = "X", ylab = "Y", zlab = "Z")

    # (2) Visualizing a simple DEM model

     data(volcano)
     z <- 2 * volcano        # Exaggerate the relief
     x <- 10 * (1:nrow(z))   # 10 meter spacing (S to N)
     y <- 10 * (1:ncol(z))   # 10 meter spacing (E to W)
     persp(x, y, z, theta = 120, phi = 15, scale = FALSE, axes = FALSE)

     # (3) Now something more complex
     #     We border the surface, to make it more "slice like"
     #     and color the top and sides of the surface differently.

     zmin <- min(z) - 20
     z <- rbind(zmin, cbind(zmin, z, zmin), zmin)
     x <- c(min(x) - 1e-10, x, max(x) + 1e-10)
     y <- c(min(y) - 1e-10, y, max(y) + 1e-10)

     fill <- matrix("green3", nr = nrow(z)-1, nc = ncol(z)-1)
     fill[,1] <- "gray"
     fill[,ncol(fill)] <- "gray"
     fill[1,] <- "gray"
     fill[nrow(fill),] <- "gray"

     par(bg = "lightblue")
     persp(x, y, z, theta = 120, phi = 15, col = fill, scale = FALSE, axes =
     FALSE)
     title(main = "Maunga Whau\nOne of 50 Volcanoes in the Auckland Region.",
           font.main = 4)

     par(bg = "slategray")
     persp(x, y, z, theta = 135, phi = 30, col = fill, scale = FALSE,
           ltheta = -120, lphi = 15, shade = 0.65, axes = FALSE)
     persp(x, y, z, theta = 135, phi = 30, col = "green3", scale = FALSE,
           ltheta = -120, shade = 0.75, border = NA, box = FALSE)

# http://www.mail-archive.com/r-help@stat.math.ethz.ch/msg13591.html
pmat <- persp(0:1, 0:1, matrix(,2,2), xlim=c(-1,1), ylim=c(-1,1), zlim=c(-1,1),
              theta=25, phi=30, expand=.9, xlab="X", ylab="Y", zlab="Z")

trans3d <- function(x,y,z, pmat) { # From the help for "persp"
  tr <- cbind(x,y,z,1) %*% pmat
  list(x = tr[,1]/tr[,4], y= tr[,2]/tr[,4])
}

theta <- seq(0, 2*pi, length=51)
phi <- seq(0, pi, length=26)
x <- cos(theta) %o% sin(phi)
y <- sin(theta) %o% sin(phi)
z <- rep(1, length(theta)) %o% cos(phi)

for (j in seq(phi)[-1]) for (i in seq(theta)[-1]) {
  idx <- rbind(c(i-1,j-1), c(i,j-1), c(i,j), c(i-1,j))
  polygon(trans3d(x[idx], y[idx], z[idx], pmat))
}