# Rotation and Inertia Tensors

February 16, 2013

Hello, I’m Glenn Fiedler and welcome to Virtual Go, my project to simulate a Go board and stones.

In the previous article we detected collision between the go stone and the go board.

We’re working up to calculating collision response so the stone bounces and wobbles before coming to rest on the board, but in order to reach this goal we first need to lay some groundwork.

It turns out that irregular objects, like go stones, are easier to rotate about certain axes than others.

This tendency is a property of the distribution of mass in a go stone. In other words, the shape of the go stone contributes to how it moves, and how it reacts to collisions. If we neglect this effect we’ll never reproduce that unique ‘wobble’ go stones have when placed on the board.

It makes sense that we’re going to need to understand this before we apply collision response.

Lets get started.

## Rotation in 3D

Physics in 3D is considerably more difficult, for me at least, than physics in 2D.

Most of this this has to do with the complexity of rotation in three dimensions.

Consider the following case in two dimensions:

It’s easy because there is only one possible axis for rotation. Clockwise or counter-clockwise, no other rotation is possible.

It follows that we can represent the orientation of an object in 2D with a single theta value, angular velocity with a scalar radians per-second, and a scalar ‘moment of inertia’ that works just like an angular equivalent of mass, eg: how hard it is to rotate the object.

When we move to three dimensions suddenly rotation can occur about any axis. Orientation becomes a quaternion, angular velocity is a vector, and now for irregular shaped objects, like go stones, we need some way to indicate that certain axes of rotation are easier to rotate about than others.

How can we represent an angular mass that depends on the shape of the object and the axis of rotation?

## Inertia Tensor

The solution is to use an inertia tensor.

An inertia tensor is a 3×3 matrix with different rules to a normal matrix. It rotates and translates differently, but otherwise, it behaves like a 3×3 matrix and you use it to transform vectors.

Specifically, you use the inertia tensor to transform angular velocity to angular momentum, and the inverse of the inertia tensor to transform angular momentum to angular velocity.

Now this becomes quite interesting because Newton’s laws guarantee that in a perfectly elastic collision angular momentum is conserved but angular velocity is not necessarily.

Why is this? It’s because angular velocity depends on the axis of rotation, so even if the angular momentum has exactly the same magnitude post-collision the angular velocity can be different if the axis of rotation changes (it usually does).

Because of this we’ll switch to angular momentum as the primary quantity and we’ll derive angular velocity from it. For consistency we’ll also switch from linear velocity to linear momentum:

    const float gravity = 9.8f * 10;
const float fps = 60.0f;
const float dt = 1 / fps;

while ( !quit )
{
stone.rigidBody.linearMomentum += vec3f(0,-gravity,0) * stone.rigidBody.mass * dt;

stone.rigidBody.linearVelocity = stone.rigidBody.linearMomentum * stone.rigidBody.inverseMass;
stone.rigidBody.angularVelocity = transformVector( inverseInertiaTensor, angularMomentum );

stone.rigidBody.position += stone.rigidBody.velocity * dt;

quat4f spin = AngularVelocityToSpin( stone.rigidBody.orientation,
stone.rigidBody.angularVelocity );

stone.rigidBody.orientation += spin * dt;
stone.rigidBody.orientation = normalize( stone.rigidBody.orientation );

RenderStone( stone );

UpdateDisplay();
}


## Calculating The Inertia Tensor

Now we need a way to calculate the inertia tensor.

The general case is quite complicated because inertia tensors are capable of representing shapes that are non-symmetrical about the axis of rotation.

$I = \begin{bmatrix} I_{xx} & I_{xy} & I_{xz} \\ I_{yx} & I_{yy} & I_{yz} \\ I_{zx} & I_{zy} & I_{zz} \end{bmatrix}$

For example, think of an oddly shaped object attached to a drill bit off-center and wobbling about crazily as the drill spins.

The good news is that due to the symmetry of the go stone and because we’re always rotating about the center of mass, our inertia tensor is much simpler:

$I = \begin{bmatrix} I_{x} & 0 & 0 \\ 0 & I_{y} & 0 \\ 0 & 0 & I_{z} \end{bmatrix}$

All we need to do in our case is to determine the Ix, Iy and Iz values.

These values represent how difficult it is to rotate the go stone about the x,y and z axes.

Interestingly, due to symmetry of the go stone, all axes on the xz plane are identical.

So really, we only need to calculate Ix and Iy because Iz = Ix.

## Numerical Integration

First lets calculate the inertia tensor via numerical integration.

To do this we just need to know is how difficult it is rotate a point about an axis.

Once we know this we can approximate the moment of inertia of a go stone by breaking it up into a discrete number of points and summing up the moments of inertia of all these points.

It turns out that the difficulty of rotating a point mass about an axis is proportional to the square of the distance of that point from the axis and the mass of the point. $I = mr^2$

This is quite interesting because it indicates that the distribution of mass has a significant effect on how difficult it is to rotate an object about an axis.

One consequence of this is that a hollow pipe is actually more difficult to rotate than a solid pipe of the same mass. Of course, this is not something we deal with in real life often, because a solid pipe of the same material would be much heavier, and therefore harder to rotate due to increased mass, but if you could find a second material of lower density such that the solid pipe was exactly the same mass as the hollow pipe, you would be able to observe this effect.

In our case we know the go stone is solid not hollow, and we can go one step further and assume that the go stone has completely uniform density throughout. This means if we know the mass of the go stone we can divide it by the volume of the go stone to find its density. Then we can divide space around the go stone into a grid, and using this density we can assign a mass to each point in the grid proportional to the density of the go stone.

Now integration is just a triple for loop summing up the moments of inertia for points that are inside the go stone:

inline bool PointInsideBiconvex( vec3f point,
const Biconvex & biconvex,
float epsilon = 0.001f )
{
const float sphereOffset = biconvex.GetSphereOffset();

if ( point.y() >= 0 )
{
// check bottom sphere (top half of biconvex)
vec3f bottomSphereCenter( 0, -sphereOffset, 0 );
const float distanceSquared = length_squared( point - bottomSphereCenter );
if ( distanceSquared <= radiusSquared )
return true;
}
else
{
// check top sphere (bottom half of biconvex)
vec3f topSphereCenter( 0, sphereOffset, 0 );
const float distanceSquared = length_squared( point - topSphereCenter );
if ( distanceSquared <= radiusSquared )
return true;
}

return false;
}

float CalculateBiconvexVolume( const Biconvex & biconvex )
{
const float h = r - biconvex.GetHeight() / 2;
return h*h + ( pi * r / 4 + pi * h / 24 );
}

void CalculateBiconvexInertiaTensor( float mass, const Biconvex & biconvex,
float & ix, float & iy, float & iz,
float resolution = 0.1f )
{
ix = 0;
iy = 0;
iz = 0;
const float width = biconvex.GetWidth();
const float height = biconvex.GetHeight();
const float xz_steps = ceil( width / resolution );
const float y_steps = ceil( height / resolution );
const float dx = width / xz_steps;
const float dy = height / y_steps;
const float dz = width / xz_steps;
const float sx = -width / 2;
const float sy = -height / 2;
const float sz = -width / 2;
const float v = CalculateBiconvexVolume( biconvex );
const float p = mass / v;
const float m = dx*dy*dz * p;
for ( int index_z = 0; index_z <= xz_steps; ++index_z )
{
for ( int index_y = 0; index_y <= y_steps; ++index_y )
{
for ( int index_x = 0; index_x <= xz_steps; ++index_x )
{
const float x = sx + index_x * dx;
const float y = sy + index_y * dy;
const float z = sz + index_z * dz;

vec3f point(x,y,z);

if ( PointInsideBiconvex( point, biconvex ) )
{
const float rx2 = z*z + y*y;
const float ry2 = x*x + z*z;
const float rz2 = x*x + y*y;

ix += rx2 * m;
iy += ry2 * m;
iz += rz2 * m;
}
}
}
}
}


## Interpreting The Inertia Tensor

A size 33 japanese go stone has width 22mm and height 9.2mm:

Using our code to calculate its inertia tensor gives the following result:

$I = \begin{bmatrix} 0.177721 & 0 & 0 \\ 0 & 0.304776 & 0 \\ 0 & 0 & 0.177721 \end{bmatrix}$

As expected, Ix = Iz due to the symmetry of the go stone.

The inertia tensor indicates that its much harder to rotate the go stone about the y axis than axes on the xz plane. Why is this?

You can see looking top-down at the go stone when rotating about the y axis a ring of mass around the edge of the stone is multiplied by a large r2 and is therefore difficult to rotate.

Contrast this with the rotation about the x axis, which has a much smaller portion of mass far away from the axis:

As you can see the distribution of mass around the axis tends to dominate the inertia tensor due to the r2 term. The same mass, twice the distance from the axis, is four times more difficult to rotate!

## Closed Form Solution

Exact equations are known for the moments of inertia of many common objects.

With a bit of math we can calculate closed form solutions for the moments of inertia of a go stone.

To calculate the exact equation for Iy we start with the moment of inertia for a solid disc: $I = 1/2mr^2$

Then we can integrate again, effectively summing up the moments of inertia of an infinite number of thin discs making up the top half of the go stone.

This leads to the following integral: $\int_0^{h/2} (r^2-(y+r-h/2)^2)^2\,dy$

With a little bit of help from Wolfram Alpha we get the following result:

float CalculateExactIy( const Biconvex & biconvex )
{
// http://wolframalpha.com
// integrate ( r^2 - ( y + r - h/2 ) ^ 2 ) ^ 2 dy from y = 0 to h/2
//  => 1/480 h^3 (3 h^2-30 h r+80 r^2)
const float h = height;
const float h2 = h * h;
const float h3 = h2 * h;
const float h4 = h3 * h;
const float h5 = h4 * h;
const float r2 = r * r;
const float r3 = r2 * r;
const float r4 = r3 * r;
return = pi * p * ( 1/480.0f * h3 * ( 3*h2 - 30*h*r + 80*r2 ) );
}


Plugging in the values for a size 33 stone, we get 0.303588 which is close to the approximate solution 0.304776.

Verifying exact solutions against numeric ones is a fantastic way to check your calculations.

Can you derive the equation for Ix?

# Next:Collision Response And Coulomb Friction

Donations offset hosting costs and encourage me to write more articles!

Hiroki Mori March 22, 2013 at 3:09 pm

Hi Glenn,
This had been what I wanted to do for long time… you just realized my dream
Keep up the good work.

Hiroki Mori
The author of “The Interactive Way to Go”

Glenn Fiedler March 22, 2013 at 3:21 pm

Thank you!