The following code defines an S-PLUS function to calculate
the power of the F test, using the method of Pearson and
Hartley (1951). The parameters (rescaled noncentrality parameter
phi, size of test alpha, degrees of freedom df1 in the numerator,
df2 in the denominator) are the same as defined by Pearson and
Hartley in constructing their charts.

To use: first run the S-PLUS code (bottom of this page) in
S-PLUS. For example, you could copy the entire sequence into
a file, and then run a 'source' command from within S-PLUS.
This creates a function 'pearsonhartley'. To run the function,
simply substitute the desired values of the four parameters.
For example, to reproduce the power calculations of the example
on page 133, here is the corresponding S-PLUS code (after
creating the function): 

> pearsonhartley(2,.05,4,10)
[1] 0.8212127
> pearsonhartley(2.5,.05,4,10)
[1] 0.957165

Thus the powers in the two examples (to two decimal places) are
respectively .82 and .96, as stated in the text.

Warning: The code is not numerically stable for large values of 
phi, and may therefore fail  in this case. However, in that
case, the power is effectively 1, as can be checked by trying 
smaller values of phi. Thus, for example, when we increase the
value of phi in the above example, we get

> pearsonhartley(3,.05,4,10)
[1] 0.9940888
> pearsonhartley(4,.05,4,10)
[1] 0.9999798
> pearsonhartley(5,.05,4,10)
[1] 1
> pearsonhartley(10,.05,4,10)
[1] 1
> pearsonhartley(20,.05,4,10)
[1] NA

For phi=20, the program fails, but we can easily see that this
corresponds to power 1 because of the results for smaller phi.

Also note that on the way, the program uses a built-in S-PLUS function

q<-qf(1-alpha,df1,df2)

which evaluates the quantile of an F distribution, with df1 and df2 degrees
of freedom, corresponding to upper tail probability alpha. The inverse
function

p<-1-pf(q,df1,df2)

evaluates the upper tail probability (p-value) associated with
an F statistic q.

The S-PLUS code:

pearsonhartley<-function(phi,alpha,df1,df2){
b<-0.5*phi*phi*(1+df1)
q<-qf(1-alpha,df1,df2)
y<-1-1/(1+df1*q/df2)
p<-pbeta(y,df1/2,df2/2)
p1<-1
jmax<-1000
for(j in 1:jmax){
p1<-p1*b/j
p<-p+pbeta(y,j+df1/2,df2/2)*p1
}
1-p*exp(-b)
}

Reference:

Pearson, E.S. and Hartley, H.O. (1951),
Charts of the power function of the analysis of variance tests,
derived from the non-central F distribution.
Biometrika 38 , 112-130.