##########################
#The Power Function of the two-sided t-test
###############################
par(mfrow=c(1,1))
alpha=0.05
d=seq(0,3.6,0.1)
nvalues=c(3,4,5,9,15,20,30,40,50,75,100)
n=nvalues[1]
t=qt(1-alpha/2,n-1,0,lower.tail=T)
beta=1-(pt(t,n-1,sqrt(n)*d,lower.tail=T)-pt(-t,n-1,sqrt(n)*d,lower.tail=T))
plot(d,beta,type="l",lty="dashed",
ylim=c(0,1),xlim=c(0,3.6),col="blue",lwd=2,
main=expression("The power function of the two sided t-test at "*alpha==0.05),
ylab=expression(beta(psi)),xlab=expression(abs(psi)/sqrt(n)))
axis(1,at=seq(0,3.6,0.1),labels=T)
axis(2,at=seq(0,1,0.1),labels=T)
for (n in nvalues[-1]){
t=qt(1-alpha/2,n-1,0,lower.tail=T)
beta=1-(pt(t,n-1,sqrt(n)*d,lower.tail=T)-pt(-t,n-1,sqrt(n)*d,lower.tail=T))
lines(d,beta,lty="dashed",col="blue",lwd=2)}
# These following x and y values to put degrees of freedom at are determined by using the locator() function 
# interactively on the graphics screen 
pts<-locator(n=length(nvalues))
text(pts$x,pts$y,nvalues-1,col="red")
text(2.8,0.9,"The numbers in red are degrees of freedom.",col="red")