t1 = 30.0 * pi / 180;
t2 = 60.0 * pi / 180;

c1 = cos(t1); s1 = sin(t1);
c2 = cos(t2); s2 = sin(t2);
T1 = [1 0 0; 0 c1 -s1;0 s1 c1]; 

T2 = [c2 0  -s2;0   1   0; s2 0 c2]; 
T12 = T1 * T2
[V12,D12] = eig(T12);
u12 = [];
for i = 1:3
    if (abs(D12(i,i) - 1.0) < 1e-3)
        u12 = V12(:,i);
    end
end

if (length(u12) == 0)
    fprintf(1, 'error cannot find eigenvalue of value 1\n');
end

fprintf(1, 'Method A (eigenvalue): Obtaining results with eigenvalue analysis\n');
fprintf(1, 'Rotation vector for T1T2 is[%g %g %g]', u12(1), u12(2), u12(3));
traceV12 = trace(T12);
rotationAngle12A = acos((traceV12 - 1) / 2.0) * 180 / pi

fprintf(1, 'Method B (skew symmetric part of T12): Obtaining results with skew symmetric part of rotation\n');
c = (traceV12 - 1) / 2.0;
s = sqrt(1 - c * c);
axis12(1) = 1 / (2.0 * s) * (T12(3, 2) - T12(2, 3));
axis12(2) = 1 / (2.0 * s) * (T12(1, 3) - T12(3, 1));
axis12(3) = 1 / (2.0 * s) * (T12(2, 1) - T12(1, 2));
rotationAngle12B = rotationAngle12A
fprintf(1, 'Method B (skew symmetric part of T12): Rotation vector for T1T2 is[%g %g %g]\n', axis12(1), axis12(2), axis12(3));
fprintf(1, 'As we can see rotation vector from eigenvector approach and its roation angle are not consistent. Either its direction must change or ange changed to its opposite value\n');

T21 = T2 * T1
fprintf(1, 'Method B (skew symmetric part of T21): Obtaining results with skew symmetric part of rotation\n');
traceV21 = traceV12;
fprintf(1, 'Why traceV21 = traceV12?\n');
fprintf(1, 'So c and s do not change\n');
axis21(1) = 1 / (2.0 * s) * (T21(3, 2) - T21(2, 3));
axis21(2) = 1 / (2.0 * s) * (T21(1, 3) - T21(3, 1));
axis21(3) = 1 / (2.0 * s) * (T21(2, 1) - T21(1, 2));
fprintf(1, 'and the rotation angles for T12 and T21 (Why?)\n');
rotationAngle21B = rotationAngle12B
fprintf(1, 'Method B (skew symmetric part of T21): Rotation vector for T2T1 is[%g %g %g]\n', axis21(1), axis21(2), axis21(3));
fprintf(1, 'Comment whether T12 and T21 represent the same rotation.\n');

fprintf(1, 'Finally we want to reconstruct the rotation tensor T21 having its axis and angle or rotation\n');
skewTensorOfAxis = [0 -axis21(3) axis21(2);axis21(3) 0 -axis21(1);-axis21(2) axis21(1) 0];
theta = rotationAngle12A / 180 * pi;
I = eye(3);
T21Reconstructed = cos(theta) * I + (1 - cos(theta)) * axis21' * axis21 + sin(theta) * skewTensorOfAxis
T21
