Back to index/Generated with Clerk from src/emmy/mechanics/rigid.cljc@f807468
(ns emmy.mechanics.rigid
(:refer-clojure :exclude [+ - * /])
(:require [emmy.calculus.derivative :refer [D]]
[emmy.function :as f]
[emmy.generic :as g :refer [sin cos + - * /]]
[emmy.matrix :as matrix]
[emmy.mechanics.lagrange :as L]
[emmy.mechanics.rotation :as r]
[emmy.quaternion :as q]
[emmy.structure :as s :refer [up]]))

Chapter 2: generalized coordinates to angular velocities

(defn three-vector-components->antisymmetric [[x y z]]
(matrix/by-rows
[0 (- z) y]
[z 0 (- x)]
[(- y) x 0]))
#object[emmy.mechanics.rigid$three_vector_components__GT_antisymmetric 0x6606c2ff "
emmy.mechanics.rigid$three_vector_components__GT_antisymmetric@6606c2ff"
]
(defn antisymmetric->column-matrix
"Given an antisymmetric matrix `a` of dimension 3, returns a column vector of
its positive components."
[a]
{:pre [(matrix/antisymmetric? a)]}
(matrix/column
(get-in a [2 1])
(get-in a [0 2])
(get-in a [1 0])))
#object[emmy.mechanics.rigid$antisymmetric__GT_column_matrix 0x5e0e4069 "
emmy.mechanics.rigid$antisymmetric__GT_column_matrix@5e0e4069"
]

Note from GJS: Suggested by Jack Wisdom on 30 Sept 2019.

(defn M-of-q->omega-of-t [M-of-q]
(fn [q]
(let [M-of-t (f/compose M-of-q q)]
(fn [t]
(antisymmetric->column-matrix
(* ((D M-of-t) t)
(matrix/transpose (M-of-t t))))))))
#object[emmy.mechanics.rigid$M_of_q__GT_omega_of_t 0x14d75fc2 "
emmy.mechanics.rigid$M_of_q__GT_omega_of_t@14d75fc2"
]
(defn M-of-q->omega-body-of-t [M-of-q]
(fn [q]
(fn [t]
(* (matrix/transpose (M-of-q (q t)))
(((M-of-q->omega-of-t M-of-q) q) t)))))
#object[emmy.mechanics.rigid$M_of_q__GT_omega_body_of_t 0x2d169aa0 "
emmy.mechanics.rigid$M_of_q__GT_omega_body_of_t@2d169aa0"
]
(defn M->omega [M-of-q]
(L/Gamma-bar
(M-of-q->omega-of-t M-of-q)))
#object[emmy.mechanics.rigid$M__GT_omega 0x7dd7a4c7 "
emmy.mechanics.rigid$M__GT_omega@7dd7a4c7"
]
(defn M->omega-body [M-of-q]
(L/Gamma-bar
(M-of-q->omega-body-of-t M-of-q)))
#object[emmy.mechanics.rigid$M__GT_omega_body 0x5176ae1c "
emmy.mechanics.rigid$M__GT_omega_body@5176ae1c"
]

Assuming omega-body is on principal axes, and A, B, C are the principal moments. Angular velocity to kinetic energy and angular momenta

(defn T-body [A B C]
(fn [[w0 w1 w2]]
(* (/ 1 2)
(+ (* A (g/square w0))
(* B (g/square w1))
(* C (g/square w2))))))
#object[emmy.mechanics.rigid$T_body 0x7f71d2bf "
emmy.mechanics.rigid$T_body@7f71d2bf"
]
(defn L-body [A B C]
(fn [[w0 w1 w2]]
(s/down (* A w0)
(* B w1)
(* C w2))))
#object[emmy.mechanics.rigid$L_body 0x66ab946b "
emmy.mechanics.rigid$L_body@66ab946b"
]
(defn L-space [M]
(fn [A B C]
(fn [omega-body]
(* ((L-body A B C) omega-body)
(g/transpose M)))))
#object[emmy.mechanics.rigid$L_space 0x269ae5e5 "
emmy.mechanics.rigid$L_space@269ae5e5"
]

Euler Angles

(defn Euler->omega [angles-path]
(fn [t]
(letfn [(M-on-path [t]
(r/Euler->M (angles-path t)))
(w-cross [t]
(* ((D M-on-path) t)
(g/transpose (M-on-path t))))]
(antisymmetric->column-matrix
(w-cross t)))))
#object[emmy.mechanics.rigid$Euler__GT_omega 0x42f8a964 "
emmy.mechanics.rigid$Euler__GT_omega@42f8a964"
]
(defn Euler->omega-body [angles-path]
(fn [t]
(* (g/transpose (r/Euler->M (angles-path t)))
((Euler->omega angles-path) t))))
#object[emmy.mechanics.rigid$Euler__GT_omega_body 0x336d6fae "
emmy.mechanics.rigid$Euler__GT_omega_body@336d6fae"
]

Assuming Euler angles rotate principal axes from reference orientation.

(defn Euler-state->omega-body
"Although this implementation appears to summarize `(M->omega-body r/Euler->M)`,
it is actually essential to prevent intermediate expression explosion."
[[_ [theta _ psi] [thetadot phidot psidot]]]
(let [omega-a (+ (* (sin psi) (sin theta) phidot)
(* (cos psi) thetadot))
omega-b (+ (* (cos psi) (sin theta) phidot)
(* -1 (sin psi) thetadot))
omega-c (+ (* (cos theta) phidot)
psidot)]
(up omega-a omega-b omega-c)))
#object[emmy.mechanics.rigid$Euler_state__GT_omega_body 0x9dfad64 "
emmy.mechanics.rigid$Euler_state__GT_omega_body@9dfad64"
]
(defn T-body-Euler [A B C]
(fn [local]
((T-body A B C)
(Euler-state->omega-body local))))
#object[emmy.mechanics.rigid$T_body_Euler 0x789bc477 "
emmy.mechanics.rigid$T_body_Euler@789bc477"
]
(def T-rigid-body
"Alias for [[T-body-Euler]]."
T-body-Euler)
#object[emmy.mechanics.rigid$T_body_Euler 0x789bc477 "
emmy.mechanics.rigid$T_body_Euler@789bc477"
]
(defn L-body-Euler [A B C]
(fn [local]
((L-body A B C)
(Euler-state->omega-body local))))
#object[emmy.mechanics.rigid$L_body_Euler 0x7faff189 "
emmy.mechanics.rigid$L_body_Euler@7faff189"
]
(def Euler-state->L-body
"Alias for [[L-body-Euler]]."
L-body-Euler)
#object[emmy.mechanics.rigid$L_body_Euler 0x7faff189 "
emmy.mechanics.rigid$L_body_Euler@7faff189"
]
(defn L-space-Euler [A B C]
(fn [local]
(let [angles (L/coordinate local)]
(* ((L-body-Euler A B C) local)
(g/transpose (r/Euler->M angles))))))
#object[emmy.mechanics.rigid$L_space_Euler 0x27e7e1cd "
emmy.mechanics.rigid$L_space_Euler@27e7e1cd"
]
(def Euler-state->L-space
"Alias for [[L-space-Euler]]."
L-space-Euler)
#object[emmy.mechanics.rigid$L_space_Euler 0x27e7e1cd "
emmy.mechanics.rigid$L_space_Euler@27e7e1cd"
]
(defn rigid-sysder [A B C]
(L/Lagrangian->state-derivative
(T-rigid-body A B C)))
#object[emmy.mechanics.rigid$rigid_sysder 0x5327ccd8 "
emmy.mechanics.rigid$rigid_sysder@5327ccd8"
]

Quaternion representation

(defn quaternion-state->omega-body [[_ q qdot]]
(let [two-q-norm (/ (* 2 q)
(g/dot-product q q))
omega**a (g/dot-product two-q-norm (* q/I-matrix qdot))
omega**b (g/dot-product two-q-norm (* q/J-matrix qdot))
omega**c (g/dot-product two-q-norm (* q/K-matrix qdot))]
(up omega**a omega**b omega**c)))
#object[emmy.mechanics.rigid$quaternion_state__GT_omega_body 0x753912b7 "
emmy.mechanics.rigid$quaternion_state__GT_omega_body@753912b7"
]

I'm not sure what these are in Quaternion land, so I'll leave them here as private for now.

(let [q:a (matrix/by-rows
[0 1 0 0]
[-1 0 0 0]
[0 0 0 1]
[0 0 -1 0])
q:b (matrix/by-rows
[0 0 1 0]
[0 0 0 -1]
[-1 0 0 0]
[0 1 0 0])
q:c (matrix/by-rows
[0 0 0 1]
[0 0 1 0]
[0 -1 0 0]
[-1 0 0 0])]
(defn quaternion-state->omega-space [[_ q qdot]]
(let [Q (matrix/up->column-matrix q)
QdotT (matrix/row* qdot)
two-m**2-inv (/ -2 (g/dot-product q q))
omega**x (* two-m**2-inv (get-in (* QdotT q:a Q) [0 0]))
omega**y (* two-m**2-inv (get-in (* QdotT q:b Q) [0 0]))
omega**z (* two-m**2-inv (get-in (* QdotT q:c Q) [0 0]))]
(up omega**x omega**y omega**z))))
#'emmy.mechanics.rigid/quaternion-state->omega-space
(defn qw-state->L-body [A B C]
(fn [[_ _ omega]]
((L-body A B C) omega)))
#object[emmy.mechanics.rigid$qw_state__GT_L_body 0x759a8620 "
emmy.mechanics.rigid$qw_state__GT_L_body@759a8620"
]
(defn qw-state->L-space [A B C]
(let [state->body (qw-state->L-body A B C)]
(fn [[_ q :as qw-state]]
(let [Lbody (state->body qw-state)
M (q/->rotation-matrix (q/make q))]
(* Lbody (g/transpose M))))))
#object[emmy.mechanics.rigid$qw_state__GT_L_space 0x6debf8d7 "
emmy.mechanics.rigid$qw_state__GT_L_space@6debf8d7"
]
(defn T-quaternion-state [A B C]
(fn [[_ q qdot]]
(let [Q (matrix/up->column-matrix q)
Qdot (matrix/up->column-matrix qdot)
m**2-inv (g/invert (get-in (* (matrix/transpose Q) Q) [0 0]))
x (* m**2-inv q/I-matrix Qdot)
y (* m**2-inv q/J-matrix Qdot)
z (* m**2-inv q/K-matrix Qdot)
M (* Q (matrix/transpose Q))]
(* 2 (+ (* A (get-in (* (matrix/transpose x) M x) [0 0]))
(* B (get-in (* (matrix/transpose y) M y) [0 0]))
(* C (get-in (* (matrix/transpose z) M z) [0 0])))))))
#object[emmy.mechanics.rigid$T_quaternion_state 0x47b8be "
emmy.mechanics.rigid$T_quaternion_state@47b8be"
]