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) ]