I have 3 parameters: A, B and C with 20000 data each of them. When I put them together in an 3D graph, they form a dots cloud being A in X axis, B in Y axis and C in Z axis. I want to fix this dots cloud in an 90% covariance error ellipsoid and determine the lenghts, i.e. width, height and strecth, in > 1000 samples individually. Width will be in X axis, height in Y axis and strecth in Z axis). So, for each sample:

1) I calculate the covariation matrix with var()

covMat <- var(cbind(A,B,C)) 
covMat
           A          B          C
A 0.08269010 0.05330448 0.05329601
B 0.05330448 0.33033311 0.14017367
C 0.05329601 0.14017367 0.18732865

2) I calculate the eigenvalues with eigen()

evals<-eigen(covMat)$values
evals
[1] 0.43162433 0.10861830 0.06010923

#I obtain 3 eigenvalues ordered from largest-to-smallest

3) I determine the lenghts of the 90% covariance error ellipsoid

ell.len <- 2*sqrt(6.251*evals)
ell.len
[1] 2.971533 1.666728 1.325708

The tricky thing: Since eigenvalues are always in a decreasing order, how can I know which number of the ell.len correspond to each lenght (width, height and strecth) since sometimes height can be larger than width and strecth, and other combinations. I want to know this because I would like to determine numerically when the ellipsoids are prolated (Height > Width > Strecth), oblate (Height > Width > Height > Strecth) and other possible combinations etc in each sample (n = 1000).

Thank you very much in advance and I hope my question is clear.