tPower <- function(alpha,n,sigma,mu,muZero) {
  c <- qt(1-alpha,n-1)
  ncp <- (sqrt(n)/sigma) * (mu - muZero)
  return(1-pt(c,n-1,ncp))
}


alpha <- 0.05
n <- 5
muZero <- 0
sigma <- 1
mu <- seq(muZero,muZero+3,0.01)

#plot(mu,tPower(alpha,n,sigma,mu,muZero),type="l",xlab="mu",ylab="power",lty=1, col=1)
#par(new=TRUE)
#sigma<- 0.5
#plot(mu,tPower(alpha,n,sigma,mu,muZero),type="l",xlab="mu",ylab="power",lty=2, col=2)
#legend(2,0.2,c("sigma=1","sigma=0.5"),lty=c(1,2), col=c(1,2))


plot(mu,tPower(alpha,n,sigma,mu,muZero),type="l",xlab="mu",ylab="power",lty=1, col=1)
par(new=TRUE)
n<- 10
plot(mu,tPower(alpha,n,sigma,mu,muZero),type="l",xlab="mu",ylab="power",lty=2, col=2)
par(new=TRUE)
n<- 30
plot(mu,tPower(alpha,n,sigma,mu,muZero),type="l",xlab="mu",ylab="power",lty=3, col=3)
legend(2,0.4,c("n=5","n=10","n=30"),lty=c(1,2,3), col=c(1,2,3))


abline(0.05,0)
abline(0.8,0)