power.ibd <- function(al,p,g0,n) 
{
   add.var <- 2*p*(1-p)*al^2   
   al.hat <- add.var/2
   Q2 <- (g0+2*p*al)^2 + al.hat 
   return(sqrt(n/2)*al.hat/Q2)
}

g0 <- 0.05
al <- seq(0,3,by=0.2)
n <- 100
p <- 0.1

plot(al,power.ibd(al,p,g0,n),type="l",xlab="alpha",ylab="ncp") 
lines(al,power.ibd(al,p,g0,30),col=2) 
legend(2,1,legend=c("n=100","n=30"),lty=c(1,1),col=1:2)