Let us consider the movement of a rigid object in space, for example the lower arm of a person with a fixed wrist and an empty hand. Once we have “somehow” measured the movement of this lower arm in space, how can we find out which forces have to act on the elbow in order to obtain the given movement in space? (This example is the 3D equivalent of the “Basic link-segment equations - the free body diagram” in Chapter 5.1 of Winter, 2009).

# Equations of Motion¶

Here the following notations is used:

 $$\vec{a}$$ acceleration of COM (= center of mass) $$\vec{R}$$ reaction forces acting at the end of a limb, e.g. a joint $$\vec{g}$$ gravity, always pulling downward with 9.81 $$m/s^2$$ m mass of the object $$\vec{M}$$ net muscle moment acting at the joints $$\Theta$$ inertial matrix $$\frac{{d\vec \omega }}{{dt}}$$ angular acceleration

The following equations describing the movement of our object must be fulfilled:

$\sum {\vec F = m\vec a}$

(This is Newton’s Second Law.)About the center of mass,

$\sum {\vec M} = {\bf{\Theta }} \cdot \frac{{d\vec \omega }}{{dt}}$

For limb movements, only two types of forces are acting on the body: mass forces, i.e. gravity $$m \vec{g}$$ ; and the proximal and distal joint reaction forces $${\vec R_{prox}}\,and\,{\vec R_{dist}}$$ , i.e. forces acting on the joints.

${\vec R_{prox}} + {\vec R_{dist}} + m\vec g = m\vec a$

Once we know the movement of the object, i.e. $$\vec{a}$$ and $$\frac{d\vec{\omega}}{dt}$$, these formulas can be used to calculate reaction forces at the joints, muscle moments, and the work required to move the limb. In this chapter we will investigate in principle how these parameters can be found. Specific tracking systems, e.g. optical or inertial based systems, will be dealt with in the next chapter.

# Rotation and Translation¶

You can combine a rotation and a translation, by introducing a $$4^{th}$$ component in your notation. If $$\vec{p}$$ indicates the starting location of an object, R the rotation of an object, $$\vec{t}$$ the translation, and $$\vec{p'}$$ the new location, then

$\begin{split}\left[ {\begin{array}{*{20}{c}} {\vec p'}\\ 1 \end{array}} \right] = \left[ {\begin{array}{*{20}{c}} {\bf{R}}&{\vec t}\\ 0&1 \end{array}} \right] \cdot \left[ {\begin{array}{*{20}{c}} {\vec p}\\ 1 \end{array}} \right]\end{split}$

The matrix $$\left[ {\begin{array}{*{20}{c}} {\bf{R}}&{\vec t}\\ 0&1 \end{array}} \right]$$ is called spatial transformation matrix, and $$\left( {\begin{array}{*{20}{c}} {\vec p}\\ 1 \end{array}} \right)$$ are sometimes referred to as homogeneous coordinates. Note that this is just a formal trick to incorporate rotation and translation into a single mathematical equation.

# Linear velocities¶

Let us start out with something simple: the translational components of movements. Translations are always commutative: I always arrive at the same position, regardless in which sequence I execute the translations:

(1)${\vec r_x} + {\vec r_y} = {\vec r_y} + {\vec r_x}$

This has important consequences. For example, I can deal with the separate movement components individually. For example, I can calculate when a bullet shot from a gun will have fallen down, without knowing how fast it is flying forward.

Another consequence is that velocity and acceleration for the translational movement components can be calculated individually:

$\begin{split}\begin{array}{l} \vec v = \frac{{d\vec x}}{{dt}}\\ \vec a = \frac{{{d^2}\vec x}}{{d{t^2}}}\\ with\,{v_i} = \frac{{d{x_i}}}{{dt}},\,and\,{a_i} = \frac{{{d^2}{x_i}}}{{d{t^2}}}\,(i = 1,2,3) \end{array}\end{split}$

This will not be the case for rotations!

# Angular Velocity¶

## Determination of Angular Velocity¶

As pointed out previously, rotation vectors describing the orientation of an object depend on the choice of the reference orientation. In contrast, the angular velocity does not depend on the reference orientation, since it only describes the movement from the current orientation to the next, which does not involve the reference orientation.

### Quaternions¶

The simplest formula describing the angular velocity $$\vec{\omega}$$ is given in the quaternion-notation :

(2)$\omega = 2*\frac{d}{{dt}}q \circ {q^{ - 1}}$

where $$\omega = (0, \vec{\omega})$$ is a quaternion, with $$\vec{\omega}$$ the common eye velocity vector. Note that the angular velocity depends not only on the time-derivative $$\frac{{dq}}{{dt}}$$ of the orientation, but also on the current object orientation q itself.

### Angular Velocity Tensor¶

Coming from rotation matrices, the angular velocity cannot be described as a simple vector. Instead we have to express it as a matrix, which is sometimes referred to as angular velocity tensor. Any vector $$\vec r$$ that rotates around an axis with an angular speed vector $$\vec \omega$$ satisfies:

$\frac {d \vec r(t)} {dt} = \vec{\omega} \times\vec{r}$

We can introduce here the emph{angular velocity tensor} associated to the angular speed $$\omega$$:

$\begin{split}\mathbf{\Omega}(t) = \begin{pmatrix} 0 & -\omega_z(t) & \omega_y(t) \\ \omega_z(t) & 0 & -\omega_x(t) \\ -\omega_y(t) & \omega_x(t) & 0 \\ \end{pmatrix}\end{split}$

This tensor $$\mathbf{\Omega}(t)$$ will act as if it were a $$(\vec \omega \times)$$ operator :

$\vec \omega(t) \times \vec{r}(t) = \mathbf{\Omega}(t) \vec{r}(t)$

Given the orientation matrix $$\mathbf{R}(t)$$ of a frame, we can obtain its instant angular velocity tensor $$\mathbf{\Omega(t)}$$ as follows. We know that:

$\frac {d \vec r(t)} {dt} = \mathbf{\Omega}(t) \cdot \vec{r}$

As angular speed must be the same for the three vectors of a rotating frame, if we have a matrix $$\mathbf{R}(t)$$ whose columns are the vectors of the if we have a matrix $$\mathbf{R}(t)$$ whose columns are the vectors of the frame, we can write for the three vectors as a whole:

$\frac {d\mathbf{R}(t)} {dt} = \mathbf{\Omega}(t) \cdot \mathbf{R} (t)$

And therefore the angular velocity tensor we are looking for is:

$\mathbf{\Omega}(t) = \frac {d\mathbf{R}(t)} {dt} \cdot \mathbf{R}^{-1}(t)$

### Rotation Vectors¶

Expressed in rotation vectors , Eqn. (2) is equivalent to

(3)$\vec \omega \,\, = \,\,2\,*\,\,\frac{{\frac{d}{{dt}}\vec r + \vec r \times \frac{d}{{dt}}\vec r}}{{1 + {{{\rm{\vec r}}}^{\rm{2}}}}}$

### Fick Angles¶

A more complex formula is required if angular velocity is expressed in Fick-angles :

(4)$\begin{split}\vec \omega = \left( {\begin{array}{*{20}{l}} {{\textstyle{d \over {dt}}}{\psi _F}\,\,*\,\,\cos ({\theta _F}{\rm{)*}}\,\,{\rm{cos(}}{\varphi _F})}& - &{\,{\textstyle{d \over {dt}}}{\varphi _F}*\,{\rm{sin(}}{\theta _F})\,}\\ {{\textstyle{d \over {dt}}}{\varphi _F}\,\,\,*\,{\rm{cos(}}{\theta _F})}& + &{\,{\textstyle{d \over {dt}}}{\psi _F}*\,\sin ({\theta _F})\,*\,{\rm{cos(}}{\varphi _F})}\\ {{\textstyle{d \over {dt}}}{\theta _F}}& - &{{\textstyle{d \over {dt}}}{\psi _F}\,\,*\,{\rm{sin(}}{\varphi _F})} \end{array}} \right)\,\end{split}$

Eqns. (2) - (4) are equivalent, as they express the same angular velocity in different coordinate systems. The time-derivatives of the orientation coordinates - $$\frac{d}{{dt}}q$$ for quaternions, $$\frac{d}{{dt}}\vec r$$ for rotation vectors, and

$\begin{split}\left[ {\begin{array}{*{20}{c}} {\frac{d}{{dt}}{\theta _{Fick}}}\\ {\frac{d}{{dt}}{\phi _{Fick}}}\\ {\frac{d}{{dt}}{\psi _{Fick}}} \end{array}} \right] = \frac{d}{{dt}}{\left[ {\begin{array}{*{20}{c}} \theta \\ \phi \\ \psi \end{array}} \right]_{Fick}}\end{split}$

for Fick-angles - are often referred to as coordinate velocity. This coordinate velocity obviously depends on the parameters chosen to describe the orientation. In contrast, the angular velocity vector $$\vec{\omega}$$ describes the actual movement of the object, with its axis given by the instantaneous axis of rotation, and its length by the angular velocity of this rotation, and does not depend on the parametrization of the orientation. The preceding formulas also show that the coordinate velocity is in general not equivalent to the angular velocity $$\vec{\omega}$$.

### Relation Angular Velocity - Angular Orientation¶

The non-commutativity of rotations can lead to seemingly counterintuitive behaviors. For example, the axis of the vector describing the rotation from one orientation (e.g. our start orientation) to the next (e.g. our end orientation) is also the axis of eye-velocity for that movement. However, due to the non-commutativity of rotations, the angular velocity axis of a movement does in general not coincide with the axes characterizing the orientation of an object (see Fig. below)!

## Conversion of Angular Velocity into Orientation¶

### Determination of Position from Velocity¶

To understand how angular velocity can be converted into orientation in space, let us first remember how linear velocity can be converted into linear position. Let us start with the position $${\vec x_0} = \vec x({t_0})$$ , and move around with a velocity $$\vec{v}(t)$$. Then the current position is

(5)$\vec x(t) = {\vec x_0} + \int\limits_{{t_0}}^t {\vec v(t)dt}$

When working with computers, we cannot perform the integral exactly, but only approximately. Splitting the time between $$t_0$$ and $$t$$ into n equal elements with width $$\Delta t$$, we get

(6)$\begin{split}\begin{array}{l} \vec x(t) \approx {{\vec x}_0} + \Delta {{\vec x}_1} + \Delta {{\vec x}_2} + .... + \Delta {{\vec x}_n}\\ with\,\Delta {{\vec x}_i} = \vec v({t_i}) \cdot \Delta t \end{array}\end{split}$

Here it is important to note that due to the commutativity of translations, the sequence of additions has no effect

${\vec x_0} + \Delta {\vec x_1} + \Delta {\vec x_2} + .... + \Delta {\vec x_n} = \Delta {\vec x_n} + \Delta {\vec x_{n - 2}} + .... + \Delta {\vec x_1} + {\vec x_0}$

### Determination of Orientation from Angular Velocity¶

Now how can we break down a continuous rotational movement with angular velocity $$\vec{\omega}(t)$$ into small steps? For a short duration $$\Delta t$$ the axis of rotation is approximately constant. During this time, the object rotates about the axis $$\vec n(t) = \frac{{\vec \omega (t)}}{{\left| {\vec \omega (t)} \right|}}$$ . The angle about which the object rotated during this period is $$\Delta \phi = \left| {\vec \omega (t)} \right| \cdot \Delta t$$ .

From Eq [eq:quatProperties] we know that a rotation about $$\vec{n}$$ by an angle of $$\Delta \phi$$ can easiest be described with a quaternion vector $$\vec{q}$$ , with its direction given by $$\vec{n}$$ , and with the length $$\sin (\frac{{\Delta \phi }}{2})$$ . So with

(7)$\Delta {\vec q_i} = \vec n(t)\sin (\frac{{\Delta \phi ({t_i})}}{2}) = \frac{{\vec \omega ({t_i})}}{{\left| {\vec \omega ({t_i})} \right|}}\sin \left( {\frac{{\left| {\vec \omega ({t_i})} \right|\Delta t}}{2}} \right)$

we get

(8)$q(t) = \Delta q({t_n}) \circ \Delta q({t_{n - 1}}) \circ ... \circ \Delta q({t_2}) \circ \Delta q({t_1}) \circ q({t_0})$

Note that here the sequence is important, and cannot be reversed!