3D coordinates rotation -- new direction for Z axis I need to rotate 3D coordinate system so Z axis points in new direction.
So, I have a direction defined by spherical coordinates ($\theta$, $\phi$), where $\theta$ (in $[0, \pi]$ range) is polar and $\phi$ (in $[0, 2\pi]$ range) is azimuthal angles. I want to transform my 3D Cartesian coordinates so that Z is now pointing in that direction.
Now, I understand that this is not a unique transformation -- I do not care how X and Y axis are going to rotate. I am interested only in having Z axis in the right place. 
Is it possible to get a transformation matrix for this?
 A: The rotation vector can be found be taking the cross product of the original $z$-axis direction with that of the desired direction, while the rotation angle can be found using the dot product. With these two values, you can obtain the required transformation matrix (using this).
A: Here's another approach using a generalized 3x3 rotational matrix.  The theory is simple.  In order to transform points for a coordinate rotate, the rotator matrix has to have all 3 new axes directions defined:
R =  [ newXaxis,  newYaxis,  newZaxis ] (each axis is a direction vector)
         (axes are mutually orthogonal and obey the right hand rule Z <-- X x Y)
Practically, if you can define 2 of the 3 axes, that's enough, since the 3rd axis
is 100% dependent on the choices of the other two. The missing axis can be computed.
R =  [ newXaxis,  newYaxis,  ------- ]    newZAxis <-- newXaxis X newYaxis
R =  [ newXaxis,  --------,  newZaxis ]   newYAxis <-- newZaxis X newXaxis
R =  [ --------,  newYaxis,  newZaxis ]   newXAxis <-- newYaxis X newZaxis
Now, if you only care where one of the 3 axes gets mapped to (as the problem statement specifies), 
R =  [ ---------,  --------,  newZaxis* ] (*see note below)
you can arbitrarily choose a pair of other axes that fulfill the axes rules for the 3x3 rotator matrix.  You use normalized cross-products to obtain them.  The normalized cross-product is the mutual-perpendicular direction finder in 3D, given two arbitrary directions (not coincident nor opposite).
So, for instance:
newXaxis <-- normalize(newZaxis x [ 1  0  0 ])
If by chance newZaxis == [ 1 0 0 ], use as an alternate newXaxis:
newXaxis <-- normalize(newZaxis x [ 0  1  0 ])  
Finally, 
newYaxis <-- newZaxis x newXaxis
Now you can populate your matrix R that does the coordinate rotation. 
=========================================
*Before beginning to assign newZaxis, convert direction angles (phi , theta) to 
equivalent 3D direction vector:
[ cos(phi) * cos(theta) ,   sin(phi) * cos(theta) , sin(theta)  ] 
