Get complex angle when computing quaternions (Satellite attitude) I am modelling a, so far uncontrolled, satellite in MATLAB , along with its attitude. I have been researching about this matter and found that quaternions are the way to go , since they don't have singularities. My goal was to get angular velocities in each axis (body frame) from the satellite motion equations :
$$I \frac{\partial ^2\theta}{\partial t^2} + \frac{\partial \theta}{\partial t} \times I \frac{\partial \theta}{\partial t} = \sum \tau$$
From which I got the angular velocity in each axis for every step of my simulation.
I then computed the attitude quaternions using a formula I saw here :
$$ q_{new} = q_{old} + \frac{dt}{2} \cdot \omega \cdot q_{old} ,$$ in which $ \omega\ $ is the angular velocity.
It seemed to work well but when I computed the Euler angles, using this : 
\begin{align}
\phi &= atan2 (2(q_0q_1 + q_2q_3), 1- 2(q_1^2+q_2^2)) \\
\theta &= asin (2(q_0q_2-q_3q_1)) \\
\psi &= atan2 (2(q_0q_3+q_1q_2), 1-2 (q_2^2+q_3q_1))
\end{align}
I got complex results in the pitch axis (outside the domain of the $asin$ function), which means something is obviously wrong but I can't think of any solution besides trying randomly until something makes sense. I even tried to divide the angular velocity by its norm but it didn't help.
Does anyone know which of my steps is wrong? I hope I was clear enough.
Thanks in advance,
Hugo

Here's a reference to how the 3D game engine people do rotations using quaternions: https://www.3dgep.com/understanding-quaternions/
In particular, for this problem, I suggest working through SLERPing. I suspect you might find there is a better modelling solution by staying within quaternion arithmetic. 
Enjoy! (It's not my stuff, I just found it useful.)
