R/Perspective Plots/Code2

# making overlapping surfaces in R using persp
# by Eli Gurarie for Ian Taylor, who deserves it, but not without minor reservations 
# that it is James Murphy who might, eventually, benefit.

# setup some basic stuff
x = seq(0,.2,length=30)
y = seq(0,20,length=25)

# create 3 surfaces pass through each other overlap
z1 = matrix(seq(0.1,14,length=30),30,25)
z2 = matrix(seq(4,6.2,length=30),30,25)
z3 = matrix(seq(0,8,length=30),30,25)

alpha <- 0.6
mycols <- c(rgb(1,0,0,alpha), rgb(0,1,0,alpha),  rgb(0,0,1,alpha))

p1 <- persp(x,y,z1,phi=25,theta=-45,ltheta=120,
		lphi=0,col=rgb(1,0,0,0),zlim=c(0,14),zlab='z',ticktype='detailed', bor="white")

for(level in c("bottom", "middle", "top"))
for(k in 1:3)
for(i in 2:length(x))
for(j in 2:length(y))
{
	x.square <- c(x[i-1],x[i],x[i],x[i-1])
	y.square <- c(y[j-1],y[j-1],y[j],y[j])

	z1.square <- c(z1[i-1,j-1], z1[i,j-1], z1[i,j], z1[i-1,j])
	z2.square <- c(z2[i-1,j-1], z2[i,j-1], z2[i,j], z2[i-1,j])
	z3.square <- c(z3[i-1,j-1], z3[i,j-1], z3[i,j], z3[i-1,j])	

	others <- (1:3)[-k]

	myz.square <- get(paste("z",k,".square", sep=""))
	otherz.square1 <- get(paste("z",others[1],".square", sep=""))
	otherz.square2 <- get(paste("z",others[2],".square", sep=""))

	max.z <- max(mean(myz.square), mean(otherz.square1), mean(otherz.square2))
	min.z <- min(mean(myz.square), mean(otherz.square1), mean(otherz.square2))

	if(level=="bottom" & mean(myz.square) == min.z |
		level=="middle" & mean(myz.square) != min.z & mean(myz.square) != max.z |
		level=="top" & mean(myz.square) == max.z)
	{	
		polygon(trans3d(x=x.square,
				y=y.square,
				z=myz.square,pmat=p1), col=mycols[k])
	}
}