Proving the Eulerian angles rotate surjectively and how to create new sets of three rotation angles I am wondering how one would create a new set of angles such that when you rotate a coordinate system with each of them about a chosen axis, you get any possible rotated coordinate system.
For example, the Eulerian angles ϕ, θ, and Ψ can create any orientation of new orthogonal axes.
The set of coordinates is rotated counterclockwise about the z-axis an angle ϕ with the transformation matrix,
$$
\begin{pmatrix}
cosϕ & sinϕ & 0 \\
-sinϕ & cosϕ & 0 \\
0 & 0 & 1 \\
\end{pmatrix}
$$
The next angle then rotates the new set of coordinates counterclockwise about the new x-axis with the transformation,
$$
\begin{pmatrix}
1 & 0 & 0 \\
0 & cosθ & sinθ \\
0 & -sinθ & cosθ \\
\end{pmatrix}
$$
The final angle then rotates counterclockwise about the new z-axis with the transformation,
$$
\begin{pmatrix}
cosΨ & sinΨ & 0 \\
-sinΨ & cosΨ & 0 \\
0 & 0 & 1 \\
\end{pmatrix}
$$
The question is, how do I prove that a set of rotation matricies will be able to give all possible orientations?
 A: To create a matrix that rotates in an arbitrary direction we will first choose a unit vector pointing in the direction of the rotation. The rotation will be in the direction $$A\hat{x}+B\hat{y}+C\hat{z}$$ such that $$A^2+B^2+C^2=1$$ We can rotate the coordinates about the x axis so that the rotation axis now sits in the xz plane. This is done with the matrix,
\begin{pmatrix}1&0&0\\0&cos\theta&sin\theta\\0&-sin\theta&cos\theta\end{pmatrix}
Where the angle between the axis and the xz plane is $$\theta = arctan(\frac{B}{C})$$ We can then rotate the coordinates such that the axis of rotation will lie along the z-axis by rotating the coordinates about the y-axis clockwise with the transformation
\begin{pmatrix}cos\phi&0&sin\phi\\0&1&0\\-sin\phi&0&cos\phi\end{pmatrix}
Where the angle between the rotation axis and the old z-axis is
$$\phi=arctan(\frac{A}{C'})$$ and C' is the new z component of the axis of rotation unit vector $$C'=-Bsin\theta + Ccos\theta$$
Since the axis lies in the z direction, we can now rotate the coordinates about the z axis then undo the two previous transformations with the inverse matricies. The angle rotated about z is the first angle in the set of three angles analogous to the Eulerian angles, call it $\gamma$. The rotation tranformation will then be
\begin{pmatrix}cos\gamma&sin\gamma&0\\-sin\gamma&cos\gamma&0\\0&0&1\end{pmatrix}
The inverses of the first two matricies are easy to calculate since proper rotation matricies have a determinant of 1. The matricies are then
\begin{pmatrix}cos\phi&0&-sin\phi\\0&1&0\\sin\phi&0&cos\phi\end{pmatrix}
for the inverse of the y-axis transformation and
\begin{pmatrix}1&0&0\\0&cos\theta&-sin\theta\\0&sin\theta&cos\theta\end{pmatrix}
for the inverse of the x-axis rotation. Combining all of these matricies gives a transformation with the elements
$$\lambda_{11} = cos\gamma\, cos^2\phi+sin^2\phi$$
$$\lambda_{12} = sin\theta\, cos\phi\, sin\phi\,(1-cos\gamma)+cos\theta\, cos\theta\, sin\gamma$$
$$\lambda_{13} = sin\theta \, cos\phi \, sin\gamma \, - cos\theta \,( cos\phi \, sin\phi \, - cos\phi \, sin\phi \, cos\gamma \, )$$
$$\lambda_{21} =  sin\theta \, cos\phi \, sin\phi \, - cos\phi \, ( cos\theta \, sin\gamma \, + sin\theta \, sin\phi \, cos\gamma \, ) $$
$$\lambda_{22} =  cos\gamma \, cos^2\theta \, + sin^2\theta \, cos^2\phi \, + cos\gamma \, sin^2\theta \, sin^2\phi \, $$
$$\lambda_{23} =  sin\theta \, ( cos\theta \, cos\gamma \, - sin\theta \, sin\phi \, sin\gamma \, ) - cos\theta \, ( sin\theta \, cos^2\phi \, + sin\phi \, ( cos\theta \, sin\gamma \, + sin\theta \, sin\phi \, cos\gamma \, ) ) $$
$$\lambda_{31} = - cos\phi \, ( sin\theta \, sin\gamma \, - cos\theta \, sin\phi \, cos\gamma \, ) - cos\theta \, cos\phi \, sin\phi \, $$
$$\lambda_{32} = cos\theta \, ( sin\theta \, cos\gamma \, + cos\theta \, sin\phi \, sin\gamma \, ) - sin\theta \, ( cos\theta \, cos^2\phi \, - sin\phi \, ( sin\theta \, sin\gamma \, - cos\theta \, sin\phi \, cos\gamma \, ) ) $$
$$\lambda_{33} = cos^2\theta \, cos^2\phi \, + cos\gamma \, cos^2\theta \, sin^2\phi \, + cos\gamma \, sin^2\theta \, $$
Which is quite a large matrix. It only gets worse since this is only the rotation for the first angle. Two more of these matricies will be applied for rotations about two other arbitrary axes. Once the three matricies are multiplied together, a new matrix with hundreds of terms is made. This problem is likely why the Eulerian angles and the Tait-Bryan angles only use rotations about the x,y,and z axes since those matricies mulitplied together gives a much simpler result. Those also have the added benefit of having physical significance as well.
The problem is still doable since the matricies can still be calculated and evaluated for the specific case of any three axes where any vector can be transformed to point in any other direction then rotated about that direction.
A simple way to check if it spans all possible orientations is If we apply the three transforms to the inverse of the Eulerian transforms put together, for any Eulerian angle $\phi, \theta, \text{ and }\psi$ there should also be a set of three angles that makes the multiplied matrix equal to the identity matrix which would mean that the transformation can redo any undone orientation shift which means that it can in turn also create that same orientation. If this condition is met, the chosen three axes span the space of possible orientations.
