plot.power = function(my0=0,my1=0,sd=1,n=30,alpha=.05)
{
  # Calculates and plots power of one sample z-test of a sample mean my1 against
  # a population mean my0 (H0: my0 = my1, H1: my0 <> my1).
  #
  # Defaults: µ0 = 0, µ1 = 0, sigma(population) = 1, n = 30,
  #           alpha(twosided) = .05
  
  gamma = (my1-my0)/sd
  delta = gamma*sqrt(n)
  se.mean = sd/sqrt(n)

  z0L = my0-qnorm(1-alpha/2)*se.mean
  z0U = my0+qnorm(1-alpha/2)*se.mean
  z0min = my0-3.5*se.mean
  z0max = my0+3.5*se.mean
  z1min = my1-3.5*se.mean
  z1max = my1+3.5*se.mean
  if (my1 > my0) mean.crit=z0U
  else mean.crit=z0L

  if (my1 != my0)
  {
    cat('\nOne-sample z-test power calculation\n\n')
    cat('          n =',n,'\n')
    cat('      sigma =',round(sd,3),'\n')
    cat('   SE(mean) =',round(se.mean,3),'\n')
    cat('         µ0 =',round(my0,3),'\n')
    cat('         µ1 =',round(my1,3),'\n')
    cat('  mean.crit =',round(mean.crit,3),'\n\n')
    cat('effect size =',round(gamma,3),'\n')
    cat('      delta =',round(delta,3),'\n')
    cat('  sig.level =',round(alpha,3),'\n')
    cat('      power =',round(1-pnorm(qnorm(1-alpha/2,mean=0,sd=se.mean),
                       mean=abs(my1-my0),sd=se.mean),3),'\n')
    cat('alternative = two-sided\n\n')
    if (n < 30) cat('Warning: N too small for z test!\n\n')
  }

  # Normal curve H0:

  curve(dnorm(x,mean=my0,sd=se.mean),from=z0min,to=z0max,
              xlim=range(z0min,z0max,z1min,z1max),ylab='density',col="blue",
              xlab=paste('alpha = ',round(alpha,3),' (two-sided)',
                         ', beta = ',round(pnorm(qnorm(1-alpha/2,mean=0,
                                     sd=se.mean),mean=abs(my1-my0),
                                     sd=se.mean),3),
                         ', power = ',round(1-pnorm(qnorm(1-alpha/2,mean=0,
                                     sd=se.mean),mean=abs(my1-my0),
                                     sd=se.mean),3),sep=""),
              main=paste('Power of z-Test of the Mean of a Single Population\n',
                         'n = ',n,', µ0 = ',my0,', µ1 = ',my1,', sigma = ',sd,
                         sep=""))
  # Acceptance region H0:
  x=seq(z0L,z0U,min(.001,1/n))
  polygon(c(z0L,x,z0U),c(0,dnorm(x,mean=my0,sd=se.mean),0),col="lightyellow")

  # Normal curve H1:
  curve(dnorm(x,mean=my1,sd=se.mean),from=z1min,to=z1max,add=T,col="red")

  text(my0,0,'µ0',pos=1,offset=.15,cex=.8)
  text(my1,0,'µ1',pos=1,offset=.15,cex=.8)

  if (my1 > my0)
  {
    text(z0U,0,round(z0U,2),pos=1,offset=.15,cex=.8)

    # Acceptance region H1:
    x=seq(z0U,z1max,min(.001,1/n))
    polygon(c(z0U,x,z1max),c(0,dnorm(x,mean=my1,sd=se.mean),0),
            col="lightyellow")

    # Rejection region H0:
    x=seq(z0U,z0max,min(.001,1/n))
    polygon(c(z0U,x,z0max),c(0,dnorm(x,mean=my0,sd=se.mean),0),density=20,
            col="blue")
    x=seq(z0min,z0L,min(.001,1/n))
    polygon(c(z0min,x,z0L),c(0,dnorm(x,mean=my0,sd=se.mean),0),density=20,
            col="blue")

    # Rejection region H1:
    if (z1min < z0U)
    {
      x=seq(z1min,z0U,min(.001,1/n))
      polygon(c(z1min,x,z0U),c(0,dnorm(x,mean=my1,sd=se.mean),0),density=20,
              angle=135,col="red")
    }
  }
  else if (my1 < my0)
  {
    text(z0L,0,round(z0L,2),pos=1,offset=.15,cex=.8)

    # Acceptance region H1:
    x=seq(z1min,z0L,min(.001,1/n))
    polygon(c(z1min,x,z0L),c(0,dnorm(x,mean=my1,sd=se.mean),0),
            col="lightyellow")

    # Rejection region H0:
    x=seq(z0U,z0max,min(.001,1/n))
    polygon(c(z0U,x,z0max),c(0,dnorm(x,mean=my0,sd=se.mean),0),density=20,
            col="blue")
    x=seq(z0min,z0L,min(.001,1/n))
    polygon(c(z0min,x,z0L),c(0,dnorm(x,mean=my0,sd=se.mean),0),density=20,
            col="blue")

    # Rejection region H1:
    if (z0L < z1max)
    {
      x=seq(z0L,z1max,min(.001,1/n))
      polygon(c(z0L,x,z1max),c(0,dnorm(x,mean=my1,sd=se.mean),0),density=20,
              angle=135,col="red")
    }
  }
  else cat('Error: µ1 must differ from µ0!\n')
}

# Example my0=68, sigma=3.1, my1=69, n=100, alpha(twosided)=.05
# (de-comment the following line):
# plot.power(my0=68,my1=69,sd=3.1,n=100)