par(mfrow=c(2,2))

M <- 2
ngrd <- 30
l <- (2*ngrd+1)
lplt <- (M+1)*l^M

plt <- scan("mono2.plt")
zinc <- plt[seq(from=M+(M*(M+1))/2+1,to=M+(M*(M+1))/2+M)]
plt <- plt[seq(from=M+(M*(M+1))/2+M+1,to=M+(M*(M+1))/2+M+lplt)]

h <- matrix(plt[seq(from=M+1,to=lplt,by=M+1)],nrow=l,ncol=l)
x <- plt[seq(from=1,to=l*(M+1),by=M+1)]     
y <- plt[seq(from=2,to=lplt,by=l*(M+1))]

persp(x,y,h,xlab="Clearance",ylab="Volume",zlab="h")
title("(a)")

Eh <- -sort(-cumsum(sort(zinc[1]*zinc[2]*h)))
Oh <- order(-h)
qh <- c(h[Oh[length(Eh[Eh>.95])]], h[Oh[length(Eh[Eh>.90])]],
        h[Oh[length(Eh[Eh>.75])]], h[Oh[length(Eh[Eh>.50])]],
        h[Oh[length(Eh[Eh>.25])]], h[Oh[length(Eh[Eh>.10])]])

contour(x,y,h,labex=0,xlab="Clearance",ylab="Volume",levels=qh)
title("(b)")

plt <- scan("mono2.p02")
x <- plt[seq(from=3+1,to=3+2*l,by=2)]
h <- plt[seq(from=3+2,to=3+2*l,by=2)]

plot(x,h,type='l',xlab="Volume",ylab="h")
title("(c)")

plt <- scan("mono2.p01")
x <- plt[seq(from=3+1,to=3+2*l,by=2)]
h <- plt[seq(from=3+2,to=3+2*l,by=2)]

plot(x,h,type='l',xlab="Clearance",ylab="h")
title("(c)")
