### Demos for  persp()  plots   -- things not in  example(persp)
### -------------------------

require(graphics); require(datasets)
## make sure a device is open
if(dev.cur() <= 1) get(getOption("device"))()
## if(dev.interactive())
## is FALSE in demo() {using source()}:
## Fix it later, copy its code for now
is.dev.interactive <- .Device %in% c("X11", "GTK", "gnome",
				"quartz", "windows", "JavaGD")
op <- par(ask = is.dev.interactive)

## (1) The Obligatory Mathematical surface.
##     Rotated sinc function.

x <- seq(-10, 10, length = 50)
y <- x
rotsinc <- function(x,y)
{
    sinc <- function(x) { y <- sin(x)/x ; y[is.na(y)] <- 1; y }
    10 * sinc( sqrt(x^2+y^2) )
}
sinc.exp <- expression(z == Sinc(sqrt(x^2 + y^2)))

z <- outer(x, y, rotsinc)


par(bg = "white")
persp(x, y, z, theta = 30, phi = 30, expand = 0.5, col = "lightblue")
title(sub=".")## work around persp+plotmath bug
title(main = sinc.exp)

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")
title(sub=".")## work around persp+plotmath bug
title(main = sinc.exp)



## (2) Visualizing a simple DEM model

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.

z0 <- min(z) - 20
z <- rbind(z0, cbind(z0, z, z0), z0)
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[ , i2 <- c(1,ncol(fill))] <- "gray"
fill[i1 <- c(1,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)
## Don't draw the grid lines :  border = NA
persp(x, y, z, theta = 135, phi = 30, col = "green3", scale = FALSE,
      ltheta = -120, shade = 0.75, border = NA, box = FALSE)

## `color gradient in the soil' :
fcol <- fill ; fcol[] <- terrain.colors(nrow(fcol))
persp(x, y, z, theta = 135, phi = 30, col = fcol, scale = FALSE,
      ltheta = -120, shade = 0.3, border = NA, box = FALSE)

## `image like' colors on top :
fcol <- fill
zi <- volcano[ -1,-1] + volcano[ -1,-61] +
           volcano[-87,-1] + volcano[-87,-61]  ## / 4
fcol[-i1,-i2] <-
    terrain.colors(20)[cut(zi, quantile(zi, seq(0,1, len = 21)),
                           include.lowest = TRUE)]
persp(x, y, 2*z, theta = 110, phi = 40, col = fcol, scale = FALSE,
      ltheta = -120, shade = 0.4, border = NA, box = FALSE)


## reset par():
par(op)