#########################################################################
## this code compute the power of an F test comparing
## the equality of several means, in 1 way ANOVA setting.
## Input:  
##   meanvec: the vector of several means, measured in the unit
##            of the common sigma.  i.e. mean[1]=2.5 really says
##            the mean[1]=2.5sigma, mean[2]=-1.2 actually says
##            mean[2]=-1.2sigma etc. (sigma can be unknown).
##            In most cases they are not all equal (alternative). 
##        
##   nvec:    the sample sizes of several groups, their means
##            are been tested for equality. Do not have to be all equal.
##   alpha:   Optional level of significance. Defalt to 0.05.
## Output:
##   power at the given alternative.
#########################################################################
power.f <- function(meanvec, nvec, alpha=0.05)
{
      if( any(is.na(meanvec)) || any(is.na(nvec)) ) 
         stop("there are NA in the input")
      if(length(meanvec) != length(nvec))
         stop("meanvec and nvec must be of same length")
      if( any( nvec<0 ) ) stop("negative sample size?")
      if( (alpha<=0) || (alpha>=1) ) stop("alpha must inside (0,1)")

      v <- length(nvec) - 1
      vv <- sum(nvec) - v - 1 

      meanbar <- sum(meanvec * nvec)/sum(nvec) 
      ncp <- sum(nvec * (meanvec - meanbar)^2) 

      1 - pf(qf(1-alpha, df1=v, df2=vv), df1=v, df2=vv, ncp=ncp)
}