mathematical model of a system of rigid bodies In my quest for getting comfortable with the mathematical modelling of a system of rigid bodies, I am posing the following problem: 
Problem:
Suppose in outer space (no gravity) there are two rigid bodies as in the following figure: 
The rigid bodies are connected with one degree of freedom rotational joint about $\vec{i}_1$ axis.  There is a coordinate system attached in the center of mass of each rigid body. Suppose the following: $\begin{bmatrix} \vec{i}_1\\ \vec{j}_1\\ \vec{k}_1\end{bmatrix} = R_{i\psi} R_{j\chi} R_{k\theta} \cdot \begin{bmatrix} \vec{i}_0\\ \vec{j}_0\\ \vec{k}_0\end{bmatrix}$ and $\begin{bmatrix} \vec{i}_2\\ \vec{j}_2\\ \vec{k}_2\end{bmatrix} = R_{i\phi} \cdot \begin{bmatrix} \vec{i}_1\\ \vec{j}_1\\ \vec{k}_1\end{bmatrix}$, where $R = R_{i\psi} R_{j\chi} R_{k\theta}$ is a rotation matrix and $R_{i\phi}$ is the rotation matrix about $\vec{i}$ axis with angle $\phi$. On the first rigid body is a attached a motor in point $S$ in the figure. The motor's shaft is attached by the second rigid body, thus the motor's bearings are forming that "1 d.o.f" rotational joint. There are no other forces or torques acting on the system. I want to model the dynamics of this system. 
Modelling:
I am using D'Alembert principle.  The system has $3 + 3 + 1$ degrees of freedom. The variables I am choosing are $x,y,z$ the coordinates (in world frame) of the center of mass of first rigid body, $\psi, \chi, \theta$ the angles giving the orientation of the first rigid body and finally the angle $\phi$. Suppose the first rigid body has the inertia matrix $I_{01}$ in it's coordinate system and the second rigid body has the inertia matrix $I_{02}$ in it's coordinate system. D'Alembert's equations for a system of $N$ rigid bodies with $q_1, ... , q_n$ degrees of freedom are: 
\begin{align}
   \begin{cases}
   \sum_{i=1}^N \left( F_i^a(t) - \dot{P}_i(t)\right) \cdot \frac{\partial \bar{v}_i}{\partial \dot{q}_1} + \sum_{i=1}^m \left( T_i^a(t) - \dot{L}_i(t)\right) \cdot \frac{\partial \omega_i}{\partial \dot{q}_1} = 0 \\
   \vdots \\
   \sum_{i=1}^N \left( F_i^a(t) - \dot{P}_i(t)\right) \cdot \frac{\partial \bar{v}_i}{\partial \dot{q}_n} + \sum_{i=1}^m \left( T_i^a(t) - \dot{L}_i(t)\right) \cdot \frac{\partial \omega_i}{\partial \dot{q}_n} = 0
   \end{cases}
   \end{align} where $\bar{v}_i$ and $\omega_i$ are the linear respectively angular velocity of each rigid body and $F_i^a$ is the total applied force on the $i'th$ rigid body whereas $T_i^a$ is the total applied torque (about the center of mass) on the $i'th$ rigid body, all in world frame.
Question:
Suppose the motor in point $S$ applies the torque $R^T \cdot \begin{bmatrix} T_m\\0\\0\end{bmatrix}$ on the second rigid body. Am want to know if in the above presented case:


*

*$F_1^a = F_2^a = \begin{bmatrix} 0\\0\\0 \end{bmatrix}$ ?

*$T_1^a = -R^T \cdot \begin{bmatrix} T_m\\0\\0\end{bmatrix} = -T_2^a$ ?


Conclusion:
Therefore the equations of the motion will result upon replacing everything in D'Alembert's equations ... am I correct? I've seen (I think) solutions using the parallel axis theorem, which I do not understand. Why is the parallel axis theorem needed? Where is it used? It seems that every thing can be modeled without it ... Am I wrong?
 A: Consider a joint connecting bodies [1] and [2]. I use following expressions for the joint kinematics
$$ \begin{align}  
\mathbf{v}_2 &  = \mathbf{v}_1 + \mathbf{x} \dot{q} \\
\boldsymbol{\omega}_2 & = \boldsymbol{\omega}_1 + \mathbf{z} \dot{q} 
\end{align} $$
where $\mathbf{x} = \left( \frac{\partial \mathbf{v}}{\partial \dot{q}} \right) $ and $\mathbf{z}  = \left( \frac{\partial \boldsymbol{\omega}}{\partial \dot{q}} \right)$, and all linear velocities are measured at the top of the joint (on body [2]).
Equal and opposite forces $\mathbf{F}$ and moments $\mathbf{T}$ act on the two bodies from the joint. If the moments are expressed on the top of the joint also (just a velocity) then the following equation is used to either find the joint acceleration or the joint force
$$ Q = \mathbf{F} \cdot \mathbf{x} + \mathbf{T} \cdot \mathbf{z} $$
For example consider a revolute joint with $\mathbf{x}=0$ and $\mathbf{z} = \hat{k}$. Then any reaction force $\mathbf{F}$ is allowed, since $\mathbf{F}\cdot 0 = 0$. The reaction moment has to have a torque along the $\hat{k}$ direction of value $Q$ since $\mathbf{T} \cdot \hat{k} = Q$.
In addition, the reaction forces are not directly specified, but a result of a free body diagram and the equations of motion
$$ \left. \begin{align} 
  \mathbf{F}_1^a - \mathbf{F} & = \dot{\mathbf{p}}_1 & \mathbf{F}_2^a + \mathbf{F} & = \dot{\mathbf{p}}_2 \\
  \mathbf{T}_1^a - \mathbf{T} - \mathbf{r}_1 \times \mathbf{F} & = \dot{\mathbf{L}}_1 & \mathbf{T}_2^a + \mathbf{T} + \mathbf{r}_2 \times \mathbf{F} & = \dot{\mathbf{L}}_2
\end{align} \right\} $$
where $\mathbf{r}_i$ is the position of the joint relative to the body center of mass, and the superscript $\square^a$ indicates applied forces/moments.
In order for all to work you have to follow these rules


*

*All vectors in an expression need to be expressed on the same point and along the same coordinate system

*Equations involving momentum (and derivatives) need to be expressed at the center of mass only and along the same coordinate system as forces (typically the global system)

*Care must be taken to include any necessary transformations when velocities, angular momentum or moments are used at a point different from where they are defined. 

*Body MMOI matrix needs to be expressed at the center of mass, but along the global coordinates in each equation of motion. 

*Joint kinematics need to be differentiated in order to find the dependency of momentum to joint acceleration.

