(defn M->Euler
"Given a 3x3 rotation matrix, returns a [[emmy.structure/up]] of the
corresponding Euler angles.
Our Euler Angle convention:
M(theta, phi, psi) = R_z(phi)*R_x(theta)*R_z(psi)"
([M]
(M->Euler M nil))
([M tolerance-in-ulps]
(let [tolerance (if (nil? tolerance-in-ulps)
u/machine-epsilon
(* tolerance-in-ulps u/machine-epsilon))
close? (us/close-enuf? tolerance)
cx (get-in M [2 2])
cx-number? (v/number? cx)]
(cond (and cx-number? (close? cx -1))
(let [theta Math/PI
phi (- (g/atan
(- (get-in M [0 1]))
(get-in M [0 0])))
psi 0]
(up theta phi psi))
(and cx-number? (close? cx +1))
(let [theta 0
phi (g/atan
(- (get-in M [0 1]))
(get-in M [0 0]))
psi 0]
(up theta phi psi))
:else
(let [theta (g/acos cx)
sx (sin theta)
phi (g/atan (/ (get-in M [0 2]) sx)
(- (/ (get-in M [1 2]) sx)))
psi (g/atan (/ (get-in M [2 0]) sx)
(/ (get-in M [2 1]) sx))]
(up theta phi psi))))))