Alright, did the matrix math for ya. Ready?
Your complete rotation matrix for your system is:
(C_C)(C_B)+(S_A)(S_C)(S_B); (S_A)(S_C)(C_B)-(S_B)(C_C); (C_A)(S_C);;
(C_A)(S_B); (C_A)(C_B); -(S_A);;
(S_A)(S_B)(C_C)-(C_B)(S_C); (S_A)(C_B)(C_C)+(S_B)(S_C); (C_A)(C_C);;;
WHERE:
S_ represents sine, C_ represents cosine
A is first angle (rotation about global/local X) - doesn't matter which since global is local for first rotation.
B is second angle (rotation about global Z)
C is final angle (rotation about local Y)
So, you need to make the .to_a matrix of your rotation be equal to this matrix (after normalizing the axis vectors), then solve for A, B, and C. This matrix is very similar to rotation of Z, then X, then Y, but still different.