Find two angles between two vectors with direction? 
So using the image above I have two vectors, u and v. Vector u is formed using P1 and P2 which are two coordinates in 3D space. Vector v is formed using P1 and P3.
I wish to first rotate vector v around the Z axis until it lines up under vector u which is shown by w in the image. After that then rotate around the x axis, so it is in line with vector u. I wish to be able to calculate those two angles. Hopefully this makes sense.
The goal is to find angle a and b
The other trick is I wish to keep track of direction as well, meaning that the given two angles can be positive or negative.
I do know how to find the angle between two vectors but am unsure how to do so like above.
 A: A counterclockwise (wrt order $(x, y)$) rotation of angle $\theta$ about the $z$-axis is given by the matrix $$M_{\theta,z}=\left(\matrix{\cos \theta& -\sin \theta&0\\\sin \theta& \cos \theta&0\\0&0&1\\}\right)$$
A counterclockwise (wrt order $(y, z)$) rotation of angle $\varphi$ about the $x$-axis is given by the matrix $$M_{\varphi,x}=\left(\matrix{1&0&0\\0&\cos \varphi& -\sin \varphi\\0&\sin \varphi& \cos \varphi\\}\right)$$
You have to find $\theta$ and $\varphi$ such that $M_{\varphi,x}\cdot M_{\theta,z}\cdot v=u$. If you know the coordinates of $u$ and $v$ this will give you a system of equations to solve.
A: You’re essentially looking for the difference in the spherical coordinates of the two vectors, so you can use the standard formulas for conversion from Cartesian coordinates: $$r=\sqrt{x^2+y^2+z^2} \\ \theta = \arccos{\frac zr} \\ \phi = \arctan{\frac yx}.$$ Remember that $\theta$ is measured from the positive $z$-axis, and you’ll need to adjust $\Delta\theta$ to be in the range $(-\pi,\pi]$ if you want the use both clockwise and counterclockwise rotations about the $z$-axis for the first rotation.  
If all you really need is the two rotations, though, they can be constructed without explicitly computing either angle. The basic idea is to build sets of orthonormal bases that represent the “before” and “after” states. For the first rotation, project $\mathbf u$ and $\mathbf v$ onto the $x$-$y$ plane and normalize to get $\mathbf u' = (u_x^2+u_y^2)^{-1/2}(u_x,u_y,0)^T$ and $\mathbf v' = (v_x^2+v_y^2)^{-1/2}(v_x,v_y,0)^T$, respectively. The required rotation matrix is then $$R_1 = \begin{bmatrix}\mathbf u'&\mathbf e_z\times\mathbf u'&\mathbf e_z\end{bmatrix} \begin{bmatrix}\mathbf v'&\mathbf e_z\times\mathbf v'&\mathbf e_z\end{bmatrix}^{-1} = \begin{bmatrix}\mathbf u'&\mathbf e_z\times\mathbf u'&\mathbf e_z\end{bmatrix} \begin{bmatrix}\mathbf v'^T \\ (\mathbf e_z\times\mathbf v')^T \\ \mathbf e_z^T\end{bmatrix},$$ with $\mathbf e_z = (0,0,1)^T$, the unit $z$-vector. (Note that computing the cross product with $\mathbf e_z$ are a matter of swapping coordinates and negating one of them.)  
The second rotation is constructed in a similar manner, except that you’re now rotating $R_1\mathbf v$ onto $\mathbf u$, so you’ll need to normalize both of those vectors and use $\mathbf w = R_1\mathbf v\times\mathbf u$, normalized as the rotation axis.  
If needed, the rotation angle $\theta$ can be extracted from each of the resulting rotation matrices using the identity $\operatorname{tr}R = 1+2\cos\theta$.
