**Introduction**

Hi, I’m Glenn Fiedler and welcome to the third article in my series on Game Physics.

In the previous article we discussed how to integrate our physics simulation forward at fixed delta time increments, regardless of display framerate.

Now are going to simulate objects that move linearly and rotate in three dimensions.

We will concentrate on a type of object called a rigid body. Rigid bodies cannot bend, compress or deform in any way. This makes their motion much easier to calculate.

To simulate the motion of rigid bodies, we must study both rigid body kinematics and rigid body dynamics. Kinematics is the study of how an object moves in the absence of forces, while dynamics describes how an object reacts to them. Together they provide all the information you need to simulate the motion of a rigid body in three dimensions.

Along the way I will show you how to integrate vector quantities, handle rotations in three dimensions and integrate to find the motion of your object as it moves and spins around the world.

Enjoy!

**Moving in the Third Dimension**

As long as we only have single floating point values for position and velocity our physics simulation is limited to motion in a single dimension, and a point moving from side to side on the screen is pretty boring!

We want our object to be able to move in three dimensions: left and right, forward and back, up and down. If we apply the equations of motion to each dimension separately, we can integrate each dimension in turn to find the motion of the object in three dimensions.

Or we could just use vectors.

Vectors are a mathematical type representing an array of numbers. A three dimensional vector has three components x, y and z. Each component corresponds to a dimension. In this article x is left and right, y is up and down, and z is forward and back.

In C++ we will implement vectors using a struct as follows:

struct Vector { float x,y,z; };

Addition of two vectors is defined as adding each component together. Multiplying a vector by a floating point number is the same as just multiplying each component. Lets add overloaded operators to the vector struct so that we can perform these operations in code as if vectors are a native type:

struct Vector { float x,y,z; Vector operator + ( const Vector &other ) { Vector result; result.x = x + other.x; result.y = y + other.y; result.z = z + other.z; return result; } Vector operator*( float scalar ) { Vector result; result.x = x * scalar; result.y = y * scalar; result.z = z * scalar; return result; } };

Now instead of maintaining completely seperate equations of motion and integrating seperately for x, y and z, we convert our position, velocity, acceleration and force to vector quantities, then integrate the vectors directly using the equations of motion from the first article:

F= madv/dt =adx/dt =v

Notice how **F**, **a**, **v** and **x** are written in bold. This is the convention used to distinguish vector quantities from single value (scalar) quantities such as mass m and time t.

Now that we have the equations of motion in vector form, how do we integrate them? The answer is exactly the same as we integrated single values. This is because we have already added overloaded operators for adding two vectors together, and multiplying a vector by a scalar, and this is all we need to be able to drop in vectors in place of floats and have everything just work.

For example, here is a simple Euler integration for vector position from velocity:

position = position + velocity*dt;

Notice how the overloaded operators make it look exactly the same as an Euler integration for a single value. But what is it really doing? Lets take a look at how we would implement vector integration without the overloaded operators:

position.x = position.x + velocity.x * dt; position.y = position.y + velocity.y * dt; position.z = position.z + velocity.z * dt;

As you can see, its exactly the same as if we integrated each component of the vector separately! This is the cool thing about vectors. Whether we integrate vectors directly, or integrate each component separately, we are doing exactly the same thing.

**Structuring for RK4**

In the example programs from previous articles we drove the simulation from acceleration assuming unit mass. This kept the code nice and simple, but from now on every object will have its own mass in kilograms so the simulation needs be driven by forces instead.

There are two ways we can do this. First, we can divide force by mass to get acceleration, then integrate this acceleration to get the velocity, and integrate velocity to get position.

The second way is to integrate force directly to get momentum, then convert this momentum to velocity by dividing it by mass, then finally integrate velocity to get position. Remember that momentum is just velocity multiplied by mass:

dp/dt =Fv=p/m dx/dt =v

Both methods work, but the second way is more consistent with the way that we must approach rotation later in the article, so we’ll use that.

When we switch to the second technique, it adds a new level of complexity to our code. Each time that momentum changes we need to make sure that the velocity is recalculated by dividing momentum by mass. Doing this manually everywhere that momentum is changed would be error prone and a chore. So we now separate all our state quantities into primary, secondary and constant values, and add a method called ‘recalculate’ to the State struct which is responsible for updating all the secondary values from the primary ones:

struct State { // primary Vector position; Vector momentum; // secondary Vector velocity; // constant float mass; float inverseMass; void recalculate() { velocity = momentum * inverseMass; } }; struct Derivative { Vector velocity; Vector force; };

If we make sure that recalculate is called whenever any of the primary values change, then our secondary values will always stay in sync. This may seem like overkill just to handle converting momentum to velocity, but as our simulation becomes more complex we will have many more secondary values, so it is important to design a system that scales well.

**Spinning around**

So far we have covered linear motion, we can simulate an object so that it moves in 3D space, but it cannot rotate.

The good news is that rotational equivalents to force, momentum, velocity, position and mass exist, and once we understand how they work, integration of rotational physics state can be performed using our RK4 integrator.

So lets start off by talking about how an object rotates. The bodies that we are simulating are rigid meaning that they cannot deform. This might not sound important, but its actually the key simplification in rigid bodies! Basically it means that we can treat the linear and rotational parts of an object’s motion as being entirely separate: a linear component (position, velocity, momentum, mass) and a rotational component rotating about the center of mass.

The center of mass of the object is the weighted average of all points making up the body by their mass. For objects of uniform density, the center of mass is always the geometric center of the object, for example the center of a sphere or a cube.

So how do we represent how the object is rotating? If you think about it a bit, you’ll realize that for a rigid body rotation can only ever be around a single axis, so the first thing we need to know is what that axis is. We can represent this axis with a unit length vector. Next we need to know how fast the object is rotating about this axis in radians per second.

If we know the center of mass of the object, the axis of rotation, and the speed of rotation then we have the all the information we need to describe how an object is rotating.

The standard way of representing rotation over time is by combining the axis and the speed of rotation into a single vector called angular velocity. The length of the angular velocity vector is the speed of rotation in radians while the direction of the vector indicates the axis of rotation. For example, an angular velocity of (2Pi,0,0) indicates a rotation about the x axis doing one revolution per second.

But what direction is this rotation in? In the example source code I use a right handed coordinate system which is standard when using OpenGL. To find the direction of rotation just take your right hand and point your thumb down the axis – your fingers will now curl in the direction of rotation. If your 3D engine uses a left handed coordinate system then just use your left hand instead.

Why do we combine the axis and rate of rotation into a single vector? Doing so gives us a single vector quantity that is easy to manipulate just like velocity for linear motion. We can easily add and subtract changes to angular velocity to change how the object is rotating just like we can add and subtract from linear velocity. If we stuck with a unit length vector and scalar for rotation speed then it would be much more complicated to apply these changes.

But there is one very important difference between linear and angular velocity. Unlike linear velocity, there is no guarantee that angular velocity will remain constant over time in the absence of forces. In other words, angular momentum is conserved while angular velocity is not. This means that we cannot trust angular velocity as a primary value and we need to use angular momentum instead.

**Angular momentum, inertia and torque**

Just as velocity and momentum are related by mass in linear motion, angular velocity and angular momentum are related by a quantity called the rotational inertia. This tensor is a measurement of how much effort it takes to spin an object around an axis. It depends on both the shape of the object and how much it weighs.

In the general case, rotational inertia is represented by a 3×3 matrix called an inertia tensor. Here we make a simplifying assumption by discussing physics in the context of simulating a cube. Because of the symmetries of the cube, we only need a single value for the rotational inertia: 1/6 x size^2 x mass, where size is the length of the sides of the cube.

Just as we integrate linear momentum from force, we integrate angular momentum directly from the rotational equivalent of force called torque. You can think of torque just like a force, except that when it is applied it induces a rotation around an axis in the direction of torque vector rather than accelerating the object linearly. For example, a torque of (1,0,0) would cause a stationary object to start rotating about the x axis.

Once we have angular momentum integrated, we multiply it by the inverse of the rotational inertia to get the angular velocity, and using this angular velocity we integrate to get the rotational equivalent of position called orientation.

However, as we will see, integrating orientation from angular velocity is a bit more complicated!

**Orientation in 3D**

This complexity is due to the difficulty of representing orientations in three dimensions.

In two dimensions orientations are easy, you just keep track of an angle in radians and you are done. In three dimensions it becomes much more complex. It turns out that you must either use 3×3 rotation matrices or quaternions to correctly represent the orientation of an object.

For reasons of simplicity and efficiency I’m going to use quaternions to represent the orientation instead of matrices. This also gives us an easy way to interpolate between the previous and current physics orientation to get smooth framerate independent animation as per the time stepping scheme outlined in the previous article.

Now there are plenty of resources on the internet which explain what quaternions are and how unit length quaternions are used to represent rotations in three dimensions. Here is a particularly nice one. What you need to know however is that, effectively, unit quaternions represent an axis of rotation and an amount of rotation about that axis. This may seem similar to our angular velocity, but quaternions are four dimensional vectors instead of three, so mathematically they are actually quite different!

We will represent quaternions in code as another struct:

struct Quaternion { float w,x,y,z; };

If we define the rotation of a quaternion as being relative to an initial orientation of the object (what we will later call body coordinates) then we can use this quaternion to represent the orientation of the object at any point in time. Now that we have decided on the representation to use for orientation, we need to integrate it over time so that the object rotates according to the angular velocity.

**Integrating orientation**

We are now presented with a problem. Orientation is a quaternion but angular velocity is a vector. How can we integrate orientation from angular velocity when the two quantities are in different mathematical forms?

The solution is to convert angular velocity into a quaternion form, then to use this quaternion to integrate orientation. For lack of a better term I will call this time derivative of orientation “spin”. Exactly how to calculate this spin quaternion is described in detail here.

Here is the final result:

dq/dt = spin = 0.5wq

Where **q** is the current orientation quaternion, and **w** is the current angular velocity in quaternion form (0,x,y,z) such that x, y, z are the components of the angular velocity vector. Note that the multiplication done between **w** and **q** is quaternion multiplication.

To implement this in code we add spin as a new secondary quantity calculated from angular velocity in the recalculate method. We also add spin to the derivatives struct as it is the derivative of orientation:

struct State { // primary Quaternion orientation; Vector angularMomentum; // secondary Quaternion spin; Vector angularVelocity; // constant float inertia; float inverseInertia; void recalculate() { angularVelocity = angularMomentum * inverseInertia; orientation.normalize(); Quaternion q( 0, angularVelocity.x, angularVelocity.y, angularVelocity.z ) spin = 0.5f * q * orientation; } }; struct Derivatives { Quaternion spin; Vector torque; };

Integrating a quaternion, just like integrating a vector, is as simple as doing the integration for each value separately. The only difference is that after integrating orientation we must renormalize the orientation quaternion to make it unit length, to ensure that it still represents a rotation.

This is required because errors in integration accumulate over time and make the quaternion ‘drift’ away from being unit length. I like to do the renormalization in the recalculate method for simplicity, but you can get away with doing it less frequently if cpu cycles are tight.

Now in order to drive the rotation of the object, we need a method that can calculate the torque applied given the current rotational state and time just like the force method we use when integrating linear motion. eg:

Vector torque( const State &state, float t ) { return Vector(1,0,0) - state.angularVelocity * 0.1f; }

This function returns an acceleration torque to induce a spin around the x axis, but also applies a damping over time so that at a certain speed the accelerating and damping will cancel each other out. This is done so that the rotation will reach a certain rate and stay constant instead of getting faster and faster over time.

**Combining linear and angular motion**

Now that we are able to integrate linear and rotational effects, how can they be combined into one simulation? The answer is to just integrate the linear and rotational physics state separately and everything works out. This is because the objects we are simulating are rigid so we can decompose their motion into separate linear and rotational components. As far as integration is concerned, you can treat linear and angular effects as being completely independent of each other.

Now that we have an object that is translating and rotating through three dimensional space, we need a way to keep track of where it is. We must now introduce the concepts of body coordinates and world coordinates.

Think of body coordinates in terms of the object in a convenient layout, for example its center of mass would be at the origin (0,0,0) and it would be oriented in the simplest way possible. In the case of the simulation that accompanies this article, in body space the cube is oriented so that it lines up with the x, y and z axes and the center of the cube is at the origin.

The important thing to understand is that the object remains stationary in body space, and is transformed into world space using a combination of translation and rotation operations which put it in the correct position and orientation for rendering. When you see the cube animating on screen it is because it is being drawn in world space using the body to world transformation.

We have the raw materials to implement this transform from body coordinates into world coordinates in the position vector and the orientation quaternion. The trick to combining the two is to convert each of them into 4×4 matrix form which is capable of representing both rotation and translation. Then we combine the two transformations into a single matrix by multiplication. This combined matrix has the effect of first rotating the cube around the origin to get the correct orientation, then translating the cube to the correct position in world space. See this article for details on how this is done.

If we then invert this matrix we get one that has the opposite effect, it transforms points in world coordinates into the body coordinates of the object. Once we have both these matrices we have the ability to convert points from body to world coordinates and back again which is very handy. These two matrices become new secondary values calculated in the ‘recalculate’ method from the orientation quaternion and position vector.

**Forces and torques**

We can apply separate forces and torques to an object individually, but we know from real life that if we push an object it usually makes it both move and rotate. So how can we break down a force applied at a point on the object into a linear force which causes a change in momentum, and a torque which changes angular momentum?

Given that our object is a rigid body, what actually happens here is that the entire force applied at the point is applied linearly, plus a torque is also generated based on the cross product of the force vector and the point on the object relative to the center of mass of the object:

F_{linear}=FF_{torque}=Fx (p-x)

Where **F** is the force being applied at point **p** in world coordinates, and **x** is the center of mass of the object.

This seems counterintuitive at first. Why is the force being applied twice? Once to linear and once to rotational motion?

What is happening here is our everyday experience with objects clouding the true behavior of an object under ideal conditions.

Remember your pushbike when you were a kid? You would have to change your tire and flip the bike upside down. You could spin the tire around by pushing on it. You don’t see any linear motion here, just rotation, so what is going on? The answer of course is that the axle of the wheel is counteracting the linear component of the force you applied, leaving only the rotational component. Not convinced? Imagine what would happen if you tried to ride your bike without an axle in your wheel…

Another example: consider a bowling ball lying on a slippery surface such as ice so that no significant friction is present. Now in your mind try to work out a way that you can apply a force at a single point on the surface of the bowling ball such that it will stay completely still while rotating on the spot. There is no way you can do this! Any point where you push would also make the bowling ball move linearly as well as rotate. To apply a pure rotation you’d have to push on both sides of the ball, canceling the linear component of your force out leaving only torque.

So remember, whenever you apply a force to an object there will always be a linear force component which causes the object to accelerate linearly, as well as, depending on the direction of the force, a rotational component that causes the object to rotate.

**Velocity at a point**

The final piece of the puzzle is how to calculate the velocity of a single point in the rigid body. To do this we start with the linear velocity of the object, because all points must move with this velocity to keep it rigid, then add the velocity at the point due to rotation.

This velocity due to rotation will not be constant for every point in the body if it is rotating, as each point in the body must be spinning around the axis of rotation. Combining the linear and angular velocities, the total velocity of a point in the rigid body is:

v_{point}=v_{linear}+v_{angular}cross (p-x)

Where **p** is the point on the rigid body and **x** is the center of mass of the object.

**Conclusion**

We have covered the techniques required to simulate linear and rotational movement of a rigid body in three dimensions. By combining the linear and rotational physics into a single physics state and integrating, we can simulate the motion of a cube in three dimensions as it moves and spins around.

Click here to download the source code for this article.

**Next:** Spring Physics

If you enjoyed this article please donate. __Donations encourage me to write more articles!__

A great article! Thank you for this detailed insight into physics and it’s implementation in games.

Looking forward to the next article.

Very well done

Hi,

Great article, however it would be good to see some more information regarding inertia tensors. Although this implementation works fine with cubes, it’s a serious limitation. (not many games that use just cubes!)

I can see that extending to include full 3×3 tensors would be simple, but for completeness it would make a great article discussing how to create the tensor for a general 3D body, and how to use it in the code.

i’ll write an article as soon as i actually understand ‘em properly — however, you’ll be interested to know that a lot of commercial engines basically use uniform tensors, irrespective of the geometry, for stability reasons… =)

It’s interesting to know most games are using fairly bad tensors!

I think a nice extension to this that would be good to see would be to generalise this article to handle oblong bodies. this would give a bit more flexability without getting over complicated….

Here are some interesting links which may be of use to people here!

Check out the paper and source for Polyhedral Mass Properties.

http://www.geometrictools.com/LibPhysics/RigidBody/RigidBody.html

This is a system for creating inertia tensors for general 3D triangulated bodies.

for a simple but usefull introduction to what inertia tensors are, try:

http://techhouse.brown.edu/~dmorris/projects/tutorials/inertia.tensor.summary.pdf

Very good,

Actually i got this hit while searching for quaternion rotation.

infact this does not help me that much in that sense,but overall it is good.

Hi there I would like to create a bowling game and wanted to know what the code would be to calculate the position of a bowling ball if all the variables are known.(ie. velocity, mass, rotaion…) any help would be great, thanks.

Hmm…I’ve recreated this code in a gaming language (Blitz3D) and my cube is moving but stuttering pretty badly. This is some great code and I am trying hard to convert it to my language. Let’s hope I can do this.

Bah…too bad there’s no edit feature. I just got this working! Very smooth. I love it. Now to get the cube bouncing around… Coffee!

The article and some comments by others provide a number of great links to to deeper. In particular, the geometrictools.com site is amazing.

Thanks for the article.

Very interesting and useful articles. I found another way, how to model multi rigid body system. http://i31www.ira.uka.de/projekte/mechanik/publications/details2003.php?publicationNumber=1_1

What do you think about it? Look’s a lot simplier then other methods solving constrained dynamics.

Dear Glenn, I’ve got a few questions and hope you can answer some!

1. In the “integration basics” example code you use … evaluate(const State &initial … for both functions, BUT in the “physics in 3d” example code you use … evaluate(State state … and directly modify the state in that function. Doesn’t that falsify the results of further calls (the 3rd and 4th actually for Derivative c and d) of evaluate()?

2. Why aren’t forces applied/calculated before or at the top of the evaluate()s?

3. In the networked physics example you have some sort of “simple” (but I don’t understand it obviously collision detection and response which works just by applying certain forces.

But how can you calculate the right force to make the cube jump back form the wall. (f = m * a) and the next evaluate call does (momentum += derivative.force * dt), but you don’t know how big that dt will be, do you?

From what I think, wouldn’t it be right to detect a collision after completely integrating with RK4, calculate the fraction of dt at which the collision happened (with the data from the previous state, but I’m not sure how this calculation would look like exactly), integrate to t+collDTime, inverse the velocity/momentum of the body and integrate from the time of collision + (dt-collDTime)?

jezz, I really seem to have no idea *hehe*

I’m no mathematician or physicist but I am a programmer so I’ll answer the (1.) question if you didn’t know already:

(const State& initial) is a constant reference to the parameter passed in, all well and good we know that won’t alter the original value.

(State state) is a pass-by-value new instantiation of the structure – it won’t effect the actual values of the object passed in, merely duplicate them for the scope of the function, one side effect of pass-by-value is that every object passed in is re-constructed, hence being slow for large blocks of data.

Saying that I can’t get his code adapted to my project, the intial “push” of the torque exists but the object begins to stop spinning!

I was puzzled by this too – it is really easy to overlook the difference between “State& state” and “State state”.

Hey glenn i would REALLY appreciate it if you could help me: Im going to try and make it as easy as possible for me to understand: Say you have a quaternion q = 1 + 2i +3j + 4k. How would you integrate this with RK4. How do you calculate the a, b, c, and d? Please help me. I need to do a INS system with different quaternions arriving every like 1ms that need to be integrated with RK4. please try explain not with C++,just with maths, cause i need to use matlab. (dont really know C++ that well). Thanx

the technique in this article for integration only works for *unit quaternions*, eg. quaternions that represent a rotation

i dont know how to integrate the quaternion example you gave, my maths is not that good

google for it! check your math book! – its what i would have to do

Nice articles! I’d like to note one thing I don’t understand though:

When updating the spin quaternion, you multiply two quaternions and a scalar. So, how do you multiply a Quaternion and a scalar? Maybe you only multiply the W component?

You are my god!

First of all great article.

Second of all is your matrix implementation compatible with D3DXMATRIX?

I have converted your stuff for use with directX but can’t for the life of me figure out why I can’t implement the calculate fustrum planes function. I have the correct code (I think).

Could this be due to an incompatibility?

nope its based on an OpenGL projection matrix

When applying a force at a point other than the center of mass, a force and a torque are generated per the equations from your “forces and torques” section:

Flinear = F

Ftorque = F x (p – x)

What about when applying a torque at a point other than the center of mass? I must imagine this also generates a force and a torque. What would the equations be for that?

first you’d have to define how that is expected to behave, for me personally i don’t exactly see what applying a torque at a point outside the center of mass would do? I’m not sure its even defined… !

So, say a hammer is floating in zero gravity. The center of mass is probably somewhere near the head. What if I want to apply a rotational force (I probably can’t actually call it torque) to the base of the handle? It seems like it would generate some torque about the center of mass as well as a linear force.

If one pushes down on a bowling ball, on a greased surface, at a point to, say, the exact north of the center of mass, the linear component of the force would entirely cancel with the normal force from the surface, and the ball would rotate without moving. New example needed

___

/ \ ____ push down at this point

\___/

no it wouldnt, because there is no rotational component to this force, by definition

“vpoint = vlinear + vangular cross (p – x)”

What exactly does your “cross” stand for?

I’m working on predictive collision detection, and knowing where points will be according to the addition of Linear Velocity and Angular Velocity is something I need a little help in understanding.

Linear is retardedly easy, but I’m not making too much sense of Angular Velocity (in the physics sense)- I understand what it represents, but how do you apply that exactly to any given point on an object?

You need to find the distance between the point and axis of rotation (or a perpendicular line), and then move it along the arc according to the velocity and the timestep (as would with 2d radial arc movements)- but what’s the simple easy to do formula that gives me what I want? What exactly is going on in this “cross(p-x)”?

If you could explain what’s going on with velocities at any given point on a rigid body that would great!

the cross product between two vectors

http://en.wikipedia.org/wiki/Cross_product

btw. in essence the cross product in the formula above is converting the angular velocity to the tangential velocity for a point, according to the distance from the axis of rotation – so it is the solution you seek

cheers

Maybe it’s just the notation that’s getting me…

Is the Point and Center of Mass being subtracted, and then cross multiplied with the AngularVelocity?

Is that the same as moving the point around the axis of rotation?

In this sense, you’re basically using the angular velocity, the position of the point, and the position of the center of mass, to rotate the point onto a new position? And then you can add Linear Velocity as normal?

yes exactly, so what the cross product gives you is the instantaneous linear velocity for the point that occurs from the rotation – you add the actual linear velocity to this, and you have the velocity of the point at the instant

study the cross product, its pretty essential – dot and cross, once you have both nailed most 3d math is pretty easy

cheers

another way of looking at this, given two unit vectors that are linearly independent (eg. they arent both along the same line…), if you cross these then you get another vector that is perpendicular to both

this is how the linear vector and the offset from the origin are transformed into the tangential velocity of the point

cheers

great.

CurrentPointPosition + LinearVelocity*(TimeStep*t) + AngularVelocity*(TimeStep*t) X (PointPosition – CenterOfMass)

t is an alpha that represents all values between 0 and 1. I need to scale AngularVelocity before it becomes linear, otherwise it will move along the wrong angle. So this expression should give me all values of a point on a rigid body as it travels from one discrete time step to another (of course, the linear velocity and angular velocity are both also dependent upon t and it’s relationship with the other forces that produce the velocity values).

That makes sense. Thanks!

Oh wait! I think I have a problem.

I’m using this alpha t variable in a distance equation between two points. By setting the distance to zero and solving for t, I know exactly when two objects would collide (if it’s between 1 and 0, then they collide in the given timestep). The equation is already pretty complicated, but there’s just one thing I need to know-

Do I NEED to apply the Scalar timestep BEFORE I cross with the difference of Point/Center?

My head tells me I do, but I’m not really sure. Since I need to isolate t in my expression, isolating t out of a dot product gets a little more complicated than I would like (I can figure that out though).

yes i think what you are doing is much more complex than this, remember of course that the linear motion of the point is just an approximation – it will be quite incorrect for large timesteps, or fast rotating objects, as rotation is non-linear

also, it does not take into account changes in the object’s linear or angular velocity over the timestep, so it will be inaccurate for bodies under acceleration or torque

in short this is a complicated problem i dont know how to solve try searching for “continuous collision detection” and you should find some resources to help you

cheers

Are the functions like integrate and evaluate for any 3d object or are they specific to the cube object? And what about the forces function?

in the source code for this article? well the collision forces are specific for a cube, everything else is general — although you’d want to add an inertia tensor specific to your object shape, if you wish to be accurate

cheers

Ok thanks i will look into inertia tensors.

“Unlike linear velocity, there is no guarantee that angular velocity will remain constant over time in the absence of forces.”

I must not fully understand, because I don’t see any problems with the following:

torque = 0

angular_accel = torque/moment_of_inertia = 0/moment_of_inertia = 0

angular_velocity = integrate(angular_accel by time) = integrate(0 by time) = 0 + constant_of_integration = initial_angular_velocity

inertia is not a single value like mass

instead it is a matrix, so when it is non-uniform (eg. the object is not a sphere…), the angular velocity will vary as the object rotates

this occurs *even* if the angular momentum is remains constant

what this means is that angular momentum is conserved (remains constant in the absence of external forces/torques…) while angular velocity is not

google for “conservation of angular momentum” for more info

cheers

I see now. Thanks for explaining that.

I did actually search for conservation of angular momentum (you mentioned it in one of the other articles), but I didn’t see how it applied to the situation.

I should have bothered to read the other posts, many people have stated things on analytical approaches, my post would have been a bit different.

About the springs, in the article i stated that you have to keep your framerate high enough to get a decent environment interaction.

I see springs as “envirenment”. In the simulation i made, i did use springs, and yes, when the framerate went down, they did get unstable.

The problem is, that the springs themselves are also objects integrated by the Runge Kutta 4 integrator. And that their forces depend on their location and their velocity.

For a correct integration, you need midframe derivative evaluation. Which is dependent on the location and the velocity of the springs.

The problem is, you don’t know these, because there is no mid frame state. You need something like a mid-frame state. When your are using functions to model the derivatives, you do have this information.

Springs are dependent on the integration function itself, you therefore do not know the mid frame position and velocity of a spring.

I don’t think you can get the full RK4 accuracy with springs.

You could enhance it by calculating a midframe state with for example a Maclaurin polynomial. I don’t know if this is a good idea.

Is there a way to correct this?

The program i wrote simulated a cartoonish character (using lots of springs to tie it together). It used RK4 integration. It ran stable at 60fps.

Is this framerate too high, did i do something wrong?

The spring forces were calculated once per frame.

That should be the cause why a simple polynomial gave the same results (visually, and by doing some testing).

I have learned a lot from this site. Springs, introduction to moment of inertia, and not to forget Runge Kutta 4. When i used it, a year ago, i was very much under the impression that i used it correctly. (The results were pleasing.)

I’m certainly planning on having another look at it.

I do however have my doubts on it’s applicability in non functionally modelled simulations.

The problem that i am having, is that there is a moment in time, a moment in which you know position, velocity, and accelleration. The logical conclusion of such an event is Newtonian physics.

In this case a Maclaurin polynomial is a very good Newton approved method of predicting the future.

My view is, is that most simulations, only have a now state, and now, needs to advance to the next state.

You don’t know anything about the future, how forces will be, just the present is known.

Is this a totally wrong way of looking at simulations?

Greetings

it can’t be fixed, you just need to choose reasonable timesteps given the tightness of the springs

for what its worth, RK4 is basically attempting via multiple sampling (4X) per-frame to discover the derivatives of the function, then reconstruct a taylor’s series approximation with 5th order error

if you know all the derivatives, of course you can just do a completely analytical solution, but here assuming the simulation is controlled by the user your best bet (for general numerical integration) is the RK4 method

of course, many other methods exist, and i’m not even going to suggest that RK4 is the best for springs, just that its a lot better than explicit euler

keep googling for more info!

cheers

A way to solve the spring problem could be to look at the last two or three frames that the step based simulation generated. By constructing a polynomial that fit’s the pasts derivatives, you can predict those of the near future.

For three states, this would lead to a quadratic polynomial (of accelleration).

This approach will smoothen the eventual path of the integrated object. It will cause some delay, if the accelleration changes a lot in an instant, the effects will be smoothened. It remains to be seen if this is desirable.

Perhaps this is a good way to fit springs, because springs behave like sine’s, a quadratic accelleration polynomial gives a better fit than a linear, or a constant one.

This is only an idea, i’d have to implement it to see if it really works.

Greetings

i’m not 100% sure, but what you are proposing sounds awfully a lot like reinventing RK4

RK4 is built on top of a taylor’s series expansion, the trick is that since you don’t know the derivatives of the function, it samples the function 4 times per-timestep to approximate them

these approximations are good enough to get the order of the error down to order 5

awesome! simply awesome! Thank you very much for such an awesome set of articles!

I think an RK4 integrator (like it is implementented in for example your “PhysicsIn3D” demo ) is a good way to approximate the integral of any function, no matter the complexity.

My understanding of RK4 is based on this picture:

http://rainman.astro.uiuc.edu/ddr/ddr-galaxy/pics/rk4.JPG

The picture clearly suggests that you need derivatives of the function. These derivatives are used in your “PhysicsIn3D demo”.

A also looked at your “Networked Physics” demo. This is a more real life example of a computer simulation (or game). Forces that are applied during the frame are unknown. The function “forces”, takes no time argument. This is because these forces are applied once per frame. Because the acceleration is constant, you cannot do better than a 2nd degree Maclaurin polynomial. I actually think your implementation in Network Physics is equivalent to a 2nd degree Maclaurin polynomial, although i haven’t proved it. If you would add (subframe) time dependent functions to the forces function, then an RK4 implementation is better.

My program looks a bit like Networked Physics. I used a Maclaurin polynomial because it is much easier to implement.

Today i created a solution to the inaccuracy of the acceleration. Normally, the acceleration is constant during a frame. There is no information about the curvature of the acceleration.

I used a simple trick, to get a better estimation of the acceleration function. In my previous post, i suggested a quadratic extrapolation. Today, i used a linear one, which works fine.

the acceleration is defined as follows:

(da/dt) = ((acceleration current frame) – (acceleration previous frame))/timestep.

a(t) = a(0) + t * (da/dt)

{Let’s assume the frame ranges from 0 to 1, timestep = 1}

An example on how this will cause the simulation to behave:

t=0:acceleration is zero

t=1:{here we begin, because the algorighm needs a prev state}

set acceleration 10 (a correct value)

End of frame: acceleration = 20 (10 + (10 – 0))

t = 2 {the value of acceleration is discarded because the correct value of acceleration is known)

acceleration = correct value (accumulated by algorithms before the integration step)

This linear prediction is better than assuming that the acceleration does not change.

The reason i think this, is because i tested it.

Test:

A sphere lies on a plane. Gravity is pulling the sphere down. The ground uses “penalty” forces to keep the sphere on the surface. In a perfect simulation, the sphere would lay still on the plane. It might sink a little bit, but it should be resting.

In an imperfect numeric simulation, the sphere would be to far inside the ground. Therefor the penalty forces would become too large (depth dependent) and the sphere would bounce.

I tested the algorithm with a constant force and with the linear estimated force.

The test statistic was the maximum height of bouncing achieved after 20 seconds.

The constant acceleration scored a value of 6.71736.

The linear estimation scored a value of 0,0162495.

This is just one way of testing the integrator. I implemented some springs, and they look fine too.

Both of the integrator’s used Maclaurin polynomials. The constant a 2nd degree, and the linear a third degree.

You can do the same trick with Runge Kutta 4 too.

Reviewing the first part of this post; when you only have constant acceleration, your function is shaped like a 2nd degree polynomial. The solution is a Maclaurin polynomial. No matter how you average subderivatives, you cannot do better than the analytical solution.

Greetings

actually, the forces are not constant over the frame because they depend on the position and velocity of the body – this is why i use RK4…

Today, i finished a 3D demo which uses a third order polynomial to integrate paths of objects.

You can find it at:

http://www.wouterbuddingh.nl/index.php?page=projecten/physics/physics.php

The simulation has an integration step of 1/100th of a second. There is no friction between objects, and no torque is applied to the spheres.

I have a question about applying forces. Do you apply a different force for the torque and the linear force?

In this equation:

Flinear = F

Ftorque = F x (p – x)

If for example you apply a force along the z-axis, the box would rotate around the z-axis and move along the z-axis. But we would only want to rotate around the z-axis. So do you have to dampen the Flinear?

It depends entirely on what you want to achieve. But yes, you generally apply damping to both linear and angular motion

ps. If you only want to rotate around the z-axis you would just apply a linear force to counteract the linear component of the force you apply, or to apply the torque directly (it has no linear component, by definition).

Thank you for the reply,

now I understand why the object moves along the axis while rotating. I also read the first ‘x’ in Ftorque = F x (p – x) as the x-position not the cross product symbol x.

I noticed that you defined torque with switched point offset and force. It should be T = p x F, not T = F x p, which leads to a torque in the reverse and wrong direction.

Looks like you are right, sorry about that

Very neatly explained!

Great tutorial & easy to understand !

Thanks

Awesome Tutorial! Thank you very much

Hi Glenn Fiedler,

I am really inspired with your work. I am an undergrad student and i am doing a project in 3d, i want to make a simulator which shows random bodies movement in 3D. I want to know which base model should i select to go further.

Thanks

Hi and thank you for an awesome set of tutorials! I have been trying to get my head around the RK4 integrator and was wondering how you would modify the example if you had some force that is dependent on velocity? Should each evaluate() call update the velocity with Euler integration so to get the curvature for this force too?

ideally, yes

I found this when googling for additional information about simulating force and momentum as well as torque and angular momentum. Most of it seems to be nearly a carbon copy of this article, and your name isn’t on it anywhere…

http://challenge.nm.org/archive/05-06/finalreports/52.pdf

No big deal

This is such a great program. I was able to the program perfectly yesterday but now I am receiving linking errors such as:

[Linker error] undefined reference to `SetPixelFormat@12′

[Linker error] undefined reference to `wglCreateContext@4′

[Linker error] undefined reference to `wglMakeCurrent@8′

[Linker error] undefined reference to `SwapBuffers@4′

which all seem to come from the “windows.h” file

I am using bloodshed dev-c++

Looks like you are not linking to OpenGL!

Glenn, thanks for a tutorial written much more clearly than just about any other.

Are the vectors such as angularVelocity expressed in the coordinate space of the spinning object, or in that of its parent space? i.e. if the spinning object is a toy airplane and the toy airplane is rolling, is its angularVelocity = (0, 0, rollRate), assuming the toy’s forward is defined as the Z axis, regardless of its orientation in world space?

Angular and linear velocity are in world space

I’m making headway here, but I’m not out of the woods.

Each tick of my simulation, I have been computing angular acceleration on my body independently for each of its own local X, Y and Z axes. That is, my inter-tick rotation vector’s three dimensions express the speed of rotation around each of three axes.

I want to integrate the small incremental rotations on the body, but I see that your use of a vector to record angular velocity is different: yours has a single axis of rotation and a magnitude for rotational rate and direction.

Is there any way I can convert my vector type to yours or is my very approach infeasible due to the non-commutative nature of rotation? If the latter, I’m sort of in a corner, and I would then not know how to compute and sum up the combined effect of the many 3D torques I apply in my tick loop (I am modeling buoyancy of a complex shape).

Tracking separate rotation amounts around x,y,z doesn’t make any sense.

You can only be rotating about a single axis at any one point in time.

That’s why the standard representation of angular velocity is a vector that has direction of the axis of rotation, and length the rate of rotation in radians per-second.

You can still apply rotational acceleration about x,y,z axis btw… just keep the representation of angular velocity standard and you’ll be fine.

Thanks, Glenn. This shoots my idea out of the water, sadly. I’d have hundreds of forces being applied at different points and in different directions each update cycle. If the forces could have been decomposed into local space, the would be only three moments of inertia to calculate (in advance), and torques could be added into three separate ones.

Lacking this means, it seems like the computational complexity would be much larger, as I’d have to calculate a separate moment of inertia each force would act against, and even then I’m not sure I’d understand how to combine these separate resultant effects into a composite axis and acceleration.

I’m not sure how the physics libraries solve this, but then again I’ll notice that they do not appear to calculate a true moment of inertia.

Just had to say… Wow these articles are very well written. I’ve been reading about this stuff for several months and these are the best articles I’ve seen on the subject. Also reading the code doesn’t make me want to shoot myself in the face which is rare.

Thanks!

Where is the normalize method for the quaternion?

const float length = sqrt( x*x + y*y + z*z + w*w );

x /= length;

y /= length;

z /= length;

w /= length;

Cheers.

I found it in the source code before you replied though:)

I’m kind of confused about using forces. I have used these articles to add physics movement to my game objects, but I don’t understand how to actually apply forces to an object.

Forces are vectors. They change the motion by applying an acceleration proportional to the mass of the object. If an object A has mass 1, and object B has mass 2, applying the same force applies 1/2 the acceleration to B. Hence, F = ma, or a = F/m. cheers

Say i have multiple force acting on an object. lets say the object needs to move towards the player and there are also explosions triggering around that apply impulses to the object, would doing something like this in the forces method be correct….

static void Forces(const MotionState &state, float t, AGKVectorExt &force, AGKVectorExt &torque) {

force = AGKVectorExt(0, 0, 0);

itterate a vector/map of forces {

force += anotherForce;

if life < 0 delete the force;

}

}

I was thinking of having a map or vector of forces for each object and maybe a seperate map or vector for global forces that act on every object. A force being its own type with a lifetime, so in your object you can do something like AddForce(vector force, float lifetime)

For moving the object to the player i could set the forces lifetime to the deltatime so it adds a new force every update to push the object towards the player, or just set the existing force to a new value.

Do you think this would work ok?

Sure. Typically you would want to accumulate forces per-frame, calculate their magnitude and apply them to the object, eg. forces due to gravity, air resistance etc.

Alternatively, you can do everything with impulses and just apply instantaneously throughout the frame.

Really it just depends on your integrator. Using RK4 and handling forces you generally want a function that calculates all forces per-timestep.

This can be annoying and is often overkill. Maybe use implicit euler or verlet integration instead?

cheers

“Alternatively, you can do everything with impulses and just apply instantaneously throughout the frame.”

Isn’t that what the forces method in your code is is doing already? Just with a single force?

Kind of, but in the case of RK4 the forces method is called multiple times per frame, not once

Yeah. By throughout the frame I thought you meant multiple times throughout the frame.

What i have is an ApplyForce method that doesn’t actually apply a force, but adds a force to a list which is applied to the object during each physics update (multiple times per frame).

I was wondering, how would I go about capping the velocity of the object to a minimum and maximum value and maybe adding friction/air friction? The enemies chasing my player need to take sharper turns (almost turn instantaneous). At the moment they are building up their velocity, whizzing past the player and then taking ages for the changed force to counter the velocity buildup.

I guess I should rename that method to AddForce instead:)

I gather, by using angular momentum, that this implementation includes gyroscopic effects naturally? I’ve been testing applying forces at an offset on a rapidly spinning object and precession seems to work inherently.

I’ve noticed similar behavior with my go stone, but I don’t know for sure.

I’m writing a ballistics simulator and while physics calculations are simple enough this is an important part. Any ideas how do confirm? I was under the impression that i would be calculating gyroscopic forces myself, based on input forces. But if this is inherent that is cool. Though it complicated somethings that i wanted to do.

I actually don’t know, ask Erin Catto or Erwin Coumans. They have the math and physics knowledge to answer such a question.

Will do.

you are the best teacher in the world…………….

When you are using inversemass isnt that supposed to be the reciprical of mass or no?

v = p/m (velocity equals momentum divided by mass)

void recalculate()

{

velocity = momentum * inverseMass;

}

I am getting very confused as I am not sure if you mean inverse when you say inverse or reciprical.

reciprocal is 1/x where as inverse is turning y = f(x) into x = f(y).

So something like log or square root is inverse.

1/mass.