rgl: ellipse3d with axes

Last week I asked about data ellipses with rgl:::ellipse3d() with lines 
showing the principal axes.
(The goal is a visual demonstration of PCA as a rotation of variable 
space to component space.)
I was trying, unsuccessfully, to use princomp() to generate the PCA axes 
and plot them using
segments3d:

PC <- princomp(trees)
sdev <- PC$sdev         # component standard deviations
sd <- sqrt(diag(cov))   # variable standard deviations

# vectors in variable space of principal components
vec <- matrix(mu,3,3, byrow=TRUE) + diag(sd) %*% PC$loadings

for (j in 1:3) {
 mat <- rbind(mu, vec[j,])
 segments3d(mat, col="red")
}

However, it occurred to me that these axes are just the orthogonal axes 
of the unit sphere
that is transformed (using chol()) in ellipse3d, so plotting the axes 
transformed in the
same way would give me what I want.

Looking at the result returned by ellipse3d, I see a normals component, 
but I'm not sure if this
represents what I want, or, if it is, how to use it to draw the ellipse 
major axes in the plot.

 > e1 <-ellipse3d(cov, centre=mu, level=0.68)
 > str(e1)
List of 6
 $ vb           : num [1:4, 1:386] 4.95 2.64 2.03 1.00 6.74 ...
 $ ib           : num [1:4, 1:384] 1 195 99 196 51 197 99 195 27 198 ...
 $ primitivetype: chr "quad"
 $ homogeneous  : logi TRUE
 $ material     : list()
 $ normals      : num [1:4, 1:386]  0.290 -0.902 -0.320  1.000  0.635 ...
 - attr(*, "class")= chr "qmesh3d"

-Michael

Michael Friendly     Email: friendly AT yorku DOT ca 
Professor, Psychology Dept.
York University      Voice: 416 736-5115 x66249 Fax: 416 736-5814
4700 Keele Street    http://www.math.yorku.ca/SCS/friendly.html
Toronto, ONT  M3J 1P3 CANADA