f.gen.sphere<-function(n=100,p=5) {
 x<-matrix(rnorm(n*p),ncol=p)
 xnew<-t(apply(x,1,norm.vec))
 xnew
 }

norm.vec<-function(x) {
 x<-x/norm(x)
 x
 }

norm<-function(x) { sqrt(sum(x^2))}

# Read in the data
d.femkite<-read.table("femkite.asc")
fm<-apply(d.femkite,2,mean) # Calculate the mean
vc<-var(d.femkite)  # Calculate the variance-covariance
# Compute var-cov^(1/2)
evc<-eigen(vc)
vc2<-(evc$vectors)%*%diag(sqrt(evc$values))%*%t(evc$vectors)
# Compute critical value
critval<-sqrt((2*44)/(45*43)*qf(0.95,2,43))
y<-f.gen.sphere(100,2) # Sample from uniform on a sphere
y<-y%*%vc2 # Transform to an ellipse
y<-y*critval # Scale proportional to 95% critical value
y<-y+matrix(rep(fm,100),ncol=2,byrow=T) # Shift to center at mean

# Plot data and confidence ellipse
plot(d.femkite[,1],d.femkite[,2],xlab="Tail",ylab="Wing",pch="x")
points(y[,1],y[,2],pch=16)