how to use quaternions for timestepping instead of matrices?

[iko79]

hi everybody,

i'm kind of new in the area of physics simulation, thus i've got a humble beginner question:

according to baraff's article "An Introduction to Physically Based Modelling: Rigid Body Simulation I -- Unconstrained Rigid Body Simulation" (which currently seems to be offline but you all surely know) i built a very basic simulation system using Quaternions. baraff starts using matrices and replaces them by quaternions in a later chapter.

his usage of these quaternions leads to my question - since the pdf's are not available i will shortly explain what he does, so skip this if you know the article:
<baraff>
he iterates through each body state Y(t) in the phase, each one consisting of

• Vector x(t)
• Quaternion q(t)
• Vector P(t)
• Vector L(t),

adding the derivative phase space vector dY/dt, scaled by the stepsize, to each of them.
</baraff>

scaling the derivative by the timestep and adding to the current phase space vector makes sense for 3d vectors, but not for quaternions i think. to just scale a quaternion by a factor and adding this scaled quaternion to the current one to get the new rotation, then just normalizing the result seems weird to me. espacaly as i read it's necessary to use slerp to interpolate beween two rotations. additionally, to apply a rotation a multiplication of quaternion is required.

however, the simulation looks fine but maybe it's because of the small timesteps the error doesn't visually attend. or maybe the rotations *should* perform a little different

does anybody of you have a answer (maybe coupled with an illuminating but due to lack of time not too bulky) paper/article/link?
thx!

[melax]

Dont think of the quaternion derivative as a unit length quaternion. Think of it as a tangent on the sphere. in this case a 4D sphere. And the derivative is used by adding it not multiplying it.Notice that the quaternion derivative is orthogonal to orientation.

Say s is the spin (sx,sy,sz,0) and q is the current orientation. Then the derivative d is s*q/2.0 where s and q follow quaternion multiplication.
Notice that:
dot(d,q) == dot(s*q/2,q)==
((s.w*q.x + s.x*q.w + s.y*q.z - s.z*q.y)*q.x +
(s.w*q.y - s.x*q.z + s.y*q.w + s.z*q.x)*q.y +
(s.w*q.z + s.x*q.y - s.y*q.x + s.z*q.w)*q.z +
(s.w*q.w - s.x*q.x - s.y*q.y - s.z*q.z)*q.w )/2
==
( + s.x*q.x*q.w + s.y*q.x*q.z - s.z*q.x*q.y
- s.x*q.y*q.z + s.y*q.y*q.w + s.z*q.x*q.y
+ s.x*q.y*q.z - s.y*q.x*q.z + s.z*q.z*q.w
- s.x*q.x*q.w - s.y*q.y*q.w - s.z*q.z*q.w )/2
==0
since everything cancels out. some terms are ommitted from the above since s.w is known to be 0.
What this means is that d is orthogonal to q. So if you start at q and move in direction d you will be travelling along the unit sphere. Now if you go too far you will start moving away from the sphere since it is curved. Hence the need to renormalize to prevent drift.