(ns emmy.mechanics.rotation
(:refer-clojure :exclude [+ - * /])
(:require [emmy.generic :as g :refer [cos sin + - * /]]
[emmy.matrix :as matrix]
[emmy.structure :as s :refer [up]]
[emmy.util :as u]
[emmy.util.stream :as us]
[emmy.value :as v]))
(defn- rotate-x-matrix-2 [c s]
(matrix/by-rows [1 0 0]
[0 c (- s)]
[0 s c]))
#object[emmy.mechanics.rotation$rotate_x_matrix_2 0x65188cd3 "
emmy.mechanics.rotation$rotate_x_matrix_2@65188cd3"
]
(defn rotate-x-matrix
"Produce the matrix of a rotation of α radians about the x axis."
[α]
(rotate-x-matrix-2 (cos α) (sin α)))
#object[emmy.mechanics.rotation$rotate_x_matrix 0x55196bda "
emmy.mechanics.rotation$rotate_x_matrix@55196bda"
]
(def Rx-matrix rotate-x-matrix)
#object[emmy.mechanics.rotation$rotate_x_matrix 0x55196bda "
emmy.mechanics.rotation$rotate_x_matrix@55196bda"
]
(defn- rotate-y-matrix-2 [c s]
(matrix/by-rows [c 0 s]
[0 1 0]
[(- s) 0 c]))
#object[emmy.mechanics.rotation$rotate_y_matrix_2 0xf30e4d0 "
emmy.mechanics.rotation$rotate_y_matrix_2@f30e4d0"
]
(defn rotate-y-matrix
"Produce the matrix of a rotation of α radians about the y axis."
[α]
(rotate-y-matrix-2 (cos α) (sin α)))
#object[emmy.mechanics.rotation$rotate_y_matrix 0x1fa0432c "
emmy.mechanics.rotation$rotate_y_matrix@1fa0432c"
]
(def Ry-matrix rotate-y-matrix)
#object[emmy.mechanics.rotation$rotate_y_matrix 0x1fa0432c "
emmy.mechanics.rotation$rotate_y_matrix@1fa0432c"
]
(defn- rotate-z-matrix-2
"Produce the matrix of a rotation of α radians about the z axis."
[c s]
(matrix/by-rows [c (- s) 0]
[s c 0]
[0 0 1]))
#object[emmy.mechanics.rotation$rotate_z_matrix_2 0x3e7f3a1f "
emmy.mechanics.rotation$rotate_z_matrix_2@3e7f3a1f"
]
(defn rotate-z-matrix
"Produce the matrix of a rotation of α radians about the z axis."
[α]
(rotate-z-matrix-2 (cos α) (sin α)))
#object[emmy.mechanics.rotation$rotate_z_matrix 0x20a6c1b1 "
emmy.mechanics.rotation$rotate_z_matrix@20a6c1b1"
]
(def Rz-matrix rotate-z-matrix)
#object[emmy.mechanics.rotation$rotate_z_matrix 0x20a6c1b1 "
emmy.mechanics.rotation$rotate_z_matrix@20a6c1b1"
]
(defn angle-axis->rotation-matrix [theta [x y z]]
(let [colatitude (g/acos z)
longitude (g/atan y x)]
(* (rotate-z-matrix longitude)
(rotate-y-matrix colatitude)
(rotate-z-matrix theta)
(matrix/transpose (rotate-y-matrix colatitude))
(matrix/transpose (rotate-z-matrix longitude)))))
#object[emmy.mechanics.rotation$angle_axis__GT_rotation_matrix 0x28b8a711 "
emmy.mechanics.rotation$angle_axis__GT_rotation_matrix@28b8a711"
]

Rotation Tuples

(defn ^:no-doc rotate-x-tuple-2 [c s]
(matrix/m->s
(s/literal-down 'l 3)
(rotate-x-matrix-2 c s)
(s/literal-up 'r 3)))
#object[emmy.mechanics.rotation$rotate_x_tuple_2 0x359c2809 "
emmy.mechanics.rotation$rotate_x_tuple_2@359c2809"
]
(defn rotate-x-tuple [α]
(rotate-x-tuple-2 (cos α)
(sin α)))
#object[emmy.mechanics.rotation$rotate_x_tuple 0x72f44367 "
emmy.mechanics.rotation$rotate_x_tuple@72f44367"
]
(defn ^:no-doc rotate-y-tuple-2 [c s]
(matrix/m->s
(s/literal-down 'l 3)
(rotate-y-matrix-2 c s)
(s/literal-up 'r 3)))
#object[emmy.mechanics.rotation$rotate_y_tuple_2 0x64c9bfc1 "
emmy.mechanics.rotation$rotate_y_tuple_2@64c9bfc1"
]
(defn rotate-y-tuple [α]
(rotate-y-tuple-2 (cos α)
(sin α)))
#object[emmy.mechanics.rotation$rotate_y_tuple 0x9241c66 "
emmy.mechanics.rotation$rotate_y_tuple@9241c66"
]
(defn ^:no-doc rotate-z-tuple-2 [c s]
(matrix/m->s
(s/literal-down 'l 3)
(rotate-z-matrix-2 c s)
(s/literal-up 'r 3)))
#object[emmy.mechanics.rotation$rotate_z_tuple_2 0x5d2ca12b "
emmy.mechanics.rotation$rotate_z_tuple_2@5d2ca12b"
]
(defn rotate-z-tuple [α]
(rotate-z-tuple-2 (cos α)
(sin α)))
#object[emmy.mechanics.rotation$rotate_z_tuple 0xd932e92 "
emmy.mechanics.rotation$rotate_z_tuple@d932e92"
]

Rotation procedures

XXX: R[xyz] should not return an up; they should return a struct of the same shape they were given. But do rotations of covectors work that way? Maybe we should assert up-ness here rather than promise to be more general than we are.

(defn Rx
"Returns a function which rotates a vector α radians about the x axis."
[α]
(fn [[x y z]]
(let [c (cos α)
s (sin α)]
(up x
(- (* c y) (* s z))
(+ (* s y) (* c z))))))
#object[emmy.mechanics.rotation$Rx 0x248dbf56 "
emmy.mechanics.rotation$Rx@248dbf56"
]
(defn Ry
"Returns a function which rotates a vector α radians about the y axis."
[α]
(fn [[x y z]]
(let [c (cos α)
s (sin α)]
(up (+ (* c x) (* s z))
y
(- (* c z) (* s x))))))
#object[emmy.mechanics.rotation$Ry 0xd63408d "
emmy.mechanics.rotation$Ry@d63408d"
]
(defn Rz
"Returns a function which rotates a vector α radians about the z axis."
[α]
(fn [[x y z]]
(let [c (cos α)
s (sin α)]
(up (- (* c x) (* s y))
(+ (* s x) (* c y))
z))))
#object[emmy.mechanics.rotation$Rz 0x35fe8df3 "
emmy.mechanics.rotation$Rz@35fe8df3"
]

Aliases to match scmutils.

(def rotate-x Rx)
#object[emmy.mechanics.rotation$Rx 0x248dbf56 "
emmy.mechanics.rotation$Rx@248dbf56"
]
(def rotate-y Ry)
#object[emmy.mechanics.rotation$Ry 0xd63408d "
emmy.mechanics.rotation$Ry@d63408d"
]
(def rotate-z Rz)
#object[emmy.mechanics.rotation$Rz 0x35fe8df3 "
emmy.mechanics.rotation$Rz@35fe8df3"
]
(defn wcross->w [A]
(up (get-in A [1 2])
(get-in A [2 0])
(get-in A [0 1])))
#object[emmy.mechanics.rotation$wcross__GT_w 0x210dd260 "
emmy.mechanics.rotation$wcross__GT_w@210dd260"
]

Rotation Matrix to Euler Angles

(defn Euler->M
"Compute the rotation matrix from a 3-vector of Euler angles.
Our Euler Angle convention:
M(theta, phi, psi) = R_z(phi)*R_x(theta)*R_z(psi)"
[[theta phi psi]]
(* (rotate-z-matrix phi)
(rotate-x-matrix theta)
(rotate-z-matrix psi)))
#object[emmy.mechanics.rotation$Euler__GT_M 0x68b3ae18 "
emmy.mechanics.rotation$Euler__GT_M@68b3ae18"
]

Ported from code added to scmutils by GJS, 28 Sept 2020.

(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)) ;; Nonunique
(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)) ;; Nonunique
(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))))))
#object[emmy.mechanics.rotation$M__GT_Euler 0x7f43fcab "
emmy.mechanics.rotation$M__GT_Euler@7f43fcab"
]