VB:Tutorials:Building A Simple Physics Engine

From GPWiki

Tutorial by Nicholas Gorski [GMan]

This article will show you the basics of making a simple physics engine. This engine will handle circle-circle collision detection. Before we start, it may be useful to read: Pool Hall Lessons[1], which is what this is based on.

1 Creating the Engine Core
2 Collision Detection
3 Collision Response
4 Almost Done
5 Done

[edit]

Creating the Engine Core

Open a new ActiveX DLL and rename the class PEngine. This is the main engine which handles everything. Now add a new module named modMain.

A common thing we need is gravity, so make two public variables called sGravityX and sGravityY, both as Single. Also we need a way to make sure our physics loop doesn’t freeze, so add a Double called dMaxTime. Place them inside modMain:

Public sGravityX As Single
Public sGravityY As Single
 
Public dMaxTime As Double

Just for this engine I will provide rudimentary world bounding box collision detection. In a real physics engine, you should use proper line collision detection: [2]. Part Two will implement that.

So add these to modMain:

Public sWorldLeft As Single
Public sWorldRight As Single
Public sWorldTop As Single
Public sWorldBottom As Single

So what do we have now? Nothing. We need a way to store circles, so make a class called PCircle. This is where a lot of code will go. We can make two public variables in PCircle named Radius to hold the radius, and Mass to hold Mass. We also need to be able to store position, velocity, and acceleration. Unfortunately, Visual Basic doesn’t let us make user-defined types public, so we have to make some functions to access those variables. Also, we need a way to store vectors, so make the PVector type. Add this code to the PCircle class:

(Just a side note here. Add a reference to DirectX to your application, use D3DVECTOR instead of 'PVector' and you're set to deal with 3D physics.. why was this done in 2D anyway?)

Public Type PVector
    X As Single
    Y As Single
End Type
 
Private Position As PVECTOR
Private Velocity As PVECTOR
Private Acceleration As PVECTOR
 
Public Restitution As Single
Public Radius As Single
 
Public Function SetPosition(ByVal X As Single, ByVal Y As Single)
Position.X = X
Position.Y = Y
End Function
 
Public Function GetPosition() As PVECTOR
GetPosition.X = Position.X
GetPosition.Y = Position.Y
End Function
 
Public Function ApplyPosition(ByVal X As Single, ByVal Y As Single)
Position.X = Position.X + X
Position.Y = Position.Y + Y
End Function
 
Public Function SetVelocity(ByVal X As Single, ByVal Y As Single)
Velocity.X = X
Velocity.Y = Y
End Function
 
Public Function GetVelocity() As PVECTOR
GetVelocity.X = Velocity.X
GetVelocity.Y = Velocity.Y
End Function
 
Public Function ApplyVelocity(ByVal X As Single, ByVal Y As Single)
Velocity.X = Velocity.X + X
Velocity.Y = Velocity.Y + Y
End Function
 
Public Function SetAcceleration(ByVal X As Single, ByVal Y As Single)
Acceleration.X = X
Acceleration.Y = Y
End Function
 
Public Function GetAcceleration() As PVECTOR
GetAcceleration.X = Acceleration.X
GetAcceleration.Y = Acceleration.Y
End Function
 
Public Function ApplyAcceleration(ByVal X As Single, ByVal Y As Single)
Acceleration.X = Acceleration.X + X
Acceleration.Y = Acceleration.Y + Y
End Function

The Apply functions are useful for things like gravity, where you want to add to the existing variables. Restitution is how much energy something keeps after a collision. So putting 0.9 for Restitution means the the object keeps 90% of its energy when it hits something. Now we get to make Visual-Basic-style pointers. Actually, Visual Basic has pointer capabilities, without Windows API, but that's not for this article.

Go to the PEngine class. Make a variable array called Circles, along with a counter saying how many circles there are, like this:

Private Circles() As PCircle
Private NumCircles As Long

Make a function to allow adding circles:

Public Function AddCircle(ByRef CircleData As PCircle)
ReDim Preserve Circles(NumCircles) As PCircle
Set Circles(NumCircles) = CircleData
NumCircles = NumCircles + 1
End Function

There. We're getting closer to a functional engine. What this function does is set the Circles(NumCircles) variable to point to the source of CircleData. So changing CircleData changes Circle(NumCircles) and vice versa.

Remember those gravity variables? We need to access them, so make a public function called SetWorldProperties:

Public Function SetWorldProperties(ByVal GravityX As Single, ByVal GravityY As Single, _ 
    ByVal MaxTime As Double, _
    ByVal WorldLeft As Single, ByVal WorldRight As Single, _ 
    ByVal WorldTop As Single, ByVal WorldBottom As Single)
sGravityX = GravityX
sGravityY = GravityY
dMaxTime = MaxTime
sWorldLeft = WorldLeft
sWorldRight = WorldRight
sWorldTop = WorldTop
sWorldBottom = WorldBottom
End Function

Closer! Now we need a way to update the position, velocity, and acceleration of the circles. However, I prefer time-based modeling, so add a class called PTimer and add this code:

Private Declare Function QueryPerformanceFrequency Lib "kernel32" (lpFrequency As Currency) As Long
Private Declare Function QueryPerformanceCounter Lib "kernel32" (lpPerformanceCount As Currency) As Long
 
Private Freq As Currency
Private StartTime As Currency
Private EndTime As Currency
Private TimeElapse As Double
Public TimeFactor As Double
 
Public Function Timing(ByVal Action As Boolean)
On Error Resume Next 'When the timer has not been initialize, but the timer
                     'is stopped, a divide by zero error occurs
Select Case Action
    Case True 'Start
        QueryPerformanceFrequency Freq
        QueryPerformanceCounter StartTime
    Case False 'Stop
        QueryPerformanceCounter EndTime
        TimeElapse = (EndTime - StartTime) / Freq
End Select
End Function
 
Public Property Get TimeElapsed() As Double
TimeElapsed = TimeElapse * TimeFactor
End Property

The TimeFactor variable allows you to go in slow-motion, or fast-motion, by setting it higher or lower than 1, while 1 runs at normal time. This is one of those things you have to tweak with. Add an object called Timer at the very top of PEngine. Don’t forget the ‘New’ keyword. The top looks like this:

Private Timer As New PTimer
 
Private Circles() As PCircle
Private NumCircles As Long

Back to circles, we need a function to update circles. This will have two parameters. One is called TimeElapsed, to allow time-based modeling. If you don't want time-based modeling, just pass 1 to this function. The second is the index to the circle object:

(Sorry to interrupt this well written text here. The segment below is erratical: The increase in Y velocity should be dependant on TimeElapsed, and the new new Y position should be calculated based on (OldVelocity * TimeElapsed) + ((NewVelocity - OldVelocity) * 0.5 * TimeElapsed). Also, ApplyAcceleration and then SetAcceleration 0,0 seems a bit unnecessary. I leave to the original author to make the necessary fixes. /Waterman)

(Actually, if you made the changes or additions, people would know what you're talking about. When your suggestions are applied to this engine, the objects go UP. There is also no referrence to mass anywhere here. How do we deal with objects with different masses)

Private Function UpdateCircle(ByVal TimeElapsed As Double, ByVal CircleIndex As Long)
With Circles(CircleIndex)
    .ApplyAcceleration sGravityX, sGravityY 'Add Gravity
    .ApplyVelocity .GetAcceleration.X, .GetAcceleration.Y 'Update Acceleration
    .ApplyPosition .GetVelocity.X * TimeElapsed, .GetVelocity.Y * TimeElapsed 'Update Position
    .SetAcceleration 0, 0 'Reset Acceleration
End With

Just one last function! The UpdateWorld function, which updates all the circles. It too takes a TimeElapsed variable:

Public Function UpdateWorld(ByVal TimeElapsed As Double)
Dim i As Long
 
For i = 0 To NumCircles - 1
    Call UpdateCircle(TimeElapsed, i) 'Update Circle
Next i
End Function

Now what do we have? A basic particle simulator. Also note that we didn’t use the Timer object. Why is this a basic particle simulator? None of the circles collide! So, whip out the 'Pool Hall Lessons' web page and follow along!(See Top)

[edit]

Collision Detection

After opening 'Pool Hall Lessons' [3], skip to the working code. Notice that we too need the vector functions. The first one you see is the Magnitude() function. So in modMain, type:

Public Function VectorMagnitude(ByRef A As PVector) As Single
VectorMagnitude = Sqr(A.X * A.X + A.Y * A.Y)
End Function

The second vector function you see is the Normalize() function, so add:

Public Function VectorNormalize(ByRef A As PVector)
Dim Mag As Single
 
Mag = Sqr(A.X * A.X + A.Y * A.Y)
 
If Mag <> 0 Then
    A.X = A.X / Mag
    A.Y = A.Y / Mag
End If
End Function

The third one you see is the Dot Product:

Public Function VectorDotProduct(ByRef A As PVector, ByRef B As PVector) As Single
VectorDotProduct = A.X * B.X + A.Y * B.Y
End Function

The last one is the Distance function:

Public Function VectorDistance(ByRef A As PVector, ByRef B As PVector) As Single
VectorDistance = Sqr((A.X - B.X) * (A.X - B.X) + (A.Y - B.Y) * (A.Y - B.Y))
End Function

Now to get started. Define the function in PEngine like this:

Public Function CircleCircleCollision(ByVal TimeElapsed As Double, _
    ByVal Index1 As Long, ByVal Index2 As Long) As Boolean

At the top of the function, type this:

Dim A As PVector, B As PVector, VA As PVector, VB As PVector
Dim MoveVec As PVector
A = Circles(Index1).GetPosition
B = Circles(Index2).GetPosition
VA = Circles(Index1).GetVelocity
VB = Circles(Index2).GetVelocity
VA.X = VA.X * TimeElapsed
VA.Y = VA.Y * TimeElapsed
VB.X = VB.X * TimeElapsed
VB.Y = VB.Y * TimeElapsed
 
MoveVec.X = VA.X - VB.X
MoveVec.Y = VA.Y - VB.Y

This function first stores the Positions and Velocities of both circles in variables with shorter names. Not only is it easier to type "A" than "Circles(Index1).GetPosition", it's faster, because the latter calls a function. The next thing the code does is it subtracts one velocity from the other. Because the test is orginally for circles where one is moving and the other is stationary, we must find their relative velocities. Now to code. The first section in 'Pool Hall Lessons' is:

// Early Escape test: if the length of the movevec is less
// than distance between the centers of these circles minus
// their radii, there's no way they can hit. 
double dist = B.center.distance(A.center);
double sumRadii = (B.radius + A.radius);
dist -= sumRadii;
if(movevec.Magnitude() < dist){
return false;
}

In Visual Basic, this is:

Dim dist As Single, sumRadii As Single
dist = VectorDistance(A, B)
sumRadii = Circles(Index1).Radius + Circles(Index2).Radius
dist = dist - sumRadii
If VectorMagnitude(MoveVec) < dist Then
    CircleCircleCollision = False
    Exit Function
End If

The next section is:

// Normalize the movevec
Vector N = movevec.copy();
N.normalize();
 
// Find C, the vector from the center of the moving 
// circle A to the center of B
Vector C = B.center.minus(A.center);
 
// D = N . C = ||C|| * cos(angle between N and C)
double D = N.dot(C);

This just initializes some variables. It does not do any tests. In Visual Basic this is:

Dim N As PVector, C As PVector, D As Single
N = MoveVec
Call VectorNormalize(N)
 
C.X = B.X - A.X
C.Y = B.Y - A.Y
 
D = VectorDotProduct(N, C)

Now this new information is put to use. This section makes sure the ball is moving towards the other ball:

// Another early escape: Make sure that A is moving 
// towards B! If the dot product between the movevec and 
// B.center - A.center is less that or equal to 0, 
// A isn't isn't moving towards B
if(D <= 0){
  return false;
}

That last part is the code, for those who don't know C++. Easy:

If D <= 0 Then
    CircleCircleCollision = False
    Exit Function
End If

Another test makes sure they can even hit each other:

// Find the length of the vector C
double lengthC = C.Magnitude();
 
double F = (lengthC * lengthC) - (D * D);
 
// Escape test: if the closest that A will get to B 
// is more than the sum of their radii, there's no 
// way they are going collide
double sumRadiiSquared = sumRadii * sumRadii;
if(F >= sumRadiiSquared){
  return false;
}

This is fairly simple to implement in Visual Basic:

Dim LengthC As Single, F As Single, sumRadiiSquared As Single
LengthC = VectorMagnitude(C)
F = (LengthC * LengthC) - (D * D)
sumRadiiSquared = sumRadii * sumRadii
 
If F >= sumRadiiSquared Then
    CircleCircleCollision = False
    Exit Function
End If

The next section gets a variable T and sees if it is less than 0:

// We now have F and sumRadii, two sides of a right triangle. 
// Use these to find the third side, sqrt(T)
double T = sumRadiiSquared - F;
 
// If there is no such right triangle with sides length of 
// sumRadii and sqrt(f), T will probably be less than 0. 
// Better to check now than perform a square root of a 
// negative number. 
if(T < 0){
  return false;
}

Also very simple in Visual Basic:

Dim T As Single
T = sumRadiiSquared - F
If T < 0 Then
    CircleCircleCollision = False
    Exit Function
End If

This last bit makes two variables and tests them:

// Therefore the distance the circle has to travel along 
// movevec is D - sqrt(T)
double distance = D - sqrt(T);
 
// Get the magnitude of the movement vector
double mag = movevec.Magnitude();
 
// Finally, make sure that the distance A has to move 
// to touch B is not greater than the magnitude of the 
// movement vector. 
if(mag < distance){
  return false;
}

In Visual Basic:

Dim Distance As Single, Mag As Single
Distance = D - Sqr(T)
Mag = VectorMagnitude(MoveVec)
 
If Mag < Distance Then
    CircleCircleCollision = False
    Exit Function
End If

Finallly! Now that we are sure the circles collide, move them with all this data we have until they just touch. Note that the article only shows code for a stationary-moving circle test, so I made code for two moving circles:

Dim UA As Single, UB As Single
 
If VectorMagnitude(VA) = 0 Then
    UA = 0
Else
    UA = VectorMagnitude(MoveVec) / VectorMagnitude(VA)
End If
If VectorMagnitude(VB) = 0 Then
    UB = 0
Else
    UB = VectorMagnitude(MoveVec) / VectorMagnitude(VB)
End If
 
VA.X = VA.X * UA
VA.Y = VA.Y * UA
VB.X = VB.X * UB
VB.Y = VB.Y * UB
 
Circles(Index1).ApplyPosition VA.X, VA.Y
Circles(Index2).ApplyPosition VB.X, VB.Y
 
VA = Circles(Index1).GetVelocity
VB = Circles(Index2).GetVelocity
A = Circles(Index1).GetPosition
B = Circles(Index2).GetPosition

This makes the circle just touch, then restores the position and velocity variables. Now that we have collision detection, time for collision response.

[edit]

Collision Response

Getting our circles to move in a realistic manner is easy. I will simply transfer the code from Pool Hall Lessons directly to Visual Basic. In C++ the code looks like this:

// First, find the normalized vector n from the center of 
// circle1 to the center of circle2
Vector n = circle1.center - circle2.center;
n.normalize(); 
// Find the length of the component of each of the movement
// vectors along n. 
// a1 = v1 . n
// a2 = v2 . n
float a1 = v1.dot(n);
float a2 = v2.dot(n);
 
// Using the optimized version, 
// optimizedP =  2(a1 - a2)
//              -----------
//                m1 + m2
float optimizedP = (2.0 * (a1 - a2)) / (circle1.mass + circle2.mass);
 
// Calculate v1', the new movement vector of circle1
// v1' = v1 - optimizedP * m2 * n
Vector v1' = v1 - optimizedP * circle2.mass * n;
 
// Calculate v1', the new movement vector of circle1
// v2' = v2 + optimizedP * m1 * n
Vector v2' = v2 + optimizedP * circle1.mass * n;
 
circle1.setMovementVector(v1');
circle2.setMovementVector(v2');

In Visual Basic, the code looks much easier:

N.X = A.X - B.X 'NOTE:
N.Y = A.Y - B.Y 'N already exists above
Call VectorNormalize(N)
 
Dim A1 As Single, A2 As Single
A1 = VectorDotProduct(VA, N)
A2 = VectorDotProduct(VB, N)
 
Dim optimizedP As Single
optimizedP = (2 * (A1 - A2)) / (Circles(Index1).Mass + Circles(Index2).Mass)
 
Dim V1 As PVector, V2 As PVector
V1.X = VA.X - (optimizedP * Circles(Index2).Mass * N.X)
V1.Y = VA.Y - (optimizedP * Circles(Index2).Mass * N.Y)
V2.X = VB.X + (optimizedP * Circles(Index1).Mass * N.X)
V2.Y = VB.Y + (optimizedP * Circles(Index1).Mass * N.Y)
 
Circles(Index1).SetVelocity V1.X * Circles(Index1).Restitution, V1.Y * Circles(Index1).Restitution
Circles(Index2).SetVelocity V2.X * Circles(Index2).Restitution, V2.Y * Circles(Index2).Restitution
 
CircleCircleCollision = True

This code, unlike 'Pool Hall Lessons', takes into account the restitution of the ball. Without it, the circles would just keep bouncing forever. Do NOT set the restitution above 1. Doing so would not only be unrealistic, it generates an overflow error.

[edit]

Almost Done

Now we just need to update the UpdateWorld function, which has taken quite a change:

Dim i As Long, j As Long, bCollided As Boolean, T As Double
 
If TimeElapsed = 0 Then Exit Function
 
Timer.TimeFactor = 1
Do 'Loop until...
    ''''''''''Circle - Circle''''''''''
    Timer.Timing True
    bCollided = False
    For i = 0 To NumCircles - 1
        For j = i To NumCircles - 1
            If i <> j Then
                If CircleCircleCollision(TimeElapsed, i, j) Then
                    bCollided = True
                End If
            End If
        Next j
    Next i
    '''''''''Circle - Line''''''''''''''
    Dim tX As Single, tY As Single, vX As Single, vY As Single
    For i = 0 To NumCircles - 1
        With Circles(i)
            tX = .GetPosition.X: tY = .GetPosition.Y
            vX = .GetVelocity.X * .Restitution: vY = .GetVelocity.Y * .Restitution
            If tX < sWorldLeft + .Radius Then
                bCollided = True
                .SetPosition sWorldLeft + .Radius, tY
                .SetVelocity -vX, vY
            ElseIf tX + .Radius > sWorldRight Then
                bCollided = True
                .SetPosition sWorldRight - .Radius, tY
                .SetVelocity -vX, vY
            ElseIf tY < sWorldTop + .Radius Then
                bCollided = True
                .SetPosition tX, sWorldTop + .Radius
                .SetVelocity vX, -vY
            ElseIf tY + .Radius > sWorldBottom Then
                bCollided = True
                .SetPosition tX, sWorldBottom - .Radius
                .SetVelocity vX, -vY
            End If
        End With
    Next i
    ''''''''''Update Timer''''''''''''''
    Timer.Timing False
    T = T + Timer.TimeElapsed
Loop Until (bCollided = False) Or (T > dMaxTime) 'no collisions are found, or the time has taken too long.
Timer.Timing True
 
'Now update circles
For i = 0 To NumCircles - 1
    Call UpdateCircle(TimeElapsed, i) 'Update Circle
Next i
 
Timer.Timing False
UpdateWorld = T + Timer.TimeElapsed

Now it makes sure the loop doesn’t freeze up. Also note that the UpdateWorld function returns a Double, telling how long the function took (in milliseconds I think). If it is in milliseconds, it is very fast. If it's in seconds, blame it on Visual Basic.

[edit]