function newP=rot3d(P,A,D,phi)
%ROT3D rotation about an arbitrary axis in 3D
%   Q = rot3d(P,A,D,phi) returns the coordinates of Q obtained by
%   rotation of point P about the axis through point A, parallel
%   to direction D, by angle phi (right hand rule about D)
%
%   3D cartesian space only!
%   P,A,D should be 3-component vectors (assumed to be cartesian comps.)
%   phi is in radians
%
%   examples:
%    Q=rot3d([1,0,0],[0,0,0],[0,0,1],pi/4)
%    Q=rot3d([1,0,0],[0,0,0],[0,1,0],pi/2)
%
%   Code is NOT FOOL-PROOF! 
%   P,A,D should be 1-by-3 row vectors of 3-by-1 column vectors but not a mix.
% 
%   See also DOT, CROSS.
%
%   Fabian Waleffe at UW for Math 321, Feb 20, 2009


D=D/norm(D);       % normalize D to unit vector
AP=P-A;            % vector AP
AP1=dot(AP,D)*D;   % AP parallel to D
AP2=AP-AP1;        % AP perpendicular to D

newP= A + AP1 + cos(phi)*AP2 + sin(phi)*cross(D,AP2);


