I have an object with certain rotation and angular velocity (α, β, γ) in 3D space, using Euler extrinsic rotation (x–y–z). Also, I have angular accelerations in global coordinate system.
When I numerically integrate this in Python, using the solve_ivp() function, the angular accelerations are interpreted as local angular accelerations. This results, for example, in exactly the wrong orientation when one of the axes is rotated 180 deg.
My solution would be to convert all rotations to quaternions. And multiply the derivative by dt, and 'quaternion add' this to the original angular velocity to get the new angular velocity.
There is however one problem with this: I do not know the time stepping as this varies. And this function requires me to return the derivatives with respect to time.
So the question is: how can I determine the local angular accelerations when I have the global rotation?