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

In previous articles we have mathematically defined the shape of a go stone and tessellated this shape so it can be drawn with 3D graphics hardware.

Now we want to make the go stone move.

I want the stone to have mass and obey Newton’s laws of motion so that the simulation is physically accurate. The stone should be accelerated by gravity and fall downwards. I also want the stone to rotate so it tumbles realistically as it falls through the air.

Now lets see how we can simulate this motion with a virtual go stone. __Lets go!__

## The Rigid Body Assumption

Try biting down on a go stone some time and you’ll agree: go stones are very, very hard.

Golf balls are pretty hard too, but if you look at a golf ball being hit by a club in super-slow motion, you’ll see that it deforms considerably during impact.

The same thing happens to all objects in the real world to some degree. Nothing is truly rigid. No real material is so hard that it never deforms.

But this is not the real world. It’s a simulation and here we are free to make whatever assumptions we want. And the smartest simplification we can make at this point is to assume that the go stone is perfectly rigid and does not deform under any circumstance.

This is known as the rigid body assumption.

## Working in Three Dimensions

Because the go stones are rigid all we need to represent the current position of the go stone is the position of its center point P.

As the center point moves so does the rest of the stone.

We’ll represent this position using a three dimensional vector **P**.

Lets define the axes so we know what the x,y,z components of P mean:

- Positive x is to the right
- Positive y is up
- Positive z is into the screen

This is what is known as a left-handed coordinate system. So called because I can use the fingers on my left hand to point out each positive axis direction without breaking my fingers.

I’ve chosen this based purely on personal preference. Also, I’m left-handed and I like my fingers.

## Linear Motion

Now we want to make the stone move.

To do this we need the concept of velocity. Velocity is also a vector but it’s not a point like P. In this case you can think of it more like a direction and a length. The direction of the velocity vector is the direction the stone is moving and the length is the speed that it is moving in some unit per-second. Here I’ll use centimeters per-second because go stones are small.

For example, if we the stone to move to the right at a rate of 5 centimeters per-second then the velocity vector is (5,0,0).

To make the stone move all we have to do is add this velocity to the stone position once per-second.

Now this isn’t very exciting obviously. We’d like to move the stone more smoothly. So instead of updating once per-second, lets update 60 frames per-second (fps).

How much should the stone move each frame? Well, we want to move the same amount of distance in 60 small steps instead of one big one. So if we take 60 steps, but each step is 1/60th the amount of movement, the stone moves at the same overall rate.

You can generalize this to any framerate with the concept of delta time or “dt”. To calculate delta time simply invert the frames per second: dt = 1/fps and you have the amount of time per-frame in seconds.

Then, simply multiply velocity by delta time and you have the change in position per-frame. Add this to position each frame and the stone moves according to the current velocity, independent of framerate.

const float fps = 60.0f; const float dt = 1 / fps; while ( !quit ) { stone.rigidBody.position += stone.rigidBody.velocity * dt; RenderStone( stone ); UpdateDisplay(); }

This is actually a type of numerical integration. Specifically, we are integrating velocity to find the change in position over time.

## Gravitational Acceleration

Next we want to add gravity.

To do this we need to change velocity each frame by some amount downwards due to gravity. Change in velocity is known as acceleration. Gravity provides a constant acceleration of 9.8 meters per-second, per-second, or 98 centimeters per-second, per-second. Acceleration due to gravity is also a vector, and since gravity pulls objects down the acceleration vector is (0,-98,0).

So how much does gravity accelerate the go stone in 1/60th of a second? Well, 98 * 1/60 = 1.633… Hey wait. This is exactly what we did with velocity to get position!

Yes it is. It’s exactly the same. Acceleration integrates to velocity just like velocity integrates to position. And both are multiplied by dt to find the amount to add per-frame, where dt = 1/fps.

float gravity = 9.8f * 10; float fps = 60.0f; float dt = 1 / fps; while ( !quit ) { stone.rigidBody.velocity += vec3f( 0, -gravity, 0 ) * dt; stone.rigidBody.position += stone.rigidBody.velocity * dt; RenderStone( stone ); UpdateDisplay(); }

Now that we’ve added acceleration due to gravity the go stone moves in a parabola just like it does in the real world.

## Angular Motion

Now lets make the stone rotate!

First we have to define how we represent the orientation of the stone. For this we’ll use a quaternion.

Next we need the angular equivalent of velocity known as… wait for it… angular velocity. It too is a vector aka a direction and a length. The direction is the axis of rotation and the length is the rate of rotation in radians per-second. One full rotation is 2*pi radians or 360 degrees so if the length of the angular velocity vector is 2*pi the object rotates around the axis once per-second.

Because we are using a left handed coordinate system the direction of rotation is clockwise about the positive axis. You can remember this by sticking your thumb of your left hand in the direction of the axis of rotation and curling your fingers. The direction your fingers curl is the direction of rotation. Notice that if you do the same thing with your right hand the rotation is the other way.

How do we integrate orientation from angular velocity? Orientation is a quaternion and angular velocity is a vector. We can’t just add them together.

The solution requires a reasonably solid understanding of quaternion math and how it relates to complex numbers. Long story short we need to convert our angular velocity into a quaternion form and then we can integrate that just like we integrate any other vector. For a full derivation of this result please refer to this excellent article.

Here is the code that I use to convert angular velocity into quaternion form:

inline quat4f AngularVelocityToSpin( quat4f orientation, vec3f angularVelocity ) { const float x = angularVelocity.x(); const float y = angularVelocity.y(); const float z = angularVelocity.z(); return 0.5f * quat4f( 0, x, y, z ) * orientation; }

Once I have this spin quaternion, I can integrate it to find the change in the orientation quaternion just like any other vector.

const float fps = 60.0f; const float dt = 1 / fps; while ( !quit ) { quat4f spin = AngularVelocityToSpin( stone.rigidBody.orientation, stone.rigidBody.angularVelocity ); stone.rigidBody.orientation += spin * iteration_dt; stone.rigidBody.orientation = normalize( stone.rigidBody.orientation ); RenderStone( stone ); UpdateDisplay(); }

The only difference is that after integration I renormalize the quaternion to ensure that it doesn’t drift from unit length. This is important because quaternions that are not unit length do not represent pure rotation.

## Why Quaternions?

Which brings us to an important point. Why are we using quaternions to represent orientation and not matrices?

Here are several good reasons:

- It’s much easier to integrate angular velocity using a quaternion than a 3×3 matrix
- It’s much faster to normalize a quaternion than it is to orthonormalize a 3×3 matrix
- It’s really easy to interpolate between two quaternions and this is very useful

We’ll still use matrices but as a secondary quantity. This means that each frame after we integrate we convert the quaternion into a 3×3 rotation matrix and combine it with the position into a 4×4 rigid body matrix and its inverse like this:

mat4f RigidBodyMatrix( vec3f position, quat4f rotation ) { mat4f matrix; rotation.toMatrix( matrix ); matrix.value.w = simd4f_create( position.x(), position.y(), position.z(), 1 ); return matrix; } mat4f RigidBodyInverse( const mat4f & matrix ) { mat4f inverse = matrix; vec4f translation = matrix.value.w; inverse.value.w = simd4f_create(0,0,0,1); simd4x4f_transpose_inplace( &inverse.value ); vec4f x = matrix.value.x; vec4f y = matrix.value.y; vec4f z = matrix.value.z; inverse.value.w = simd4f_create( -dot( x, translation ), -dot( y, translation ), -dot( z, translation ), 1.0f ); return inverse; }

Now whenever we transform vectors want to go go in/out of go stone body space we’ll use this matrix and its inverse. It’s the best of both worlds.

## Bringing It All Together

The best thing about rigid body motion is that you can calculate linear and angular motion separately and combine them together and it just works.

Here is the final code:

const float gravity = 9.8f * 10; const float fps = 60.0f; const float dt = 1 / fps; while ( !quit ) { stone.rigidBody.velocity += vec3f( 0, -gravity, 0 ) * dt; 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(); }

And here is the end result:

http://youtu.be/gmZGktbTAGA&w=704

**Next:** Collision Detection: Stone vs. Board

If you enjoyed this article please consider making a small donation. __Donations encourage me to write more articles!__

IIUC quat4f( float w, float x, float y, float z ) is

cos(w / 2) + (x * i + y * j + z * k) * sin(w / 2)

cos(0 / 2) = 1;

sin(0 / 2) = 0;

so quat4f(0, x, y, z) is (1, 0, 0, 0) for all values of angular velocity, no?

What am I missing?

That the quat4f(0,x,y,z) is immediately multiplied by the orientation quaternion rather than being used directly. Read the linked article for the derivation of this result.

I realized that quat4f(0, x, y, z) shouldn’t be interpreted as a quaternion that rotates by angle 0 around axis x,y,z.

It’s a quaternion q with values:

q.w = 0;

q.x = x;

q.y = y;

q.z = z;

Yep.

Thank you for the article!

A small note: I think AngularVelocityToSpin() function has a bug.

It should return 0.5f * quat4f( length(angularVelocity ), normalize(angularVelocity) ) * orientation

quat4(0, x, y, z) would be always constant quaternion and doesn’t reflect the angular velocity

That’s not correct. Note that the quat4f ctor is quat4f( float w, float x, float y, float z ). cheers