function [T,Tbar] = RotationMatrix3D(t);
    R = eye(6,6); R(4,4)=2; R(5,5)=2; R(6,6)=2; % Reuter matrix
    a = [cosd(t), sind(t), 0;sind(t), cosd(t), 0; 0, 0, 1];
    T(1:6,1:6) = 0;
    for i=1:1:3
    for j=1:1:3
    if i==j; alpha = j; else alpha = 9-i-j; end
    for p=1:1:3
    for q=1:1:3
    if p==q beta = p; else beta = 9-p-q; end
    T(alpha,beta) = 0;
    if alpha<=3 & beta<= 3; T(alpha,beta)=a(i,p)*a(i,p); end
    if alpha> 3 & beta<= 3; T(alpha,beta)=a(i,p)*a(j,p); end
    if alpha<=3 & beta>3; T(alpha,beta)=a(i,q)*a(i,p)+a(i,p)*a(i,q);end
    if alpha>3 & beta>3; T(alpha,beta)=a(i,p)*a(j,q)+a(i,q)*a(j,p);end
    end
    end
    end
    end
    Tbar = R*T*R^(-1);
end
