Extracting a rotation axis from a rotation matrix

I am working with 3 x 3 rotation matrices, R(ij), of cubic bicrystal misorientations, so there are up to 24 symmetry related descriptions (i.e. different permutations of the rotation matrix) for any starting misorientation. Usually I extract the axis/angle description from each one by standard methods - acos (theta) = acos(((traceR)-1)/2), and the axis as [R32-R23, R13-R31, R21-R12]. This works fine, except when theta = 180 degrees, for which this formula gives the axis as [0 0 0], as we know. So to get the axes for the 180 case, it seems to me I can do the following:

1) Use [V,D] = eig® to get a matrix D whose columns are the eigenvalues of R. One of these will be equal to +1, the other two will be -1.

2) Identify which column of D contains the +1, and then take the equivalent column of V as the eigenvector, i.e. the rotation axis. So far so good.

I am trying to use a for-(else)-end loop to do this, by using a counter i to examine each column of D in turn, and then summing the three elements of that column. If the sum comes to +1, I want to use that corresponding i to obtain the eigenvector as:

[V(1,i), V(2,1), V(3,1)], and then come out of the loop. I am relatively new to MATLAB, and I cannot get a loop to work. Can anyone help point me to a working code to do this loop?

Thanks.

Answer 1

[Andrew Newell]

tol = 100*eps;
i = find(abs(diag(D)-1)<tol);
rotAxis = V(:,i);

(Edited to take rounding error into account. Increase tol if necessary.)

Answer 2

[William]

Andrew, I tried this, but it comes back with the response i = Empty matrix 0-by-1 for some of the 24 symmetry operations, and then these return "Empty matrix 3-by-0" for the rotation axis, and produces no axis.

This is the loop for getting the degenerate descriptions:

for i=1:1:24

    RotMat2 = RotMat*e(:, :, i);

    [V,D] = eig(RotMat2);

    RotAng2 = acosd((trace(RotMat2)-1)/2);

    j = find(diag(D)==1)

    RotAxis2 = V(:,j)

"Then come some printing commands, to print RotAng2, RotAxis2, and RotMat2"

end

Where RotMat is the starting rotation matrix. I am sure I am missing the obvious, but I can't get this to work. Hoping a more experienced user can :-)

[William]

I should have added that the matrices e are the 24 symmetry operations of the cubic point group.

[Walter Roberson]

Take one of the cases where you are expecting 1, and subtract 1 from it and look at the result. I suspect that what you will see a small residue, indicating that the value is not exactly 1 due to round-off error. It is seldom a good idea to use "==" for a floating point calculation.

[Andrew Newell]

Right. That was careless of me. I have corrected my answer.

Answer 3

[William]

Gentlemen, thank you both for the help. So far, this seems to work.

Much appreciated, from a newbie :-)