//I make no claims to the validity or accuracy of this program. Use at your own risk.
//Don’t blame / sue me if it ends up causing you or your system damage in some way.
//
//*****************************************************************************************************
//Name : Simple2DCollisionsWithFriction.c
//Author: Rahil Baber
//Date : 3rd March 2011
//
//This program simulates a single 2D object colliding against the floor, with friction.
//Note that this program does not handle the object at rest properly.
//
//During the simulation press any key to restart the simulation.
//
//The shape of the object which we are colliding and its properties are set in the function
//InitializeSimulation(). Modify the variables in this function to create your own simulations.
//
//This program was written for the website EuclideanSpace check out
//http://euclideanspace.com/physics/dynamics/collision/practical/RahilBaberCorrectionToBrianMirtich.pdf
//for a description of the algorithm.
//
//The program uses OpenGL to render the simulation and GLUT to create and manage the windows.
//
//To compile the program in Visual Studio you may need to add the libraries
//opengl32.lib glu32.lib glut32.lib
//to "Project>Properties>Configuration Properties>Linker>Input>Additional Dependencies".
//
//To compile the program using gcc type "gcc Simple2DCollisionsWithFriction.c -lglut -lGLU -lGL -lm"
//at the command line.
//
//If you've found this program useful or interesting or are creating your own physics simulation
//application based on any part of it, then drop me an e-mail (see the pdf for the address).
//*****************************************************************************************************
#include
#include
#include
#define PI 3.141592653589793
//The time in milliseconds between frames of the animation.
#define TIMERMILLISECS 33
//*** The global variables ***
//Used to measure time (in milliseconds) to create a consistent animation speed.
int prevTime;
//The width and height of the simulation.
//These values are set when the window is resized (see the Resize() function).
double simulationWidth; //=(window width)/(window height)
double simulationHeight; //=1.0
//The scene is rendered as a line drawing, lineThickness indicates the thickness of those lines.
double lineThickness = 0.01;
//The remaining global variables describe the simulation and are initialized in InitializeSimulation().
double groundHeight; //The height of the ground.
double g; //The acceleration due to gravity.
double cF; //The coefficient of friction.
double cR; //The coefficient of restitution.
double I; //The moment of inertia of the object.
double cx, cy; //The position of the object.
double vx, vy; //The velocity of the object.
double a; //The angle the object is orientated at.
double w; //The angular velocity of the object.
double rx, ry; //The relative position of the point of contact during a collision.
//The following variables describe the shape of the object. The shape is described by an array of
//coordinates (boundary[0][] to boundary[noOfPoints-1][]) which are joined by line segments. The point
//boundary[0][] is joined to boundary[1][], which in turn is joined to boundary[2][], etc. We take the
//boundary to be a closed loop and so we join boundary[noOfPoints-1][] to boundary[0][].
int noOfPoints; //The number of points used to describe the boundary of the object.
double boundary[1000][2]; //The coordinates of the boundary of the object.
//*** The function declarations ***
//This function initializes the simulation. Modify the variables in this function to create your own
//simulations.
void InitializeSimulation();
//These functions calculate the centre of mass and the moment of inertia of the object assuming the
//boundary describes a uniformly distributed solid or hollow object (respectively). They translate the
//coordinates of the boundary so that the centre of mass lies at the origin.
void CalculatePropertiesSolid();
void CalculatePropertiesHollow();
//Calculates the lowest point of the object in its current orientation.
//Used to determine if the object has collided with the ground and the relative position of the point
//of contact.
double yMin();
//Calculates the new linear and angular velocity of the object after a collision with the ground.
void ResolveCollision();
//Updates the simulation and calls Display() every "TIMERMILLISECS" milliseconds.
void Timer(int value);
//Used by GLUT to handle window resizing.
//Sets the simulationHeight to 1 and the simulationWidth to width/height.
void Resize(int width, int height);
//Used by GLUT to handle keyboard input.
//Any keyboard input is interpreted as a command to reset the simulation.
void Keyboard(unsigned char key, int x, int y);
//Draws a line of thickness "lineThickness". Used by Display().
void DrawLine(double x1, double y1, double x2, double y2);
//Used by GLUT to draw graphics to the window.
void Display(void);
//*** The function definitions ***
//This function initializes the simulation. Modify the variables in this function to create your own
//simulations.
void InitializeSimulation()
{
int i;
//A few notes before we begin:
//The object does not have to be convex.
//We are assuming the mass of the object is 1 (this should make no difference to the simulation).
//The lower left corner of the window is (0,0).
//The upper right corner of the window is (simulationWidth, simulationHeight).
//simulationHeight is always 1. (This should give you a rough idea of the scale of the simulation.)
//Define the boundary of the object by setting noOfPoints and the boundary[][] array.
//The centre of mass is assumed to be at the origin. We can use CalculatePropertiesSolid() and
//CalculatePropertiesHollow() to calculate the centre of mass and translate the boundary
//appropriately at a later time if we so wish.
//Increase noOfPoints to get a circle, and change the 0.1 coefficients to get an ellipse.
noOfPoints = 4;
for(i=0;i0)
return; //The object is not colliding.
//Calculate the collision matrix K and its inverse invK.
//We are assuming the mass of the object is 1.
K11 = 1+ry*ry/I; K12 = -rx*ry/I;
K21 = -rx*ry/I; K22 = 1+rx*rx/I;
D = K11*K22-K12*K21; //The determinant of K.
invK11 = K22/D; invK12 = -K12/D;
invK21 = -K21/D; invK22 = K11/D;
//Initialize (ux,uy), Wc, and Wd.
ux = u0x;
uy = u0y;
Wc = 0;
Wd = 0;
//Step 3. Determining the type of ray we are on.
bSticking = 0;
bConverging = 0;
bDiverging = 0;
if(-tolerance0)
{
//Step 3b.
kx = -K11*cF+K12;
ky = -K21*cF+K22;
if(-tolerance0)
bDiverging = 1;
else
bConverging = 1;
}
}
else
{
//Step 3c.
kx = K11*cF+K12;
ky = K21*cF+K22;
if(-tolerance0)
bConverging = 1;
else
bDiverging = 1;
}
}
//Step 4. We are on a converging ray.
if(bConverging==1)
{
//Step 4a.
uOrigin = uy-ky*ux/kx;
if(uOrigin<=0)
{
//Step 4b.
Wc = -ux*uy/kx+ky*ux*ux/(2*kx*kx);
ux = 0;
uy = uOrigin;
bSticking = 1;
}
else
if(0=-tolerance)
{
//Just to be safe.
kx = 0;
ky = 1/invK22;
}
}
bDiverging = 1;
}
}
//Step 8. We are on a diverging/stationary ray (or stable sticking has occurred).
if(bDiverging==1)
{
//Step 8a.
uOld = uy;
if(uy<0)
{
//Step 8b. We are in a compression phase.
Wc += -uy*uy/(2*ky);
uy = 0;
}
//Step 8c. We are in a decompression phase.
uy = sqrt(2*ky*(-cR*cR*Wc-Wd)+uy*uy);
ux = kx*(uy-uOld)/ky+ux;
}
//Step 9. Calculate the impulse and the new velocities.
px = invK11*(ux-u0x)+invK12*(uy-u0y);
py = invK21*(ux-u0x)+invK22*(uy-u0y);
vx = v0x+px;
vy = v0y+py;
w = w0+(rx*py-ry*px)/I;
//Step 10. Rotate the space (is not necessary).
return;
}
//Updates the simulation and calls Display() every "TIMERMILLISECS" milliseconds.
void Timer(int value)
{
int newTime;
int elapsedTime;
double dt;
//Ask GLUT to re-call this function after "TIMERMILLISECS" millieconds.
glutTimerFunc(TIMERMILLISECS,Timer,0);
//Work out how much time has passed since the last time we updated the simulation.
newTime = glutGet(GLUT_ELAPSED_TIME);
elapsedTime = (newTime-prevTime);
prevTime = newTime;
//Update the simulation in intervals of 6 milliseconds.
while(elapsedTime>=6)
{
dt = 0.006;
elapsedTime -= 6;
//Update the position, orientation, and velocity of the object.
cx += vx*dt;
cy += vy*dt+g*dt*dt/2;
vy += g*dt;
a += w*dt;
//Check if a collision has occurred.
if(cy+yMin()