struct Point {
  Vec2 Position;
  Vec2 OldPosition;
  Vec2 Acceleration;
};

class Physics {
  int PointCount;
  Point* Points[ MAX_VERTICES ];

  float Timestep;

public:
  void  UpdateVerlet();

  //Constructors, getters/setters etc. omitted
};

void Physics::UpdateVerlet() {
  for( int I = 0; I < PointCount; I++ ) {
    Point& P = *Points[ I ];

    Vec2 Temp = P.Position;
    P.Position += P.Position - P.OldPosition + P.Acceleration*Timestep*Timestep;
    P.OldPosition = Temp;
  }
}

struct Edge {
  Vertex* V1;
  Vertex* V2;

  float OriginalLength; //The length of the edge when it was created

  //Constructors etc. omitted
};

void Physics::UpdateEdges() {
  for( int I = 0; I < EdgeCount; I++ ) {
    Edge& E = *Edges[ I ];

    //Calculate the vector mentioned above
    Vec2 V1V2 = E.V2->Position - E.V1->Position; 

    //Calculate the current distance
    float V1V2Length = V1V2.Length(); 
    
    //Calculate the difference from the original length
    float Diff       = V1V2Length - E.OriginalLength; 
    
    V1V2.Normalize();

    //Push both vertices apart by half of the difference respectively 
    //so the distance between them equals the original length
    E.V1->Position += V1V2*Diff*0.5f; 
    E.V2->Position -= V1V2*Diff*0.5f;
  }
}

struct PhysicsBody {
  int VertexCount;
  int EdgeCount;

  Vertex* Vertices[ MAX_BODY_VERTICES ];
  Edge*   Edges   [ MAX_BODY_EDGES    ];

  void ProjectToAxis( Vec2& Axis, float& Min, float& Max );

  //Again, constructors etc. omitted
};

void PhysicsBody::ProjectToAxis( Vec2& Axis, float& Min, float& Max ) {
  float DotP = Axis*Vertices[ 0 ]->Position;

  //Set the minimum and maximum values to the projection of the first vertex
  Min = Max = DotP; 

  for( int I = 1; I < VertexCount; I++ ) {
    //Project the rest of the vertices onto the axis and extend 
    //the interval to the left/right if necessary
    DotP = Axis*Vertices[ I ]->Position; 

    Min = MIN( DotP, Min );
    Max = MAX( DotP, Max );
  }
}

bool Physics::DetectCollision( PhysicsBody* B1, PhysicsBody* B2 ) {
  //Just a fancy way of iterating through all of the edges of both bodies at once
  for( int I = 0; I < B1->EdgeCount + B2->EdgeCount; I++ ) { 
    Edge* E;

    if( I < B1->EdgeCount )
      E = B1->Edges[ I ];
    else
      E = B2->Edges[ I - B1->EdgeCount ];
    
    //Calculate the axis perpendicular to this edge and normalize it
    Vec2 Axis( E->V1->Position.Y - E->V2->Position.Y, E->V2->Position.X - E->V1->Position.X ); 
    Axis.Normalize();

    float MinA, MinB, MaxA, MaxB; //Project both bodies onto the perpendicular axis
    B1->ProjectToAxis( Axis, MinA, MaxA );
    B2->ProjectToAxis( Axis, MinB, MaxB );

    //Calculate the distance between the two intervals - see below
    float Distance = IntervalDistance( MinA, MaxA, MinB, MaxB ); 

    if( Distance > 0.0f ) //If the intervals don't overlap, return, since there is no collision
      return false;
  }

  return true; //There is no separating axis. Report a collision!
}

float Physics::IntervalDistance( float MinA, float MaxA, float MinB, float MaxB ) {
  if( MinA < MinB )
    return MinB - MaxA;
  else
    return MinA - MaxB;
}

class Physics {
  struct {
    float Depth;
    Vec2  Normal;
  } CollisionInfo;

  //Everything else omitted
}

bool Physics::DetectCollision( PhysicsBody* B1, PhysicsBody* B2 ) {
  float MinLength = 10000.0f; //Initialize the length of the collision vector to a relatively large value
  for( int I = 0; I < B1->EdgeCount + B2->EdgeCount; I++ ) {
    Edge* E;

    if( I < B1->EdgeCount )
      E = B1->Edges[ I ];
    else
      E = B2->Edges[ I - B1->EdgeCount ];

    Vec2 Axis( E->V1->Position.Y - E->V2->Position.Y, E->V2->Position.X - E->V1->Position.X );
    Axis.Normalize();

    float MinA, MinB, MaxA, MaxB;
    B1->ProjectToAxis( Axis, MinA, MaxA );
    B2->ProjectToAxis( Axis, MinB, MaxB );

    float Distance = IntervalDistance( MinA, MaxA, MinB, MaxB );

    if( Distance > 0.0f )
      return false;
      
    //If the intervals overlap, check, whether the vector length on this 
    //edge is smaller than the smallest length that has been reported so far
    else if( abs( Distance ) < MinDistance ) { 
      MinDistance = abs( Distance );

      CollisionInfo.Normal = Axis; //Save collision information for later
    }
  }

  CollisionInfo.Depth = MinDistance;

  return true; //There is no separating axis. Report a collision!
}

struct {
  float Depth;
  Vec2  Normal;

  Edge*   E;
  Vertex* V;
} CollisionInfo;

bool Physics::DetectCollision( PhysicsBody* B1, PhysicsBody* B2 ) {
  float MinDistance = 10000.0f;
  for( int I = 0; I < B1->EdgeCount + B2->EdgeCount; I++ ) { //Same old
    Edge* E;

    if( I < B1->EdgeCount )
      E = B1->Edges[ I ];
    else
      E = B2->Edges[ I - B1->EdgeCount ];

    Vec2 Axis( E->V1->Position.Y - E->V2->Position.Y, E->V2->Position.X - E->V1->Position.X );
    Axis.Normalize();

    float MinA, MinB, MaxA, MaxB;
    B1->ProjectToAxis( Axis, MinA, MaxA );
    B2->ProjectToAxis( Axis, MinB, MaxB );

    float Distance = IntervalDistance( MinA, MaxA, MinB, MaxB );

    if( Distance > 0.0f )
      return false;
    else if( abs( Distance ) < MinDistance ) {
      MinDistance = abs( Distance );

      CollisionInfo.Normal = Axis;
      CollisionInfo.E      = E; //Store the edge, as it is the collision edge
    }
  }

  CollisionInfo.Depth = MinDistance;

  //Ensure that the body containing the collision edge lies in 
  //B2 and the one containing the collision vertex in B1
  if( CollisionInfo.E->Parent != B2 ) { 
    PhysicsBody* Temp = B2;
    B2 = B1;
    B1 = Temp;
  }

  //This is needed to make sure that the collision normal is pointing at B1
  int Sign = SGN( CollisionInfo.Normal*( B1->Center - B2->Center ) ); 

  //Remember that the line equation is N*( R - R0 ). We choose B2->Center 
  //as R0; the normal N is given by the collision normal

  if( Sign != 1 )
    CollisionInfo.Normal = -CollisionInfo.Normal; //Revert the collision normal if it points away from B1


  float SmallestD = 10000.0f; //Initialize the smallest distance to a high value
  for( int I = 0; I < B1->VertexCount; I++ ) {
    //Measure the distance of the vertex from the line using the line equation
    float Distance = CollisionInfo.Normal*( B1->Vertices[ I ]->Position - B2->Center ); 

    //If the measured distance is smaller than the smallest distance reported 
    //so far, set the smallest distance and the collision vertex
    if( Distance < SmallestD ) { 
      SmallestD = Distance;
      CollisionInfo.V = B1->Vertices[ I ];
    }
  }

  return true;
}

void Physics::ProcessCollision() {
  Vec2 CollisionVector = CollisionInfo.Normal*CollisionInfo.Depth;

  CollisionInfo.V->Position += CollisionVector*0.5f;
}

Vertex* E1 = CollisionInfo.E->V1;
Vertex* E2 = CollisionInfo.E->V2;

float T = ( CollisionInfo.V->Position.X - CollisionVector.X - E1->Position.X )/(  E2->Position.X - E1->Position.X );

float T;
if( abs( E1->Position.X - E2->Position.X ) > abs( E1->Position.Y - E2->Position.Y ) )
  T = ( CollisionInfo.V->Position.X - CollisionVector.X - E1->Position.X )/(  E2->Position.X - E1->Position.X );
else
  T = ( CollisionInfo.V->Position.Y - CollisionVector.Y - E1->Position.Y )/(  E2->Position.Y - E1->Position.Y );

void Physics::ProcessCollision() {
  Vec2 CollisionVector = CollisionInfo.Normal*CollisionInfo.Depth;

  Vertex* E1 = CollisionInfo.E->V1;
  Vertex* E2 = CollisionInfo.E->V2;

  float T;
  if( abs( E1->Position.X - E2->Position.X ) > abs( E1->Position.Y - E2->Position.Y ) )
    T = ( CollisionInfo.V->Position.X - CollisionVector.X - E1->Position.X )/(  E2->Position.X - E1->Position.X);
  else
    T = ( CollisionInfo.V->Position.Y - CollisionVector.Y - E1->Position.Y )/(  E2->Position.Y - E1->Position.Y);

  float Lambda = 1.0f/( T*T + ( 1 - T )*( 1 - T ) );

  E1->Position -= CollisionVector*( 1 - T )*0.5f*Lambda;
  E2->Position -= CollisionVector*      T  *0.5f*Lambda;

  CollisionInfo.V->Position += CollisionVector*0.5f;
}

void Physics::Update() {
  UpdateForces();
  UpdateVerlet();
  IterateCollisions();
}

Someone asked me in the #gamedev IRC channel about how to make a 2d vehicle simulator. Instead of spending all day trying to explain the concepts to him, I decided just to write this tutorial. Please bear with me, this is my first tutorial.

So as I mentioned we’re going to be learning how to make a basic 2d vehicle simulator. We’re going to do it in C# and try do use as few hacks as possible. I’ve broken the process down into three steps. First, we will learn how to setup a basic game application in C#.NET and how to draw some basic graphics (emphasis on basic.) Next, we will learn how to create a rigid body simulator using a simple Euler integrator with a variable time step. And last but not least, we will calculate vehicle forces simulating the tire patch contacting the road. And that’s all there is to it! Let’s get started.

Math Requirements

Phase One

Phase one, as I mentioned, is to create the renderer; something graphical so we can actually see what our simulation is doing. This will make it a lot easier to debug. Create a windows form project in C# and place a picturebox control on it (name it "screen"). This control is where we will display our simulation. We could just start drawing to this screen but we’re going to be using double buffering as well to avoid flicker, so we need to create the back buffer now. That bit of code looks like this.

Graphics graphics; //gdi+
Bitmap backbuffer;
Size buffersize;

//intialize rendering
private void Init(Size size)
{
  //setup rendering device
  buffersize = size;
  backbuffer = new Bitmap(buffersize.Width, buffersize.Height);
  graphics = Graphics.FromImage(backbuffer);
}

//main rendering function
private void Render(Graphics g)
{
  //clear back buffer
  graphics.Clear(Color.Black);

  //draw to back buffer
  graphics.DrawLine(new Pen(Color.Yellow), 1, 0, 1, 5);

  //present back buffer
  g.DrawImage(backbuffer, new Rectangle(0, 0, buffersize.Width, buffersize.Height), 0, 0, buffersize.Width, buffersize.Height, GraphicsUnit.Pixel);
}

Now by default, the graphics of a picturebox control has the origin in the top-left corner, and extends downward for +y and to the right for +x. This is highly unnatural for most cases. In addition to that, it has extremely large units. Since we will be simulating in the metric system, I recommend introducing a scale factor to scale up the simulation and make it much more visible. The transformation looks like this and takes place after Graphics.Clear() is called.

graphics.ResetTransform();
graphics.ScaleTransform(screenScale, -screenScale);
graphics.TranslateTransform(buffersize.Width / 2.0f / screenScale, -buffersize.Height / 2.0f / screenScale);

Now the line should draw starting right at the center of the screen.

Forms Wiring

Up until now, I havn’t explained how to connect all the functions. The first thing you’ll need to do is call the Render function from your on_paint event. Next, you’ll need to create a function that gets called continously to update the simulation. It is preferred to call this function on the Application_Idle event. So create an event handler for Application_Idle and have it call your DoFrame function. Inside this function you’ll need to

1. Process input
2. Update the simulation
3. Invalidate the screen control

The Timer

Since we don’t know how often our DoFrame function will be getting called, we need to code everything to handle a variable time step. To utilize this we must measure the time between DoFrame calls. So I’ll introduce the timer which, very simply, queries the number of milliseconds that have passed since the computer was turned on. So we store this number every frame and on a subsequent frame we compute the difference, which gives us the amount of time that has passed since the last frame. Here is my very simple timer object. Note: you will need to call GetETime in your intialize function in order to clear the timer, otherwise the first call to it will return the amount of time that has passed since the computer was turned on.

class Timer
{
  //store last time sample
  private int lastTime = Environment.TickCount;
  private float etime;

  //calculate and return elapsed time since last call
  public float GetETime()
  {
    etime = (Environment.TickCount - lastTime) / 1000.0f;
    lastTime = Environment.TickCount;

    return etime;
  }
}

Conclusion of Phase One

So up until now we’ve covered: setting up a rendering surface using GDI, wiring a form to process a game loop and draw it to the screen, and computing the time that has passed since the last frame. Our application looks like this:

using System;
using System.Collections.Generic;
using System.ComponentModel;
using System.Data;
using System.Drawing;
using System.Drawing.Drawing2D;
using System.Text;
using System.Windows.Forms;

namespace racing_simulation_2d
{
  //our main application form
  public partial class frmMain : Form
  {
    //graphics
    Graphics graphics; //gdi+
    Bitmap backbuffer;
    Size buffersize;
    const float screenScale = 3.0f;
    Timer timer = new Timer();

    //keyboard controls
    bool leftHeld = false, rightHeld = false;
    bool upHeld = false, downHeld = false;

    //vehicle controls
    float steering = 0; //-1 is left, 0 is center, 1 is right
    float throttle = 0; //0 is coasting, 1 is full throttle
    float brakes = 0; //0 is no brakes, 1 is full brakes

    public frmMain()
    {
      InitializeComponent();
      Application.Idle += new EventHandler(ApplicationIdle);

      screen.Paint += new PaintEventHandler(screen_Paint);
      this.KeyDown += new KeyEventHandler(onkeydown);
      this.KeyUp += new KeyEventHandler(onkeyup);

      Init(screen.Size);
    }

    //intialize rendering
    private void Init(Size size)
    {
      //setup rendering device
      buffersize = size;
      backbuffer = new Bitmap(buffersize.Width, buffersize.Height);
      graphics = Graphics.FromImage(backbuffer);

      timer.GetETime(); //reset timer
    }

    //main rendering function
    private void Render(Graphics g)
    {
      //clear back buffer
      graphics.Clear(Color.Black);
      graphics.ResetTransform();
      graphics.ScaleTransform(screenScale, -screenScale);
      graphics.TranslateTransform(buffersize.Width / 2.0f /
          screenScale, -buffersize.Height / 2.0f / screenScale);

      //draw to back buffer
      DrawScreen();

      //present back buffer
      g.DrawImage(backbuffer, new Rectangle(0, 0, buffersize.Width,
       buffersize.Height), 0, 0, buffersize.Width, 
       buffersize.Height, GraphicsUnit.Pixel);
    }

    //draw the screen
    private void DrawScreen()
    {
      //draw our simulation here
    }

    //process game logic
    private void DoFrame()
    {
      //get elapsed time since last frame
      float etime = timer.GetETime();

      //process input
      ProcessInput();

      ////////////////////////////////
      //integrate our simulation here
      ////////////////////////////////

      //redraw our screen
      screen.Invalidate();
    }


    //process keyboard input
    private void ProcessInput()
    {
      if (leftHeld)
        steering = -1;
      else if (rightHeld)
        steering = 1;
      else
        steering = 0;

      if (upHeld)
        throttle = 1;
      else
        throttle = 0;

      if (downHeld)
        brakes = 1;
      else
        brakes = 0;
    }

    private void onkeydown(object sender, KeyEventArgs e)
    {
      switch (e.KeyCode)
      {
        case Keys.Left:
          leftHeld = true;
          break;
        case Keys.Right:
          rightHeld = true;
          break;
        case Keys.Up:
          upHeld = true;
          break;
        case Keys.Down:
          downHeld = true;
          break;
        default: //no match found
          return; //return so handled dosnt get set
      }

      //match found
      e.Handled = true;
    }

    private void onkeyup(object sender, KeyEventArgs e)
    {
      switch (e.KeyCode)
      {
        case Keys.Left:
          leftHeld = false;
          break;
        case Keys.Right:
          rightHeld = false;
          break;
        case Keys.Up:
          upHeld = false;
          break;
        case Keys.Down:
          downHeld = false;
          break;
        default: //no match found
          return; //return so handled dosnt get set
      }

      //match found
      e.Handled = true;
    }

    //rendering - only when screen is invalidated
    private void screen_Paint(object sender, PaintEventArgs e)
    {
      Render(e.Graphics);
    }

    //when the os gives us time, run the game
    private void ApplicationIdle(object sender, EventArgs e)
    {
      // While the application is still idle, run frame routine.
      DoFrame();
    }

    private void MenuExit_Click(object sender, EventArgs e)
    {
      this.Close();
    }
  }

  //keep track of time between frames
  class Timer
  {
    //store last time sample
    private int lastTime = Environment.TickCount;
    private float etime;

    //calculate and return elapsed time since last call
    public float GetETime()
    {
      etime = (Environment.TickCount - lastTime) / 1000.0f;
      lastTime = Environment.TickCount;

      return etime;
    }
  }
}

Phase Two - Rigid Body Simulation

Ok, now we’re getting into some good stuff here. Let's put everything we just covered on the back burner now and talk about some physics. We’re going to be using a very simple Euler integration method. Basically, each frame we accumulate a bunch of forces (in our case from each wheel of the vehicle) and calculate the resultant acceleration, which is in the form of A=F/M (the same as F=MA, Newton’s second law of motion). We use this to modify Newton’s first law of motion, “an object in motion stays in motion…” So we calculate our A, and we integrate it into our V. Without an A, V would be constant, hence staying in motion, if no forces should act on it. Newton's third law gets applied in the form that any potential force the vehicle is applying to the ground, gets applied in the opposite direction to the vehicle (I'll explain this in the vehicle section). This topic is much easier to explain with symbols. So, P is our vehicle position, V is its linear velocity, F is the net force acting on it, M is its mass, A is the resultant acceleration, and T is the time step (the value our timer gave us from the last frame).

  A = F / M
  V = V + A * T
  P = P + V * T

This is a basic linear rigid body simulator. Each frame, we total up some F, integrate it, and then zero out F to restart the accumulation the next frame. Now let’s talk about rotation. The angular case is nearly identical to the linear case (especially in 2D). Instead of P we have an Angle, instead of V we have an Angular Velocity, instead of F we have a torque, and instead of M we have inertia. So the angular model looks like this

  AngA = Torque / Inertia
  AngV = AngV + AngA * T
  Angle = Angle + AngV * T

This is what that bit of code looks like. % is the cross product operator for my vector class.

public void AddForce(Vector worldForce, Vector worldOffset)
{
  //add linar force
  m_forces += worldForce;
  //and it's associated torque
  m_torque += worldOffset % worldForce;
}

Code Dump

Here is my rigid body object. You’ll notice all the properties I mentioned above. It has a Draw function which will draw its rectangle to the provided graphics object. It has an AddForce function, a space conversion method, to and from world space (very handy), and a function that returns the velocity of a point on the body (in world space). This point velocity is a combination of the linear velocity and the angular velocity. But the angular velocity is multiplied by the distance the point is from the center of rotation and perpendicular to its offset direction. So to kill two birds with one stone, I simply find the orthogonal vector to the point offset and multiply it by the angular velocity (then add the linear velocity.)

One thing you may be curious about is how I calculate the inertia value. That is a generalized formula I found at this link.

//our simulation object
class RigidBody
{
  //linear properties
  private Vector m_position = new Vector();
  private Vector m_velocity = new Vector();
  private Vector m_forces = new Vector();
  private float m_mass;

  //angular properties
  private float m_angle;
  private float m_angularVelocity;
  private float m_torque;
  private float m_inertia;

  //graphical properties
  private Vector m_halfSize = new Vector();
  Rectangle rect = new Rectangle();
  private Color m_color;

  public RigidBody()
  { 
    //set these defaults so we dont get divide by zeros
    m_mass = 1.0f; 
    m_inertia = 1.0f; 
  }

  //intialize out parameters
  public void Setup(Vector halfSize, float mass, Color color)
  {
    //store physical parameters
    m_halfSize = halfSize;
    m_mass = mass;
    m_color = color;
    m_inertia = (1.0f / 12.0f) * (halfSize.X * halfSize.X)
       * (halfSize.Y * halfSize.Y) * mass;

    //generate our viewable rectangle
    rect.X = (int)-m_halfSize.X;
    rect.Y = (int)-m_halfSize.Y;
    rect.Width = (int)(m_halfSize.X * 2.0f);
    rect.Height = (int)(m_halfSize.Y * 2.0f);
  }

  public void SetLocation(Vector position, float angle)
  {
    m_position = position;
    m_angle = angle;
  }

  public Vector GetPosition()
  {
    return m_position;
  }

  public void Update(float timeStep)
  {
    //integrate physics
    //linear
    Vector acceleration = m_forces / m_mass;
    m_velocity += acceleration * timeStep;
    m_position += m_velocity * timeStep;
    m_forces = new Vector(0,0); //clear forces

    //angular
    float angAcc = m_torque / m_inertia;
    m_angularVelocity += angAcc * timeStep;
    m_angle += m_angularVelocity * timeStep;
    m_torque = 0; //clear torque
  }

  public void Draw(Graphics graphics, Size buffersize)
  {
    //store transform, (like opengl's glPushMatrix())
    Matrix mat1 = graphics.Transform;

    //transform into position
    graphics.TranslateTransform(m_position.X, m_position.Y);
    graphics.RotateTransform(m_angle/(float)Math.PI * 180.0f);

    try
    {
      //draw body
      graphics.DrawRectangle(new Pen(m_color), rect);

      //draw line in the "forward direction"
      graphics.DrawLine(new Pen(Color.Yellow), 1, 0, 1, 5);
    }
    catch(OverflowException exc)
    {
      //physics overflow :(
    }  

    //restore transform
    graphics.Transform = mat1;
  }

  //take a relative vector and make it a world vector
  public Vector RelativeToWorld(Vector relative)
  {
    Matrix mat = new Matrix();
    PointF[] vectors = new PointF[1];

    vectors[0].X = relative.X;
    vectors[0].Y = relative.Y;

    mat.Rotate(m_angle / (float)Math.PI * 180.0f);
    mat.TransformVectors(vectors);

    return new Vector(vectors[0].X, vectors[0].Y);
  }

  //take a world vector and make it a relative vector
  public Vector WorldToRelative(Vector world)
  {
    Matrix mat = new Matrix();
    PointF[] vectors = new PointF[1];

    vectors[0].X = world.X;
    vectors[0].Y = world.Y;

    mat.Rotate(-m_angle / (float)Math.PI * 180.0f);
    mat.TransformVectors(vectors);

    return new Vector(vectors[0].X, vectors[0].Y);
  }

  //velocity of a point on body
  public Vector PointVel(Vector worldOffset)
  {
    Vector tangent = new Vector(-worldOffset.Y, worldOffset.X);
    return tangent * m_angularVelocity + m_velocity;
  }

  public void AddForce(Vector worldForce, Vector worldOffset)
  {
    //add linar force
    m_forces += worldForce;
    //and it's associated torque
    m_torque += worldOffset % worldForce;
  }
}

Testing

To make sure your rigid body works, instantiate one in your Init() function and apply a force with some offset in the DoFrame function. If you apply a constant worldOffset, the body will continue to accelerate its angular velocity. If you take your offset and run it through the RelativeToWorld function, the body will angularly accelerate in one direction and then come back the other way, like a pendulum as the point the force is applied to changes. Play around with this for a while, this has to work and make sense in order for the next section to work.

Phase Three - The Vehicle

Assuming everything has gone well above, you should have a rigid body actor in your scene that you can apply forces to and watch move around. Now all that’s left is to calculate these forces in a way that will simulate a vehicle. For that we are going to need a vehicle object. I recommend deriving directly from you rigid body object since the chassis is essentially a rigid body. In addition to that we will need to construct a “wheel" object. This wheel will handle the steering direction of each wheel, the velocity the wheel is spinning, and calculate the forces that that particular wheel applies to the chassis (all in vehicle space). Since our wheel is known to be constrained to the vehicle, we don’t need to simulate it as another rigid body (though you could, but not in the 2D case.) We will simply duplicate the angular properties of the rigid body in the wheel object.

So we’ll need: Wheel Velocity, Wheel Inertia, and and Wheel Torque. We’ll also need the relative offset of the wheel in the vehicle space, and the angle the wheel is facing (this is constant for the back wheels, unless you want 4 wheel steering.) Just like the rigid body, the wheel's torque function acts as an accumulator, we add torques to it and after it gets integrated the torque is zeroed out. The AddTorque function is where you will apply a wheel torque from either the transmission (to make you go) or from the brakes (to make you stop). Internally the wheel will generate a torque caused by the friction on the road.

The wheel object also needs a SetSteering function. This function calculates two vectors: an effective Side Direction, and an effective Forward Direction (both in vehicle space) that the tire patch will act on. The force applied on the tire by the ground acting in the side direction will directly translate into the chassis. Meanwhile the force acting in the forward direction will not only act on the chassis, but it will induce a rotation of the tire. Here is the SetSteering function; you will see I used the Drawing2D.Matrix to transform the initial forward and side vectors by the steering angle (I had to convert the vectors to “points" in order to transform them by the matrix.)

public void SetSteeringAngle(float newAngle)
{
  Matrix mat = new Matrix();
  PointF[] vectors = new PointF[2];

  //foward vector
  vectors[0].X = 0;
  vectors[0].Y = 1;
  //side vector
  vectors[1].X = -1;
  vectors[1].Y = 0;

  mat.Rotate(newAngle / (float)Math.PI * 180.0f);
  mat.TransformVectors(vectors);

  m_forwardAxis = new Vector(vectors[0].X, vectors[0].Y);
  m_sideAxis = new Vector(vectors[1].X, vectors[1].Y);
}

Force Calculation

So, if the vehicle is sitting there not moving with its front wheels turned, and you push it, a force will be generated in the opposite direction you push. This force gets projected onto these two directions. If the wheels were straight there would be no side force. So the vehicle would simply roll forward. But since the wheels are turned, there is a bit of the force that acts in the “effective side direction” so we apply an opposite force to the chassis. This is what causes you to turn when you steer the wheels. To get this force that gets projected onto the two directions, we need to first determine the velocity difference between the tire patch and the road. If the wheel is spinning at the same speed the ground is wizzing by, then there is effectively no force acting on the vehicle. But as soon as you slam on the brakes and stop the wheel, there is a huge velocity difference and this is what causes the force that stops your car.

So here is the process broken down into 6 steps, for each wheel.

Step 1, calculate the effective direction vectors (with steering function).

Step 2, calculate velocity difference. The ground speed is determined via the “PointVel” function on the rigidbody, given the current wheel’s world offset.

Step 3, project this velocity onto the two effective directions.

Step 4, generate an equal and opposite force for the two direction and call this the “response force”. This is what gets added to the chassis for each wheel.

Step 5, calculate the torque that the forward response force created on the wheel, and add this to the wheel torque.

Step 6, integrate the wheel torques into the wheel velocity.

That bit of code looks like this:

public Vector CalculateForce(Vector relativeGroundSpeed, float timeStep)
{
  //calculate speed of tire patch at ground
  Vector patchSpeed = -m_forwardAxis * m_wheelSpeed *
    m_wheelRadius;

  //get velocity difference between ground and patch
  Vector velDifference = relativeGroundSpeed + patchSpeed;

  //project ground speed onto side axis
  float forwardMag = 0;
  Vector sideVel = velDifference.Project(m_sideAxis);
  Vector forwardVel = velDifference.Project(m_forwardAxis, out forwardMag);

  //calculate super fake friction forces
  //calculate response force
  Vector responseForce = -sideVel * 2.0f;
  responseForce -= forwardVel;

  //calculate torque on wheel
  m_wheelTorque += forwardMag * m_wheelRadius;

  //integrate total torque into wheel
  m_wheelSpeed += m_wheelTorque / m_wheelInertia * timeStep;

  //clear our transmission torque accumulator
  m_wheelTorque = 0;

  //return force acting on body
  return responseForce;
}

Almost Done!

We’re in the home stretch here now. Now we have a way to calculate the force each wheel generates on the chassis. Every frame, all we have to do is set our transmission and brake torques, our steering angle, calculate each wheel force, add these to the chassis, and integrate the rigid body. Badaboom badabing, vehicle done! :)

Conclusion

Here is the entire source code for the project. If you have any questions or comments please feel free to post them here and either I can make things more clear or maybe someone else could offer some better expertise. If you’d like you can email me at Kincaid05 on google's fine emailing service.

Thanks for reading and I hope this was informative.
-Matt Kincaid

• The objects that the algorithm is dealing with are convex³ 3D objects only. It is true that the algorithm will work for some concave⁴ objects but care must be taken⁵.
• Detection of any pair of faces colliding is a well-defined and solved problem.
• Both objects have a pre-calculated center point.
• The accuracy constant k is set to be (min/max estimates)
• Object A has exactly n faces while object B has m faces.

• Optimization: check collision between object A's bounding sphere and object B's bounding sphere.
• Find the closest k faces of object A relative to object B.
• Calculate the geometric center of the new selection and the bounding sphere radius.
• Find the closest k faces of object B relative to object A's new selection of k faces.
• Calculate the geometric center of object B's new selection of faces and their maximal bounding sphere radius.
• Optimization: check collision between spheres to determine if there is even a chance for face collisions.
• Sort the two sets of faces by increasing distance (optional, might be replaces by a good insertion algorithm).
• Test the two sets of faces against each other, starting with the closest pairs of faces.

• Object B can be considered as a point object located at its center (object centers have been shown in step 1).

• Using both object's positions, calculate the relative position vector. Assuming is the center of object A (x,y,z - a 3D vector) and the center of object B, the vector connecting the center points is then:
• Start looping through all faces of object A, find each face center relative to object A's center. Total time complexity is O(n). It takes less operations to calculate the faces centers than to do it in the beginning and apply transformations on them.
• For each face, find the vector connecting its center with object B's center, using the simple relation: , where is the vector connecting the center of face i with the center of object B, the relative positions of the two object's centers and is the position of the face relative to object A's center (the object it belongs to). Of course, i = 0,1,2,3...,n.
• For each , find the vector's size: and store into selection*. Of course, the size itself is irrelevant, as is the direction of the vector. The only important quantity is the vector's size squared. Note that there is no need to take the square root because if a certain face's squared vector size is larger than another face's squared vector size, the same holds for the actual sizes themselves.
• Depending on the insertion technique used to insert the distances (squared sizes of the relative position vectors), find the smallest k distances which should be straightforward if the insertion was done in an efficient way, and remember their corresponding faces.

• Calculate the relative position vector between object B's center and the center of the new selection of k faces found earlier. Assuming is the center of object B and the center of the selection in object A, the vector connecting the center points is then:
• Start looping through all faces of object B, find each face center relative to the center of the selection in object A. Total time complexity is O(m). It takes less operations to calculate the faces centers than to do it in the beginning and apply transformations on them.
• For each face, find the vector connecting its center with the center of the selection in object A, using the simple relation: , where is the vector connecting the center of face i with the center of the selection in object A, the relative positions of the center of object B and the center of the selection in object A and is the position of the face relative to object B's center (the object it belongs to). Of course, i = 0,1,2,3...,m.
• For each , find the vector's size: and store into selection*.
• Depending on the insertion technique used to insert the distances (squared sizes of the relative position vectors), find the smallest k distances which should be straightforward if the insertion was done in an efficient way, and remember their corresponding faces.

• The algorithm supports detection of mostly convex objects. Some concave objects will work as well but there are bound to be degenerate cases that cause the algorithm to fail. It hasn't been proven mathematically that there are any cases that will cause failure but assuming there aren't is not a good idea.
• Moving objects in computer simulations are by nature problematic because time in computer systems is discrete. In numerical approximations, time is discrete and therefore object positions are calculated in time steps and are really "teleporting" in small steps to create the illusion of motion. The problem arises when the velocity of objects is very large and/or objects are very small. It is quite possible that an object will be close to another object while not colliding, yet at the next time step will be half intersecting with the other object. The problem is then what faces to test. After all, the "closest" faces the object previously had are now inside the other close object. Even more so, it's possible that if the velocity is large enough, the object might pass right through the other object without any way of us detecting this.
• In a similar fashion, rotating objects pose another problem. If objects have an angular velocity (or momentum, whichever you prefer) it's possible that before the collision some faces are the closest selection while after a small time step, other faces should be the closest selection. For example, consider a normal cube floating above a table at a small height. The closest faces to the table are naturally the two faces (triangles) that comprise the base of the cube. If the cube is about to rotate 45 degrees in one time step (very fast rotation) it wouldn't be correct checking the base of the cube against the table.

• Find the closest object in collection A relative to the entire collection B. This implies each object in A has to have a center (this is mandatory for OFG anyway) and the entire collection has to have a center as well (averaging the centers of the objects, should be a part of initialization).
• Find the closest object in collection B relative to the object found earlier in collection A.
• Feed the two objects found to the OFG algorithm.

• Calculate the translation of object A after a time unit elapses.
• Feed the two objects into the OFG algorithm, forgoing step 6 completely and without the bounding sphere calculations in steps 3 and 5. Also, there is no need to perform step 8 (the final step) of the OFG algorithm just yet.
• Connect each vertex in selection A with the same vertex on the virtual selection of A that is translated by the amount calculated in step 1. This gives a vector that represents the edge to test against object B.
• Test the newly created edges against the selection in object B found by the OFG algorithm. They are the most likely to be intersecting object B.
• Test the original faces in selection A against the selection in B, as was intended in step 8 of the OFG algorithm.
• Test the translated virtual selection of A against the selection of B using step 8 of the OFG algorithm.

• Using the translation found in step 1 earlier and given that it's easily possible to calculate the relative velocity of the two objects, it's possible to estimate the time of collision (therefore obtaining the delta time needed for the objects to collide).
• Moving half of the distance (half a time unit) and checking for collisions. If collision occurred, great. If not, move another half, so on and so forth. This has a time complexity of log(n) times of performing all the collision tests.

• Solve the plane equation of the second face for each vertex in the first face. If all vertices produce the same sign, the first face definitely doesn't pierce the plane in question. Creating the plane equation is easy given two vectors on that plane (three vertices of one face can give two vectors). If any vertex gives the other sign, there is a possibility of collision.

• Calculate the intersection points of the edges found piercing the plane in stage one using the second face's plane equation.
• For each edge in the second face, take the cross product of the edge vector with the normal to the face, giving a perpendicular plane.
• The points must lie either on the plane or on either side. Using the same logic as in stage one can determine this for each edge plane.
• If the point lies in the same part of space ("side") for all planes corresponding to the edges of the face, the point is within the face and therefore a collision has occurred
• If both points lie outside the planes of the edges (not all signs are the same for each point checked against the edge planes), there is definitely no collision.

• Use the 1024 new samples taken in a[n] and b[n] to compute the instant energy 'e', using the following formula (i0 is the position of the 1024 samples to process): /reference/programming/features/beatdetection/formula1.png
  (R1)
• Compute the local average energy '' on the 44100 samples of a history buffer (B): /reference/programming/features/beatdetection/formula2.png
  (R2)
• Shift the 44032 history buffer (B) of 1024 indexes to the right so that we make room for the 1024 new samples and evacuate the oldest 1024 samples.
• Move the 1024 new samples on top of the history buffer.
• Compare 'e' to 'C * ' where C is a constant which will determine the sensibility of the algorithm to beats. A good value for this constant is 1.3. If 'e' is superior to 'C * ' then we have a beat!

• Compute the instant sound energy 'e' on the 1024 new sample values taken in (an) and (bn) using the formula (R1)
• Compute the average local energy with (E) sound energy history buffer: /reference/programming/features/beatdetection/formula3.png
  (R3)
• Shift the sound energy history buffer (E) of 1 index to the right. We make room for the new energy value and flush the oldest.
• Pile in the new energy value 'e' at E[0].
• Compare 'e' to 'C*'.

• Translation
• Scale
• Rotation

Patrick Borthelow, "Sensing Tire Slip Angles At the Racetrack," Sensors, September 1991.

Carrol Smith, Engineer to Win, Prepare to Win, Build to Win, from Classic Motorbooks, P.O. Box 1/RT021, Osceola, WI, 54020.

concerns adhesion, grip angle, and elastic modeling. This has been covered in Parts 2, 4, 7, and 10, and will be covered in several later parts.

concerns handling, suspension movement, and motion of a car around a course; has been covered in Parts 1, 4, 5, and 8 and will continue.

concerns modeling of engine performance and acceleration. Has been covered in Parts 3, 6, and 9 and will also continue.

concerns the design of a working program that captures all the physics. This is the ultimate goal of the series. It was begun in Part 8 and will eventually dominate discussion.

• Weight Transfer
• Keeping Your Tires Stuck to the Ground
• Basic Calculations
• There is No Such Thing as Centrifugal Force
• Introduction to the Racing Line
• Speed and Horsepower
• The Circle of Traction
• Simulating Car Dynamics with a Computer Program
• Straights
• Grip Angle

presenting a detailed model of suspension movement (suggested by Bob Mosso)

presenting the dynamics of entering and leaving corners, chicanes, and slaloms (this one suggested by Karen Babb)

explaining why spins and other losses of control occur

exploring what, exactly, is meant by smoothness

in a computer program; in several articles

in a computer program; also in several articles

(new-race-car)

(define car-161 (new-race-car))

(define car-161 (new-race-car)) 
(define car-240 (new-race-car))
(define car-70  (new-race-car))

(race-car-x car-161)
(race-car-y car-161)
(race-car-z car-161)

(race-car-vx car-161)
(race-car-vy car-161)
(race-car-vz car-161)

(define x1 (race-car-x car-161))
(define y1 (race-car-y car-161))
(define z1 (race-car-z car-161))

(define x2 (race-car-x car-161))
(define y2 (race-car-y car-161))
(define z2 (race-car-z car-161))

(define vx (/ (- x2 x1) (- t2 t1)))
(define vy (/ (- y2 y1) (- t2 t1)))
(define vz (/ (- z2 z1) (- t2 t1)))

(/ (- x2 x1) (- t2 t1))

(x2 - x1) / (t2 - t1)	

(/ ...)	

(- ...) (- ...)	

(/ (- ...) (- ...))

(/ (- x2 x1) (- t2 t1))

(set-race-car-vx! car-161 vx)
(set-race-car-vy! car-161 vy)
(set-race-car-vz! car-161 vz)