Introduction to Physical Simulations

Introduction To Physical Simulations

If you are familiar to physics and want to start implementing code of physical simulations, this tutorial could help you. In order to benefit from this tutorial you should be familiar to vectoral operations in 3D as well as physical concepts such as force and velocity.

In this tutorial you will find a very simple physical simulation engine. Contents of the tutorial is as follows:

Contents:

The Design:
* class Vector3D	--->	An Object To Represent A 3D Vector Or A 3D Point In Space.
Force and Motion:
* class Mass	--->	An Object To Represent A Mass.
How A Simulation Should Operate:
* class Simulation	--->	A Container Object For Simulating Masses.
Operating A Simulation By An Application:
* class ConstantVelocity : public Simulation	--->	A Simulation Object That Creates A Mass With A Constant Velocity.
Applying Force:
* class MotionUnderGravitation : public Simulation	--->	A Simulation Object That Creates A Mass Which Moves Under Gravitation.
* class MassConnectedWithSpring : public Simulation	--->	A Simulation Object That Creates A Mass Connected To A Point By A Spring.

The Design:

Design of physical simulation engines is not always simple. But there is a simple order of dependency; application depends on simulation toolkit and the simulation toolkit depends on math libraries. Here, we will make use of this simple order. Our purpose is to obtain a container to simulate motion of masses. The simulation toolkit will include an object "class Mass" and an object "class Simulation". "class Simulation" will be our container. When we obtain the Simulation class, we will be able to develop applications. But before that, we need a math library. The library includes only one class "class Vector3D". Vector3D class will be used to represent points, vectors, position, velocity, and force in 3 dimensional space.

* class Vector3D ---> An Object To Represent A 3D Vector Or A 3D Point In Space.

Vector3D class is the only member of our modest math library. Vector3D holds x, y, and z values and it implements operators for vector arithmetics in 3D. Addition, subtraction, multiplication, and division operators are coded in Vector3D. Since this tutorial focuses on physics, I will not go into details of Vector3D. If you have a look at Physics1.h you will see how simple Vector3D is.

Force And Motion:

For implementing physical simulations, we should know what a mass is. A mass has a position and a velocity. A mass has weight on Earth, Moon, Mars, and at any place where gravitation exists. Weight is different on different gravitations of different places. But there is one common value for a mass, which is the same in all conditions. This value is also called mass. Mass value of a mass! Mass value represents "how much a mass exists in space"! For example a book is a mass with a weight of say 1 kg on The Earth and with a weight of 0.17 kg on The Moon and has a mass value of 1 kg everywhere. The mass value is designated to be equal to its mass on The Earth.

After having understood the mass of a mass, we should go on with force and motion. A mass, with a non-zero velocity in space, moves in the direction of the velocity. Therefore, one reason of the change in position is the velocity. Passing of time is another reason. Change in position depends on how fast a mass moves and how much time has passed. You should have understood until here to go on to the next paragraph. If not, spend some time thinking on the relation between position, velocity and time.

Velocity of a mass changes if there is force acting on it. Velocity tends to the direction of the force. This tending is proportional to the force and inversly proportional to mass. The change in velocity per unit time is called acceleration. More the force on a mass, more its acceleration. More the mass value of a mass, less its acceleration. When acceleration is formulated it is:

acceleration = force / mass

From here we obtain the famous equation:

force = mass * acceleration

(We will mostly use the acceleration formula)

For preparing a physical medium to simulate, you should be aware of the environment that the simulation takes place. The environment in this tutorial is simply empty space waiting to be filled by masses we create. The units of the values to represent masses and time shall be decided firstly. I have decided to use the time unit as seconds and units of position values as meters. Accordingly, unit of velocity becomes meters per second (m/s). And the unit of acceleration becomes meters per second per second (m/s/s) = ((m/s)/s) = (m / (s*s))! (this means velocity per second since velocity is meters per second) I have decided the unit of mass values as kilograms (kg).

* class Mass ---> An Object To Represent A Mass.

Now we are starting to use the theory! We have to write a class to represent a mass. It should hold the mass value, the position, the velocity, and the force applied at an instance.

class Mass
{
public:
    float m;                                // The Mass Value.
    Vector3D pos;                               // Position In Space.
    Vector3D vel;                               // Velocity.
    Vector3D force;                             // Force Applied On This Mass At An Instance.
 
    Mass(float m)                               // Constructor.
    {
        this->m = m;
    }
 
    ...

We want to apply force to this mass. At an instance in time, there might be several sources of external forces acting on the mass. The vector sum of these forces gives the net force on the mass at that instance. Before starting to apply forces, we should reset the force on the mass. Then we can add external forces on the mass.

(class Mass continued)
 
void applyForce(Vector3D force)
{
    this->force += force;                        // The External Force Is Added To The Force On The Mass.
}
 
void init()                             // This Method Sets The Force Values To Zero.
{
    force.x = 0;
    force.y = 0;
    force.z = 0;
}
 
...

There is a simple order of things to do in a simulation:

Reset the force (see the init() method())
Apply external forces
Iterate time by "the change in time"

Here, iterating the time is implemented with "The Euler Method". The Euler Method is a simple simulation method. There are more sophisticated methods for simulations. But Euler is good enough for lots of applications. Most of computer and video games use The Euler Method. What this method does is that it calculates the next velocity and next position of a mass according to the force applied and time passed. The iteration is done in void simulate(float dt):

(class Mass continued)
 
void simulate(float dt)
{
    vel += (force / m) * dt;                    // Change In Velocity Is Added To The Velocity.
                                    // The Change Is Proportinal With The Acceleration (force / m) And Change In Time.
 
    pos += vel * dt;                        // Change In Position Is Added To The Position.
                                    // Change In Position Is Velocity Times The Change In Time.
}

How A Simulation Should Operate:

In a physical simulation, at every iteration, the same process takes place. Forces are set to zero, forces are applied, new positions and new velocities are calculated. This process cycles as long as we want the time to pass. This process is implemented in "class Simulation".

* class Simulation ---> A Container Object For Simulating Masses.

Simulation class holds masses as its members. The role of the class is to create and delete masses, and maintain the simulation procedure.

class Simulation
{
public:
    int numOfMasses;                            // Number Of Masses In This Container.
    Mass** masses;                              // Masses Are Held By Array Of Pointers. (Here Mass** Represents A 1 Dimensional Array).
     
    Simulation(int numOfMasses, float m)                    // Constructor Creates Some Masses With Mass Values m.
    {
        this->numOfMasses = numOfMasses;
         
        masses = new Mass*[numOfMasses];                // Create An Array Of Pointers.
 
        for (int a = 0; a < numOfMasses; ++a)                // We Will Step To Every Pointer In The Array.
            masses[a] = new Mass(m);                // Create A Mass As A Pointer And Put It In The Array.
    }
 
    virtual void release()                          // Delete The Masses Created.
    {
        for (int a = 0; a < numOfMasses; ++a)                // We Will Delete All Of Them.
        {
            delete(masses[a]);
            masses[a] = NULL;
        }
             
        delete(masses);
        masses = NULL;
    }
 
    Mass* getMass(int index)
    {
        if (index < 0 || index >= numOfMasses)                // If The index Is Not In The Array.
            return NULL;                        // Then Return NULL.
 
        return masses[index];                       // Get The Mass At The index.
    }
 
...

The simulation procedure has three steps:

init() to set forces to zero
solve() to apply forces
simulate(float dt) to iterate masses by the change in time

(class Simulation continued)
 
virtual void init()                         // This Method Will Call The init() Method Of Every Mass.
{
    for (int a = 0; a < numOfMasses; ++a)                // We Will init() Every Mass.
        masses[a]->init();                   // Call init() Method Of The Mass.
}
 
virtual void solve()                            // No Implementation Because No Forces Are Wanted In This Basic Container.
{
                                    // In Advanced Containers, This Method Will Be Overridden And Some Forces Will Act On Masses.
}
 
virtual void simulate(float dt)                     // Iterate The Masses By The Change In Time.
{
    for (int a = 0; a < numOfMasses; ++a)                // We Will Iterate Every Mass.
        masses[a]->simulate(dt);             // Iterate The Mass And Obtain New Position And New Velocity.
}
 
...

The simulation procedure is packed into one method:

    (class Simulation continued)
 
    virtual void operate(float dt)                      // The Complete Procedure Of Simulation.
    {
        init();                             // Step 1: Reset Forces To Zero.
        solve();                            // Step 2: Apply Forces.
        simulate(dt);                           // Step 3: Iterate The Masses By The Change In Time.
    }
};

By now, we have a simple physical simulation engine. It is based on a math library. It contains Mass and Simulation classes. It uses a very common pattern of simulation procedure and it uses The Euler Method. Now we are ready to develop applications. The applications that we will develop are:

Mass with constant velocity
Mass under gravitational force
Mass connected to a still point by a spring

Operating A Simulation By An Application:

Before we write a specific simulation, we should know how to operate simulations by applications. In this tutorial, the simulation engine and the application to operate the simulations are seperated in two files. In the application file there is a function as:

1	`void` `Update (DWORD` `milliseconds)` `// Perform Motion Updates Here.`

This function is called repeatedly at every frame update. The "DWORD milliseconds" is the time period from the previous frame to the current frame. From this, we can say that we should iterate simulations according to the "milliseconds". If the simulations follow this time period, they should go parallel with the real world's time. To iterate a simulation, we simply call the "void operate(float dt)" method. To call this method, we should know "dt". Since we take the time unit as seconds we firstly convert milliseconds to seconds (see below in the code). Then we use a value "slowMotionRatio" which means, how slow we want to run the simulation relative to the real world time. We divide dt by this value and we obtain a new dt. Now we can add dt to "timeElapsed". "timeElapsed" is the time of the simulation, not the time of the real world.

void Update (DWORD milliseconds)
{
    ...
    ...
    ...
 
    float dt = milliseconds / 1000.0f;                  // Let's Convert milliseconds To Seconds.
 
    dt /= slowMotionRatio;                          // Divide dt By slowMotionRatio And Obtain The New dt.
 
    timeElapsed += dt;                          // Iterate Elapsed Time.
 
    ...

Now dt is almost ready for operating the simulation. But... There is one very important thing that we should know: dt is the detail of precision. If dt is not small enough, your simulation would show instability and the motion would not be calculated precisely. Stability analysis is used for physical simulations to find the maximum dt value that a simulation can handle. In this tutorial we will not go into the details and if you are just developing a game but not a scientific application, it is always valid to find the value of maximum dt by trial and error.

As an example, in a car racing game, it is convenient to use dt as about 2 to 5 milliseconds for regular car, and 1 to 3 milliseconds for a formula car. In an arcade car simulation it is possible to use dt as about 10 to 200 milliseconds. Less the value of dt, more the CPU ticks we need to catch up the real world time. That is why physical simulations are rarely used in older games.

In the code below we define the maximum possible dt as 0.1 seconds (100 milliseconds). With this value we will calculate the number iterations to be made at the current update. We write a formula:

int numOfIterations = (int)(dt / maxPossible_dt) + 1;

numOfIterations is the number of iterations to be made for a simulation. Say that the application is running with 20 frames per second, which gives dt = 0.05 seconds. Then numOfIterations becomes 1. The simulation will be iterated once by 0.05 seconds. Say dt was 0.12 seconds. Then numOfIterations is 2. Below, just after "int numOfIterations = (int)(dt / maxPossible_dt) + 1;", you should see that dt is calculated once again. There, dt is divided by numOfIterations and it becomes dt = 0.12 / 2 = 0.06. dt was originally more than the maximum possible value 0.1. Now we have dt as 0.06 but we will operate the simulation twice so that we catch up 0.12 seconds as a result. Examine the code below and be sure that you understand the paragraph above.

    ...
 
    float maxPossible_dt = 0.1f;                        // Say That The Maximum Possible dt Is 0.1 Seconds.
                                        // This Is Needed So We Do Not Pass Over A Non Precise dt Value.
 
    int numOfIterations = (int)(dt / maxPossible_dt) + 1;           // Calculate Number Of Iterations To Be Made At This Update Depending On maxPossible_dt And dt.
    if (numOfIterations != 0)                       // Avoid Division By Zero.
        dt = dt / numOfIterations;                  // dt Should Be Updated According To numOfIterations.
 
    for (int a = 0; a < numOfIterations; ++a)                // We Need To Iterate Simulations "numOfIterations" Times.
    {
        constantVelocity->operate(dt);                   // Iterate constantVelocity Simulation By dt Seconds.
        motionUnderGravitation->operate(dt);             // Iterate motionUnderGravitation Simulation By dt Seconds.
        massConnectedWithSpring->operate(dt);                // Iterate massConnectedWithSpring Simulation By dt Seconds.
    }
}

Let's begin to write the applications:

1. Mass with constant velocity
* class ConstantVelocity : public Simulation ---> A Simulation Object That Creates A Mass With A Constant Velocity.

Mass with constant velocity does not need any external force. We just have to create 1 mass and set its velocity to (1.0f, 0.0f, 0.0f) so that it moves in the x direction with a velocity of 1 m/s. We will derive a class from Simulation. This class is "class ConstantVelocity":

class ConstantVelocity : public Simulation
{
public:
    ConstantVelocity() : Simulation(1, 1.0f)                // Constructor Firstly Constructs Its Super Class With 1 Mass And 1 Kg.
    {
        masses[0]->pos = Vector3D(0.0f, 0.0f, 0.0f);         // A Mass Was Created And We Set Its Position To The Origin.
        masses[0]->vel = Vector3D(1.0f, 0.0f, 0.0f);         // We Set The Mass's Velocity To (1.0f, 0.0f, 0.0f) m/s.
    }
};

When the operate(float dt) method of ConstantVelocity class is called, it calculates the next state of the mass. This method is called by the main application, before every re-draw of the window. Say that your application is running with 10 frames per second. If the exact time change was passed to operate(float dt) at every frame, then dt would be 0.1 seconds. When the simulation calls simulate(float dt) of the mass, its new position will be incremented by velocity * dt which is

Vector3D(1.0f, 0.0f, 0.0f) * 0.1 = Vector3D(0.1f, 0.0f, 0.0f)

At every iteration, the mass moves 0.1 meters to the right. After 10 frames, it will have moved 1.0 meter to the right. The velocity was 1.0 m/s and it moves 1.0 meter in one second. Was that a coincidental or a logical result? If you can't answer this question spend some time thinking about the relation above.

When you run the application, you see the mass with constant velocity moving in the x direction. The application provides two modes of motion flow. By pressing F2 you get time flow parallel to the real world. And by pressing F3 you get time flow 10 times slower than the real world time. On the screen you will see lines to represent the coordinate plane. The spacing between these lines is 1 meter. By the use of the lines, you can observe that the mass moves 1 meter in a second when it is set to the real world time mode. And in the slow mode it moves 1 meter in ten seconds. The techique described above is a common one to make the real-time simulations flow parallel to the real world time. To use this technique, you must have strictly decided on the units of your simulation.

Applying Force:

In the constant velocity simulation, we did not apply any force to the mass. Because we know that if a force acts on a body, it accelerates. When we want accelerating motion, we apply forces. At one operation of a simulation, we apply forces in the "solve" method. When the operation comes to the "simulate" phase, the net force results as the total of forces. The net force determines the motion.

Say that you want to apply 1 N force on a mass in the x direction. Then you should write:

mass->applyForce(Vector3D(1.0f, 0.0f, 0.0f));

in the "solve" method. If you want to add another force say 2 N in the y direction you should write:

mass->applyForce(Vector3D(1.0f, 0.0f, 0.0f));
mass->applyForce(Vector3D(0.0f, 2.0f, 0.0f));

in the "solve" method. You can add any force in any formulated way and you obtain the motion. In the next application you will see a single force in a formulated way.

2. Mass under gravitational force
* class MotionUnderGravitation : public Simulation ---> A Simulation Object That Creates A Mass Which Moves Under Gravitation.

MotionUnderGravitation class creates a mass and it applies a force to it. This force is the force of gravitation. Force of gravitation is equal to the mass times the gravitational acceleration:

F = m * g

Gravitational acceleration is the acceleration of free body. On the earth, when you drop an object, it gains 9.81 m/s velocity every second, unless it experiences a force other than the gravitational force. Therefore the gravitational acceleration is constant for all masses on the earth and it is 9.81 m/s/s. (This is independent of the mass. All masses fall with the same acceleration.)

MotionUnderGravitation class has such a constructor:

class MotionUnderGravitation : public Simulation
{
    Vector3D gravitation;                           // The Gravitational Acceleration.
 
    MotionUnderGravitation(Vector3D gravitation) : Simulation(1, 1.0f)  // Constructor Firstly Constructs Its Super Class With 1 Mass And 1 Kg.
    {                                   // Vector3D Gravitation, Is The Gravitational Acceleration.
        this->gravitation = gravitation;             // Set This Class's Gravitation.
        masses[0]->pos = Vector3D(-10.0f, 0.0f, 0.0f);           // Set The Position Of The Mass.
        masses[0]->vel = Vector3D(10.0f, 15.0f, 0.0f);           // Set The Velocity Of The Mass.
    }
 
    ...

The constructor gets a Vector3D gravitation, which is the gravitational acceleration and the simulation uses that in the force to be applied.

virtual void solve()                            // Gravitational Force Will Be Applied Therefore We Need A "Solve" Method.
{
    for (int a = 0; a < numOfMasses; ++a)                // We Will Apply Force To All Masses (Actually We Have 1 Mass, But We Can Extend It In The Future).
        masses[a]->applyForce(gravitation * masses[a]->m);    // Gravitational Force Is As F = m * g. (Mass Times The Gravitational Acceleration).
}

Above, in the code, you should see the force formula as F = m * g. The application creates MotionUnderGravitation with a Vector3D value as "Vector3D(0.0f, -9.81f, 0.0f)". -9.81 means an acceleration in the negative y direction so that the mass falls down when displayed on the window. Run the application and observe what's happening. Write 9.81 instead of -9.81 and observe the difference.

3. Mass connected to a still point by a spring * class MassConnectedWithSpring : public Simulation ---> A Simulation Object That Creates A Mass Connected To A Point By A Spring.

In this example, we want to connect a mass to a still point by a spring. The spring should pull the mass towards the connection position, so that it oscillates. In the constructor, MassConnectedWithSpring class sets the connection position and sets the position of the mass.

class MassConnectedWithSpring : public Simulation
{
public:
    float springConstant;                           // The More springConstant, The Stiffer The Spring Force.
    Vector3D connectionPos;                         // The Arbitrary Still Point That The Mass Is Connected.
 
    MassConnectedWithSpring(float springConstant) : Simulation(1, 1.0f) // Constructor Firstly Constructs Its Super Class With 1 Mass And 1 Kg.
    {
        this->springConstant = springConstant;               // Set The springConstant.
 
        connectionPos = Vector3D(0.0f, -5.0f, 0.0f);            // Set The connectionPos.
 
        masses[0]->pos = connectionPos + Vector3D(10.0f, 0.0f, 0.0f);    // Set The Position Of The Mass 10 Meters To The Right Side Of The connectionPos.
        masses[0]->vel = Vector3D(0.0f, 0.0f, 0.0f);         // Set The Velocity Of The Mass To Zero.
    }
 
    ...

Velocity of the mass is set to zero and its position is set to 10 meters to the right side of the connectionPos so that it can be pulled towards left at the beginning. Force of a spring is formulated as

F = -k * x

k is a value to represent how stiff the spring should be. And x is the distance from the mass to the connection position. The negative sign at the formula indicates that this is an attractive force. If it was positive, the spring would push the mass which is not something that we would expect to see.

virtual void solve()                                // The Spring Force Will Be Applied.
{
    for (int a = 0; a < numOfMasses; ++a)                    // We Will Apply Force To All Masses (Actually We Have 1 Mass, But We Can Extend It In The Future).
    {
        Vector3D springVector = masses[a]->pos - connectionPos;      // Find A Vector From The Position Of The Mass To The connectionPos.
        masses[a]->applyForce(-springVector * springConstant);       // Apply The Force According To The Famous Spring Force Formulation.
    }
}

The spring force in the code above is the same as the spring formula (F = -k * x). Here instead of x we use a Vector3D because we want to use 3D space. "springVector" gives the difference of the mass position and the connectionPos and springConstant stands for k. More the value of springConstant, more the force, and faster the mass oscillates.

In this tutorial, I tried to point to the key concepts of physical simulations. If you are interested in physics, you will not have hard time to create new simulations of your own. You can try complicated interactions and obtain very attractive demos and games. The next steps to take will be rigid body simulations, simple mechanisms and advanced simulation methods.

For any comments or questions please contact me:

Erkin Tunca (erkintunca@icqmail.com)

Jeff Molofee (NeHe)

* DOWNLOAD Visual C++ Code For This Lesson.

* DOWNLOAD Borland C++ Builder 6 Code For This Lesson. ( Conversion by Le Thanh Cong )
* DOWNLOAD Code Warrior 5.3 Code For This Lesson. ( Conversion by Scott Lupton )
* DOWNLOAD Delphi Code For This Lesson. ( Conversion by Michal Tucek )
* DOWNLOAD Dev C++ Code For This Lesson. ( Conversion by Warren Moore )
* DOWNLOAD Linux/GLut Code For This Lesson. ( Conversion by Laks Raghupathi )
* DOWNLOAD Visual Studio .NET Code For This Lesson. ( Conversion by Grant James )

< Lesson 38Lesson 40 >