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
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.)