##### R-code for making Figure 4.1. It reads data from
##### "example41.dat" generated by the R-code "example41.r"

dat <- matrix(scan("example41.dat"),ncol=5,byrow=T)

xbar <- numeric()
rbar <- numeric()

for(i in 1:20){
   xbar[i] <- mean(dat[i,])
   rbar[i] <- max(dat[i,])-min(dat[i,])
}

i <- seq(1,20)

postscript("fig41.ps",width=7.2,height=3.6,horizontal=F)

par(mfrow=c(1,2), mar=c(4,4,1,2))

U=0.1005963+3/(2.326*sqrt(5))*2.025616
C=0.1005963
L=0.1005963-3/(2.326*sqrt(5))*2.025616

plot(i[1:10],xbar[1:10],type="o",lty=1,pch=16,xlab="i",
     ylab=expression(bar(X)[i]),mgp=c(2,1,0),xlim=c(0,21), 
     ylim=c(-1.7,1.7), cex=0.8)
lines(i[11:20],xbar[11:20],type="o",lty=3,pch=1,cex=0.8)
lines(i,rep(U,length(i)),lty=2,cex=0.8)
lines(i,rep(C,length(i)),lty=2,cex=0.8)
lines(i,rep(L,length(i)),lty=2,cex=0.8)
axis(4,at=c(L,C,U),labels=c("L","C","U"),cex=0.8,las=1)
title(xlab="(a)",cex=0.9)

U=3/sqrt(5)
C=0
L=-3/sqrt(5)

plot(i[1:10],xbar[1:10],type="o",lty=1,pch=16,xlab="i",
     ylab=expression(bar(X)[i]),mgp=c(2,1,0),xlim=c(0,21), 
     ylim=c(-1.7,1.7), cex=0.8)
lines(i[11:20],xbar[11:20],type="o",lty=3,pch=1,cex=0.8)
lines(i,rep(U,length(i)),lty=2,cex=0.8)
lines(i,rep(C,length(i)),lty=2,cex=0.8)
lines(i,rep(L,length(i)),lty=2,cex=0.8)
axis(4,at=c(L,C,U),labels=c("L","C","U"),cex=0.8,las=1)
title(xlab="(b)",cex=0.9)

graphics.off()