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
 A: I'm somewhat leery of taking the statements in the source of the equation at face value. I don't think they're explained very well.
Extracting a factor of $q_\text{old}$ on the right-hand side of the equation you used, we get
$$ q_\text{new} = q_\text{old} + \frac t2 \mathbf\omega \,q_\text{old}
= \left(1 + \frac t2 \mathbf\omega\right) q_\text{old}. $$
The composition of rotations in quaternions (what you're trying to do here) is accomplished by multiplication of one unit quaternion (in this case $q_\text{old}$)
by another unit quaternion.
In the equation above we multiply $q_\text{old}$ by $1 + \frac t2 \mathbf\omega,$
which is a quaternion but not a unit quaternion (unless $t=0$ or $\mathbf\omega=0$).
I think the exact formula for a constant rotation $\mathbf\omega$ for time $t$ is actually
$$ q_\text{new} = \exp\left(\frac t2 \mathbf\omega\right) q_\text{old}, $$
using the exponential function $\exp(\cdot).$
The expression $1 + \frac t2 \mathbf\omega$ is a truncated Taylor series of 
$\exp\left(\frac t2 \mathbf\omega\right),$ that is, a first-order approximation.
If $\frac t2 \mathbf\omega$ is very small then it's a reasonably good approximation,
so if you take very small steps and renormalize the result (scale it down to make it a unit quaternion again) after every step, you'll get better results.

An alternative approximation of the exponential function
(adapted from the comments) is
$$
\frac{1 + (t/4) \mathbf\omega}{1 − (t/4) \mathbf\omega}
=\frac{(1 + (t/4) \mathbf\omega)^2}
      {1+\lvert \mathbf\omega \rvert^2\,t^2/16}.
$$
This is mathematically a unit quaternion, so you should not have to "normalize" the result. Be aware, however, that due to the inherent inaccuracies of floating-point arithmetic, when you do this on a computer you may get a result for a computation such as $2(q_0q_2 - q_3q_1)$ that is ever so slightly greater than $1$ or less than $-1.$ You can set $\theta$ to $\pm\frac\pi2$ in those cases. (On the other hand if you get something like $1.01$ then there is probably still some error in the formulas or their implementation.)
A: You may have to normalize your quaternion, dividing it by its absolute value at each stage so you always have $q_0^2+q_1^2+q_2^2+q_3^2=1$ after this division.
Then the argument of the arcsin function should stay in range.  The constrained extremal values of 
$f=2(q_0q_2-q_1q_3)$
under the normalization constraint $q_0^2+q_1^2+q_2^2+q_3^2=1$ is found via Lagrange multipliers to occur when all $q_2=\pm q_0, q_3=\pm q_1$, and this forces the extrema of $f=\pm1$.
A: 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.)
