# c++ How to use angular velocity that derived from inertia and force(torque) in 3d

• I am relatively new to game development. May my terminology and description are not appropriate. Please excuse my poor phrasing and help me by giving advice on how to question better if this question seems less fitting. I really appreciate your efforts.

Hi. I am having hard time interpreting the set of values I have. I have inertia and force(torque) in terms of x y z. FYI I used x and y coordinates as my ground, flat coordinates and z as my up/down.

I am assuming that since f = ma, that angular acceleration must be a = f / m. So I divide my torque by inertia. Then I add those x y z values to my angular velocity variable's x y z.

However these x y z values confuse me. Don't I need angle/sec or radian/sec sort of values in order to apply rotation? The x y z values I have seemed to not say anything about radians or angular movement.

• Question : If I have ( 1, 2, 3 ) or any ( x, y, z ) as my angular velocity, how do I actually apply it as angular movement?

FYI Here I am pasting my code :

float mass = 100;
float devidedMass = 1.0/12 * mass;
Vec3 innertia(
devidedMass* (_box._size.z*_box._size.z + _box._size.x*_box._size.x),
devidedMass* (_box._size.y*_box._size.y + _box._size.x*_box._size.x),
devidedMass* (_box._size.y*_box._size.y + _box._size.z*_box._size.z ));
box._angAccel += forceAng/innertia;
box._angVelo += box._angAccel;
box._angAccel.allZero();


source of my inertia calculation http://www.health.uottawa.ca/biomech/courses/apa4311/solids.pdf

Trying to tackle solid mechanics are we?

First let me clear up a bit of your confusion regarding various concepts of Newtonian mechanics:

I have inertia and force(torque) in terms of x y z.

Inertia is not a tuple (x,y,z), it is in fact a 3x3 matrix but in the case of a simple box, calculating said matrix within the set of coordinates naturally associated with the box yields a diagonal matrix (i.e. whose only non-zero entries are on the diagonal,so there are 3).

The pdf you linked to gives Ix, Iy and Iz which are the three non-zero entries of the inertia matrix I (but I'm going to call it J because I is often used for the identity matrix). But if J had been calculated relative to a randomly chosen set of coordinates (if the axes in the pictures to the right of the table in the pdf hadn't been in these positions priviledging the geometry of the solid shapes in question) then the matrix J would have been much more complex and, in addition to Ix, Iy and Iz there would have been three other real numbers necessary to define inertia: Ixy = Iyx, Ixz = Izx and Iyz = Izy. Three additional numbers only and not six because J is symmetric.

I don't know what Vec3 innertia(..., ..., ...)); was supposed to mean, but if inertia is going to be a vector, then you'd better make sure your solids always have nice symmetries and invariances with respect to the axes of your axis system.

Classical mechanics tells us that J.(dω/dt) = T where T is the torque, ω is the angular velocity and dt is a small amount of time (your timestep). Both are vectors. To get only , i.e. the small amount of angular velocity that you will add to your current angular velocity, you need to multiply both sides of the equation by the inverse of J and then by dt:

box._angVelo += dt * inverse(J) * T;


If you are sure that J will always be a diagonal matrix then you can avoid calculating a matrix inverse, store the three real numbers making up J in a vector instead, like you did, and use a weird component-wise "vector division" which I guess is what you originally wanted to do with forceAng/innertia, and it will work. But it will harm the realism of the physics of many solids that aren't shaped in a an axes-friendly way.

Next:

Don't I need angle/sec or radian/sec sort of values in order to apply rotation? The x y z values I have seemed to not say anything about radians or angular movement.

Angular velocity is, conceptually, the first derivative of some kind of "orientation", that is to say, it is the amount by which this orientation varies with respect to time. The angular velocity's components are indeed in radians per second. Ideally, you'd like to be able to multiply the angular velocity of a solid with a timestep dt and get a small variation of some kind of "orientation vector" which you could then add to the current orientation vector, pretty much like you did with velocity and position. But it isn't that easy.

What we call an "orientation" is in fact a rotation operation. Asking "How is my solid orientated?" in fact amounts to asking "What rotation should I apply to my solid so that its orientation is the one I want?".

There are three ways to encode rotations:

• A rotation can be defined as a vector, or a tuple of three real numbers representing axis and angle, but only at the price of relatively heavy computation can we use this formalism to compute a rotated vector.

• Another formalism used to encode rotations is that of quaternions, but while still better than the axis/angle representation complexity-wise, they're not the best at computing rotated vectors. In fact, they're the best solution when you must chain rotations, apply several in sequence.

• Finally, the solution we're going to use here: good ol' matrices. Multiplying a rotation matrix with the vector you want to rotate is the fastest way to rotate said vector (although chaining matrix products is less efficient than chaining quaternion products, and a 3x3 matrix is 9 real numbers as opposed to 4 for a quaternion and 3 in the angle/axis representation).

That being said, let's answer your

Question : If I have ( 1, 2, 3 ) or any ( x, y, z ) as my angular velocity, how do I actually apply it as angular movement?

We want to change the matrix encoding our solid's orientation based on the angular velocity we got from the equations above.

1. But how do we even use that matrix to get our solid in position anyway? Well, let's imagine that the solid is at the center of the 3D space, and in its original orientation. If G is its center of mass, then G is exactly at the origin of 3D space. The solid has not ever been rotated. Let's call the solid, when it is in this position, the "base version" of the solid.
2. Now, suppose that our solid has moved since the time it was in its base version. Let's call the rotation matrix that will give us its position R. Since this matrix depends on time, R(t). G has a position vector OG which we have already computed with little effort at this point of the computation (G is now away from the origin O). With all this, at any time t, for any given point P belonging to our solid in its new position,

OP(t) = OG(t) + R(t).OP'

where P' is the point of the base version of the solid corresponding to P. Like, it's the same point, but in the base version. If I have a pointy hat and I move it upside down and onto that shelf over there, if the tip of the hat is P, then P' is the tip of the hat when upright and centered. Making this distinction between the base version of the solid and the real, displaced solid allows us to compute the position of the solid by storing only one known, clean, easy to manipulate version of the solid as well as the transformation that allows us to move it to where we want it to be. You don't store the current position of the solid and apply a small transformation to it. You store a base version of the solid and refresh the transformation that takes it from its base version to its actual version.

This has a ton of other advantages. But i'm not going to cover them, because that's not the question.
3. Now let's modify R(t) to get R(t+dt). R maps the vectors (1,0,0), (0,1,0) and (0,0,1) onto the vectors contained in R's first, second and third column respectively, per the definition of a matrix. Knowing these three column vectors is necessary and sufficient to know R. In fact, from any two of them you can deduce the third. That's a property of rotation matrices. If we can compute the new values of R.(1,0,0) and R.(0,1,0), we'll be set.
4. Let's take the first derivative of the important equation above. From the "positions" equation we deduce the "velocities" equation: the velocity of any point P belonging to our moving solid is given by

Vp = Vg + R'(t).OP'

where R'(t) is the first derivative of R(t). But the laws of kinematics also tell us that

Vp = Vg + ω×GP

where × denotes the cross product. So there you go:

R'(t).OP' = ω×GP

So if OP' is to be (1,0,0) then what is P? It's the point you get by applying the equation of positions to OP' = (1,0,0). Let's now define

R1 = R.(1,0,0),
R2 = R.(0,1,0)
and R3 = R.(0,0,1)

as well as

R'1 = R'.(1,0,0),
R'2 = R'.(0,1,0)
and R'3 = R'.(0,0,1).

We have

OP = OG + R(t).OP'
= OG + R1

and also OP = OG + GP so that

GP = R1.

Therefore,

ω×GP = ω×R1 and since, as we said, OP' = (1,0,0) and R'(t).OP' = ω×GP = R'1, we finally have

R'1 = ω×R1.

The same applies to (0,1,0) so that R'2 = ω×R2. R'3 is obtained from R'1 and R'2 with the relation R'3 = R'1×R'2 but we could also get it using ω. At any rate, we now fully know R'(t).

5. Finally, R(t+dt) = R(t) + R'(t).dt. That's your new "orientation". With the positions equation, you now have the position of the solid at time t+dt.

I have a bit of work to do so I'll come back later and add pseudocode that translates this, with a few necessary remarks regarding the loss of precision due to floating-point operations, etc. In the end, the code is very short.

Edit : I forgot to mention that R is also required to compute ω(t+dt). That is because of the very definition of J. To quote Wikipedia on the moment of inertia,

The use of the inertia matrix in Newton's second law assumes its components are computed relative to axes parallel to the inertial frame and not relative to a body-fixed reference frame. This means that as the body moves the components of the inertia matrix change with time. In contrast, the components of the inertia matrix measured in a body-fixed frame are constant.

To use Newton's second law and compute the angular velocity of our solid, we must either bring our problem back to the solid's own frame (base version of the solid, centered at the origin, using the usual axes of space) and use the one J we know, or we can re-compute J every time the solid rotates. Bringing the problem back within the solid's frame is much easier as it requires only linear transformations (understand rotations) of the variables of the problem (ω, T, etc) whereas computing J again would mean integral calculus over the rotated solid, etc, it would be a mess.

The linear transformations in question are simple: since what separates our inertial frame (relative to the inertial frame, the solid has moved) from the solid's own frame (in the solid's own frame, the center of mass is at the origin, etc) is a rotation (+ a translation but it does not interest us as we are looking at torque and rotation stuff), all you have to do is undo this rotation (apply the inverse of the orientation matrix) to the torque T and to ω(t), do your computations in the solid's frame where J is the constant matrix you know and then rotate everything back:

J.inverse(R).(dω)/(dt) = inverse(R).T

therefore

dω = dt.R.inverse(J).inverse(R).T

i.e.

ω(t+dt) = ω(t) + dt.R.inverse(J).inverse(R).T

but R is an orthogonal matrix because it is a rotation matrix so we have transpose(R) = inverse(R), and a matrix transpose is much easier to compute than an inverse, so finally,

ω(t+dt) = ω(t) + dt.R.inverse(J).transpose(R).T

Notice that the formula for the inverse of the inertia matrix in the inertial frame has appeared: we do have

R.inverse(J).transpose(R) = inverse( R.J.transpose(R) )

Finally, here's the (unbelievably simple) code:

// Notice that you NEVER use J itself, only its inverse, in calculations,
// so you should have a stored copy of inverse(J) as an object attribute
// for your solid for instance. Here I call it "invJ".

// We can use either omega(t) for the next step or omega(t+dt), it's practically
// the same thing. We now want to compute the new R(t+dt).
// We get the three columns of R as an array of Vec3 objects
Vec3[] rColumns = R.getColumns();
// We will store the columns of R(t+dt) in another array
Vec3[] nextrColumns;
// We apply our formulas
for (int i=0; i<3; i++) {

You do have a measurement that is in r/s.
This is because torque = moment of inertia * (delta angular velocity / delta time).
So we can derive delta angular velocity = delta time * (torque / moment of inertia). Since you're probably working with time steps, you can throw out the delta components. The resulting vector is equal to the angular velocity in radians over seconds. Or, in other words, a quaternion.
If you were to write out the vector in terms of variables it would be (change in x angle, change in y angle, change in z angle). They are also unit vectors in this case, so all you need to do is add the quaternion to the object's angle, component-wise.