Calculating torque to set rotation in one frame

[pHySiQuE]

I made a function that uses physics forces to position an object at a specified point:

Code: Select all

   void Entity::PhysicsSetPosition(const float& x, const float& y, const float& z, const float& strength)
   {
      if (body)
      {
         Vec3 d;
         d.x = x - body->nextmat[3][0];
         d.y = y - body->nextmat[3][1];
         d.z = z - body->nextmat[3][2];
         d -= GetVelocity()/60.0;
         d *= 3600.0 * GetMass() * strength;// multiply by 60*60 because we want the movement to occur in one frame, and we are trying to set the velocity in meters per second
         body->force = Vec3(0);
         body->AddForce(d);
      }
   }

I'm trying to do the same for rotation. My function works perfectly, except I can't figure out how to negate the body's existing angular velocity:

Code: Select all

   void Entity::PhysicsSetRotation(const float& pitch, const float& yaw, const float& roll, const float& strength)
   {
      if (body)
      {
         Quat currentrotation = body->nextquat;
         Quat desiredrotation = Quat(Vec3(pitch,yaw,roll));
         Quat omega = GetOmega();//what to do with this?
         Vec3 torque = desiredrotation.GetAverageOmega(currentrotation,strength);
         torque *= 60.0 * GetMass();
         body->AddTorque(torque);
      }
   }

Any ideas how this is done?

[JoeJ]

If you have the wanted angular velocity to reach your target - i think you get that by:
Vec3 targetAngVel = desiredrotation.GetAverageOmega(currentrotation,strength), but i#m not sure what this function does -
you can convert it to torque like this:

Code: Select all

vector targetAngVel = the angular velocity you want to reach the targe orientation;
vector curAngVel = BodyGetOmega()

targetAngVel -= curAngVel * 1.0; // take actual angvel in account - something between 0.0 and 1.0; recommended >0.5

vector torque = targetAngVel / timestep;
torque = matrix.Unrotate (torque); // matrix = BodyGetMatrix
torque[0] *= Ixx; // Inertia from BodyGetMassMatrix
torque[1] *= Iyy;
torque[2] *= Izz;
torque = matrix.Rotate (torque);

torque *= 0.1; // something between 0.0 and 1.0; recommended <0.5, high value could generate overspin



[JoeJ]

...there must be something wrong with my very first assumption, because i miss the timestep.
So this is how i get the targetAngVel, which basicaly only is a axis / angle convertion of the rotational error and dividing angle with timestep:

sVec3 targetAngVel = AngVelFromAToB (bodyCurrentRotation, bodyTargetRotation) / timestep;

Note that my Quat is (x,y,z,w), while newton gives (w,x,y,z)

Code: Select all

inline sVec3 AngVelFromAToB (const sQuat &qA, const sQuat &qB)
{

   const float matchTolerance = 0.01f;
   const float faultTolerance = 0.05f;

   sQuat q = QuatFromAToB (qA, qB);
   sVec3 *omegaDir = (sVec3*)&q;
   float sql = omegaDir->SqL();
   if (sql   < (matchTolerance * matchTolerance))
      return sVec3 (0, 0, 0);

   float length = sqrt (sql);
   if (q[3] < -1) q[3] = -1;
   if (q[3] >  1) q[3] =  1;
   sVec3 angVel = (*omegaDir / length) * (2.0 * acos(q[3]));
   if (length   < faultTolerance) angVel *= (length - matchTolerance) / (faultTolerance - matchTolerance);
   return angVel;

}

// ensure to take the shorter arc...
inline sQuat QuatFromAToB (const sQuat &qA, const sQuat &qB)
{
   sQuat q;
   if (qA.Dot(qB) < 0.0f)
   {
      q[0] = qA[0]; q[1] = qA[1]; q[2] = qA[2];
      q[3] = -qA[3];
   }
   else
   {
      q[0] = -qA[0]; q[1] = -qA[1]; q[2] = -qA[2];
      q[3] = qA[3];
   }
         
   return qB * q;
}


[pHySiQuE]

I don't think this is going to negate the angular velocity:

Code: Select all

targetAngVel -= curAngVel * 1.0; // take actual angvel in account - something between 0.0 and 1.0; recommended >0.5

[JoeJ]

I don't think this is going to negate the angular velocity:

I don't understand what you mean, that's no neagtion - if you're no C guy it's the same as:

targetAngVel = targetAngVel - curAngVel * 1.0;

curAngVel is the current angular velocity of the body.
Example: If current angular velocity would be exactly the velocity needed to reach the target, the calculation leads to zero additional torque, because in that case there's nothing to do.

[pHySiQuE]

Angles aren't like linear vectors. For example, if you want to interpolate halfway between two euler rotations, you can't just use a linear interpolation. You have to convert to Quaternions and perform a spherical linear interpolation.

[JoeJ]

But in this case you do not want to do a slerp - you want to subtract current from wanted angular velocity, which is as easy as that.
I still don't understand your doupts - the only thing you missed is to use inertia for rotation like mass for translation.

I've figured out that code myself, but use it in lots of situations and answered 'control orientation' questions like yours
at least 10 times with it, it works for me an others.
I suggest you try it out and if there's a problem post your code...
If i do not completely misunderstand your goal i'm sure it is exactly what you want and correct (if there's no typo).

[JoeJ]

...hmm, the more i think on it, the more interesting your doupts become.
I've always assumed that interpolating 2 rotations in vector format (like angular velocity used in physics) is the same as
linear (not spherical) interpolation between those 2 rotations in quaternion format.

Now i've found that this is not true using this code:

Code: Select all


float t = 0.25; // interpolation factor

   sVec3 d0 (1,0,0); float a0 = 0.3;
   sVec3 d1 (0,1,0); float a1 = 1.3;

   sVec3 da0 = d0 * a0;
   sVec3 da1 = d1 * a1;

   sVec3 proofV = da0 * (1.0-t) + da1 * t; // linear vector interpolation


   sQuat q0, q1;
   q0.FromAxisAndAngle (d0, a0);
   q1.FromAxisAndAngle (d1, a1);

   sQuat qLerp; qLerp.Lerp (q0, q1, t); qLerp.Normalize(); // linear quat interpolation
   sVec3 axis; float angle;
   qLerp.ToAxisAndAngle (axis, angle);
   sVec3 proofQ = axis * angle;




At the end proofQ does not exactly match proofV - the difference is too big to say it's a precission issue (max error around 0.01 - same with double prec.).
I'm confused... maybe off-topic, but interesting.

May some math guru enlighten us