shade.dnorm=function(m,s,a,b,Col="red"){
#add to plot of normal density, shading area between a and b
x = seq(a,b,length=100)
y=dnorm(x,m,s)
n = length(x)
X = c(x[1],x,x[n])
Y =c(0,y,0)
polygon(X,Y,col=Col)
}

PlotZTest = function(mu,sigma,n,alpha=.05,C0="red"){
se = sigma/sqrt(n)
x0 = mu - 4*se
x1 = mu + 4*se
X = seq(x0,x1,length=400)
Y0 = dnorm(X,mu,se)
plot(X,Y0,type="l",lwd=3, ylab="Density")
title(paste("alpha = ",alpha,sep="")) 
abline(h=0)
shade.dnorm(mu,se,x0,qnorm(alpha/2,mu,se),Col=C0)
shade.dnorm(mu,se,qnorm(1-alpha/2,mu,se),x1,Col=C0)
}

PlotZPower = function(mu0,sigma,n,mu1,alpha=.05,C0="red",C1="green"){
if(mu1 < mu0) mu1 = mu0+abs(mu0-mu1)
se = sigma/sqrt(n)
x0 = min(mu0,mu1)-4*se
x1 = max(mu0,mu1) + 4*se
X = seq(x0,x1,length=400)
Y0 = dnorm(X,mu0,se)
Y1 = dnorm(X,mu1,se)
plot(X,Y0,type="l",lwd=2)
title(paste("Power = ",
    round(pnorm(qnorm(1-alpha/2,mu0,se),mu1,se,lower.tail=FALSE),3),sep="")) 
lines(X,Y1,lty=2,lwd=2)
abline(h=0)
shade.dnorm(mu0,se,x0,qnorm(alpha/2,mu0,se),Col=C0)
shade.dnorm(mu1,se,qnorm(1-alpha/2,mu0,se),x1,Col=C1)
shade.dnorm(mu0,se,qnorm(1-alpha/2,mu0,se),x1,Col=C0)
}