Solving Rigid Body ODEs in rotating frame I'm trying to solve numerically the motion of a rigid body in a fluid. To this end I have to solve the fluid PDEs, and at every time step the 6 ODEs of the rigid motion degrees of freedom:
$$(m\mathcal{I} + \mathcal{A}) \frac{d\mathbf{U}}{dt} + \boldsymbol{\Omega} \times ((m\mathcal{I}+\mathcal{A}) \mathbf{U}) = m_1 \mathbf{g} + \mathbf{F}\\
(\mathcal{J}+\mathcal{D}) \frac{d\Omega}{dt} + \boldsymbol{\Omega} \times ((\mathcal{J}+\mathcal{D})\boldsymbol{\Omega}) + \mathbf{U}\times(\mathcal{A}\mathbf{U}) = \mathbf{M} $$
where $\mathcal{A}, \mathcal{D}$ model the added mass effects of the body. 
These ODEs are expressed in a system of axes that rotates over time the same way as the body axes, but has a fixed origin. Therefore I need the components of the gravity vector $\mathbf{g}$ in this relative system. My first attempt is to say that since
$$\mathbf{0} = \frac{d}{dt}\mathbf{g} = \left( \frac{d \mathbf{g}}{dt}\right)_r + \boldsymbol{\Omega}\times \mathbf{g}$$
is the relation between the time derivatives in the fixed and rotating frame, I can simply add to the 6 ODEs of the body above the other 3 ODEs
$$ \dot{\mathbf{g}}_r = - \boldsymbol{\Omega}\times \mathbf{g}.  $$
Therefore, every time step, I update the values of the force $\mathbf{F}$ and torque $\mathbf{M}$ and advance the 9 ODE system with a Runge Kutta method, and I get the $\mathbf{g}$ components in the rotated reference for free with the new $\mathbf{\Omega}$ and $\mathbf{U}$. 
Is it correct though? What is the relation that gives the transformation of the components of $\mathbf{g}$ from the laboratory reference to the rotated one at every time step? What are the angles that I'm using here? Because I know for example that Euler angles are generally used in these situations, but I haven't introduced any convention here, and I'm not worrying about the order in which I do the three successive rotations. 
Finally, is it better to use quaternions in such problems, to avoid singularities?
 A: How to simulate
$$(m\mathcal{I} + \mathcal{A}) \frac{d\mathbf{U}}{dt} + \boldsymbol{\Omega} \times ((m\mathcal{I}+\mathcal{A}) \mathbf{U}) = m_1 \mathbf{g} + \mathbf{F}\\
(\mathcal{J}+\mathcal{D}) \frac{d\Omega}{dt} + \boldsymbol{\Omega} \times ((\mathcal{J}+\mathcal{D})\boldsymbol{\Omega}) + \mathbf{U}\times(\mathcal{A}\mathbf{U}) = \mathbf{M} $$
we want to simulate those equations in body fixed system (B-System), thus all vectors components must be given in B-system.
the transformation matrix $R$ between B-system and  inertial system (I-System)  can be build out of three matrices.
$$R_x(\phi)= \left[ \begin {array}{ccc} 1&0&0\\ 0&\cos \left( 
\phi \right) &-\sin \left( \phi \right) \\ 0&\sin
 \left( \phi \right) &\cos \left( \phi \right) \end {array} \right] 
$$
$$R_y(\theta)=\left[ \begin {array}{ccc} \cos \left( \theta \right) &0&\sin \left( 
\theta \right) \\ 0&1&0\\ -\sin
 \left( \theta \right) &0&\cos \left( \theta \right) \end {array}
 \right] 
$$
and
$$R_z(\psi)=\left[ \begin {array}{ccc} \cos \left( \psi \right) &-\sin \left( 
\psi \right) &0\\ \sin \left( \psi \right) &\cos
 \left( \psi \right) &0\\0&0&1\end {array} \right] 
$$
where $\phi\,,\theta\,,\psi$ are the Euler angles
for example
$$R=R_z(\phi)\,R_y(\theta)\,R_z(\psi)\tag 1$$
or any other combination, each rotation matrix has singularity in one of the Euler angle. 
The vectors components in B-system:
$$\vec{g}\mapsto R^T\,\vec{g}$$
$$\vec{F}\mapsto R^T\,\vec{F}$$
$$\vec{M}\mapsto R^T\,\vec{M}$$
with:
$$\dot{R}=R\,\left[ \begin {array}{ccc} 0&-\Omega_{{z}}&\Omega_{{y}}
\\ \Omega_{{z}}&0&-\Omega_{{x}}\\ 
-\Omega_{{y}}&\Omega_{{x}}&0\end {array} \right] 
$$
thus:
$$\vec{\Omega}=\underbrace{\left[ \begin {array}{ccc} -\sin \left( \theta \right) \cos \left( 
\psi \right) &\sin \left( \psi \right) &0\\ \sin
 \left( \theta \right) \sin \left( \psi \right) &\cos \left( \psi
 \right) &0\\ \cos \left( \theta \right) &0&1
\end {array} \right]}_{J_R} 
\,\underbrace{\begin{bmatrix}
  \dot{\phi} \\
  \dot{\theta} \\
  \dot{\psi}\\
\end{bmatrix}}_{\vec{\dot{\varphi}}}\tag 2$$
to see where the singularity is, you invert the matrix $J_R$
$$J_R^{-1}=\left[ \begin {array}{ccc} -{\frac {\cos \left( \psi \right) }{\sin
 \left( \theta \right) }}&{\frac {\sin \left( \psi \right) }{\sin
 \left( \theta \right) }}&0\\\sin \left( \psi
 \right) &\cos \left( \psi \right) &0\\ {\frac {\cos
 \left( \theta \right) \cos \left( \psi \right) }{\sin \left( \theta
 \right) }}&-{\frac {\cos \left( \theta \right) \sin \left( \psi
 \right) }{\sin \left( \theta \right) }}&1\end {array} \right] 
$$
so the singularity in this case is for $\theta=0$.
with equation (2) you obtain:
$$\vec{\dot{\Omega}}=J_R\vec{\ddot{\varphi}}+\dot{J}_R\,\vec{\dot{\varphi}}\tag 3$$
put  equation (2) and (3) in your ODE's , and multiply the second equation by $J_R^T$ you get 6 differential equations 
$$\frac{d\vec U}{dt}=\ldots$$
$$\frac{d^2\vec{\varphi}}{dt^2}=\ldots$$
to do the numerical simulation, you must transfer those ODE's to first order differential equations $\vec{\dot{y}}=\vec{f}(\vec{y})$
remarks:
if you want other singularity, you can change the combination of your transformation matrix $R$, for example $R=R_x(\phi)\,R_y(\theta)\,R_z(\psi)$
Edit 
other why to simulate :
from equation (2) you get:
$$\vec{\dot{\varphi}}=J_R^{-1}\,\vec{\Omega}$$
now all your ODE's are first order
$$\vec{\dot{U}}=\ldots$$
$$\vec{\dot{\Omega}}=\ldots$$
$$\vec{\dot{\varphi}}=\ldots$$
