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

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

xbar <- numeric()
cn   <- numeric()

xbar[1] <- mean(dat[1,])
cn[1]   <- xbar[1]

for(i in 2:20){
   xbar[i] <- mean(dat[i,])
   cn[i] <- cn[i-1]+xbar[i]
}

ii <- seq(1,20)

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

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

h=0.75
mask1 <- rep(0,20)

plot(ii[1:20],cn[1:20],type="o",lty=1,pch=16,xlab="n",
     ylab=expression(C[n]),mgp=c(2,1,0),xlim=c(0,21), 
     ylim=c(-3,4),cex=0.8)

for (i in 11:20){
mask1[i] <- cn[11]-h+(i-11)*0.1
}

lines(ii[11:20],mask1[11:20],type="l",lty=1,cex=0.8)
title(xlab="(a)",cex=0.9)

mask2 = rep(0,20)

plot(ii[1:20],cn[1:20],type="o",lty=1,pch=16,xlab="n",
     ylab=expression(C[n]),mgp=c(2,1,0),xlim=c(0,21), 
     ylim=c(-3,4),cex=0.8)

for (i in 11:20){
mask1[i] = cn[20]-h+(i-20)*0.1
mask2[i] = cn[20]+h+(20-i)*0.1
}

lines(ii[11:20],mask1[11:20],type="l",lty=1,cex=2.5)
lines(ii[11:20],mask2[11:20],type="l",lty=2,cex=2.5)
lines(rep(20,2),c(cn[20]-h,cn[20]+h),lty=3,cex=2.5)
lines(rep(11,2),c(mask1[11],mask2[11]),lty=3,cex=2.5)
cord.x <- c(11,11,20,20)
cord.y <- c(mask1[11],mask2[11],cn[20]+h,cn[20]-h)
polygon(cord.x,cord.y,density=7,angle=135,border=NA,lty=1,cex=0.6)

title(xlab="(b)",cex=0.9)

graphics.off()