library(fields)

## persp example code
par(bg = "ivory")
n=3
switch(n, 
  {f <-  function(a, b) 2*exp(-(a^2+b^2))
   fm<- "z= 2*exp(-(x^2+y^2))"},
  {f <-  function(a, b) sqrt(a^2+b^2)
   fm<- "z= sqrt(x^2+y^2)"},
  {f <-  function(a, b) a^2-b^2
   fm<- "z= x^2-y^2"},
  {f <-  function(a, b) -3/sqrt((a^2+b^2))
   fm<- "z= -/sqrt(-(x^2+y^2)"},
 {f <-  function(a, b) sin((a^2+b^2)/1.5)
   fm<- "z= sin(x^2+y^2)/1.5"}



)

x <- seq(-3, 3, length = 60)
y <- seq(-3, 3, length = 60)
z <- outer(x, y, f)
nrz <- nrow(z)
ncz <- ncol(z)
# Create a function interpolating colors in the range of specified colors
jet.colors <- colorRampPalette( c("black","blue", "green","cyan","yellow","orange","red") )
# Generate the desired number of colors from this palette
nbcol <- 100
color <- jet.colors(nbcol)
# Compute the z-value at the facet centres
zfacet <- (z[-1, -1] + z[-1, -ncz] + z[-nrz, -1] + z[-nrz, -ncz])/4
# Recode facet z-values into color indices
facetcol <- cut(zfacet, nbcol)
persp(x, y, z, col = color[facetcol], theta = 60, phi = 20, axes=T, ticktype='detailed', border=NA)

## add color bar
image.plot(legend.only=T, zlim=range(zfacet), col=color)
#mtext(fm,2,0,col=4)
title(fm)

# 拡大縮小回転などのマウス操作ができる
library(rgl)
persp3d(x, y, z,  expand = 1.0, col = rainbow(20))

 #########################
par3d(windowRect=c(200,200,600,600));
rgl.bg(color=c("lightcyan"))
rgl.spheres(0,0,0,0.8,color=rainbow(20))
rgl.spheres(rnorm(15), rnorm(15), rnorm(16), radius=0.4, color=rainbow(16))
rgl.bbox(color="#aa55cc")
#########################
  rgl.bg(color=c("lightcyan"))
  x <- sort(rnorm(1000))
  y <- rnorm(1000)
  z <- rnorm(1000) 
   par3d(windowRect=c(100,100,600,600))  
  plot3d(x, y, z, type="s",col=rainbow(1000))

######## 3次元空間内曲線　　x=x(t)  y=y(t)  z=z(t)
  library(scatterplot3d)
  library(rgl)
  t <- seq(0, 4*pi, 0.01)  ;s <- seq(0, 4*pi, 0.01)
  #x <- sin(t);  y <- cos(t);  z<- t   　　　　　　　　　 # 螺旋
  #x <- sin(t);  y <- sin(2*t);  z<- sin(3*t)    # リサージュ
  x <- t-sin(t);  y <- 1-cos(t);  z<- t/5        #　サイクロイド

  par(cex=0.1)
  scatterplot3d(x, y, z, highlight.3d=TRUE, col.axis="blue",
                col.grid="lightblue", main="scatterplot3d - 1", pch=20)
  pp <- recordPlot()
  #png("cycloid0.png", width = 400, height = 400)
  replayPlot(pp)
  dev.off()
par3d(windowRect=c(100,100,700,700))
plot3d(x, y, z,  expand = 1.0, col = 4)