function rot = getRotationMatrix3d(axis, a)
% obtain rotation matrix based on angle and coordinate axis OR a vector
%
% Syntax: rot = getRotationMatrix3d(axis, angle)
%   or    rot = getRotationMatrix3d('vector', v)
%
% The axis is 'x', 'y' or 'z' and the angle must be provided in radian.
%   With the case 'vector', the the rotation matrix is calculated so that
%   the normalized vector v matches the unit vector if the z-axis.
%
% See also: getRotationMatrix2d, getRotationMatrix4d, rotatePolytope

% The elk-library: convex geometry applied to crystallization modeling.
%   Copyright (C) 2012 Alexander Reinhold
%
% This program is free software: you can redistribute it and/or modify it 
%   under the terms of the GNU General Public License as published by the 
%   Free Software Foundation, either version 3 of the License, or (at your 
%   option) any later version.
%
% This program is distributed in the hope that it will be useful, but 
%   WITHOUT ANY WARRANTY; without even the implied warranty of 
%   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU 
%   General Public License for more details.
%
% You should have received a copy of the GNU General Public License along 
%   with this program.  If not, see <http://www.gnu.org/licenses/>.

switch lower(axis)
    case 'x'
        rot = [1       0       0;
               0       cos(a) -sin(a);
               0       sin(a)  cos(a)];
    case 'y'
        rot = [cos(a)  0       sin(a);
               0       1       0;
              -sin(a)  0       cos(a)];
    case 'z'
        rot = [cos(a) -sin(a)  0;
               sin(a)  cos(a)  0;
               0       0       1];
    case 'vector'
        % correct v to be a column vector
        if size(a, 1) == 1
            a = a';
        end
        
        % rotate vector so that it lies in xz-plane (eliminates y-direction)
        if a(2) ~=0
            anglez = acos(dot(a(1:2) / norm(a(1:2)), [1,0])) * ...
                -sign(a(2));
            matrix_z = getRotationMatrix3d('z', anglez);
            a = matrix_z*a;
        else
            matrix_z = eye(3);
        end
        
        % rotate vector so that it lies in yz-plane (eliminates z-direction)
        % normal is then colinear with z-axis
        if a(1) ~= 0
            angley = acos(dot(a([1,3]) / norm(a([1,3])),[0,1])) * ...
                -sign(a(1));
            matrix_y = getRotationMatrix3d('y', angley);
        else
            matrix_y = eye(3);
        end
        
        % assign rotation matrix
        rot = matrix_y*matrix_z;
        
    otherwise
        error(['Axis ''' axis ''' unknown, use ''x'', ''y'', ''z'' or ''vector''.'])
end
