4. Simulation of Billiards
You may complete this portion of the assignment either individually or with a partner, following the
paradigm of pair programming. Assignments submitted by individuals and pairs will be graded identically,
so you should choose the programming arrangement you prefer.
Objective Your team's task for this week is to complete a MATLAB simulation of the game of billiards.
Your simulation will contain two billiard balls moving on a flat, perfectly frictionless table, able to collide
with each other and with the walls. The below left image shows the trajectories taken by the billiard balls in
a sample simulation. The below right image instantaneously depicts the linear momentum of the two balls
separately and together.
Billiards Simulation by KJK (Solutions)
Horizontal Position (m)
-0.05 0.05 0.1
Horizontal Momentum (kg'mls)
The system has already been set up for you in MATLAB so that you can focus on the dynamics; the
main script (simulate_billiards_starter. m) defines the system's parameters, runs the simulation, draws
the ball trajectories, animates their momentum vectors, and graphs their energy over time; all you need to
do is add the code that handles collisions. This script uses two subsidiary functions that you will not need
to modify ( compute_billiards_derivatives. m and check_billiards_collisions. m).
As you can see at the top of the main simulation file, the pool table consists of four walls located at
leftWallX, rightWallX, lowerWallY, and upperWallY. All distances in the simulation are measured in
meters, and velocities are measured in meters per second.
One of the billiard balls is white, and the other is black. Each ball has its own coefficient of restitution
e for impacts with the walls (whiteBallWallE and blackBallWallE). They have a separate coefficient of
restitution for collisions with one another (ballBallE). All collisions are frictionless.
Each ball is defined by its radius (whiteBallRadius and blackBallRadius) and its mass (whiteBallMass and
blackBallMass). The ball masses are currently programmed to be proportional to the volume of a
cylinder of the defined radius and fixed height, with both balls having the same density; thus the visual size
of a ball also indicates its mass.
You control the initial billiard ball positions and velocities by setting the appropriate variables at the
top of the code, e.g., whiteBallinitialPosition. Each position and velocity vector has two components
to reflect the planar nature of this simulation. In all such vectors, the first element is the x-component
(horizontal), and the second is the y-component (vertical).
Your job is to use your knowledge of particle impact dynamics to program the interactions between the
balls and the walls, and then to program the interaction between the two balls themselves. Your simulation
should work correctly for all possible combinations of reasonable parameter values. To help improve your
intuition for particle impacts, you should test a wide variety of scenarios as you develop your code.
Distinction from Previous Assignments It is important to note that the contacts you are programming
this week do not behave like springs, so they are very different from the wall you programmed in the puck
simulation assignment. Instead of operating on deflections, these particle impacts are instantaneous. As you
can read in Kasdin and Paley Chapter 4, instantaneous collisions create discontinuous velocities. We are
using the event detection capabilities of ode45 to identify the instants where collisions occur so that we can
handle them (calculate how they affect the ball states) before restarting the simulation. We could not create
such a simulation by writing only the state derivatives, as the derivative of the velocity vector would need
to have infinite magnitude for an infinitesimal duration during collisions, which is not compatible with our
regular simulation approach.
Dynamic Debugging In addition to setting up the problem, running the simulation, and graphing the
resulting ball trajectories, the main simulation code also graphs each ball's kinetic energy and the system's
total kinetic energy over time. You should use this energy graph to help judge whether your simulation is
behaving as it should. A sample energy graph is shown below; this is the energy graph for the simulation
shown in the figure on the previous page. You should also watch the momentum animation shown to the
right of the billiards table to check your simulation and deepen your understanding of collisions.
0.025 u:ro:o::,I I 111 1111 I 111 1111 I
I-Total Energy I --+-White Ball Energy
Black Ball Energy
f-ffi--ffif fff b
2 3 4 5
6 7 8 9
As mentioned above, the graphics and the simulation engine have already been programmed
for you so that you can focus on the dynamics. Your tasks are:
2. Set the value of the studentNames variable near the top of the code.
3. Run the code to make sure the simulation is set up correctly. The balls should both move
according to their initial conditions, but they will not react to collisions with the walls or each
other. If you watch the command line, you will see messages appear when collisions are
detected; your job is to write the code that modifies the ball states when a ball collision is
4. Look through the provided starter code (all three files) for at least fifteen minutes to see how
it is written. Select one aspect of the code that doesn't fully make sense to you, and figure it out.
Use MATLAB's documentation (doc), breakpoints, and the debugger to figure out how this
aspect of the code works; you can even make small modifications to the code to test out how it
5. Find the two commented lines of asterisks (**********) All of the code you write should be
between these two lines, which are lines 161 and 219 in the starter code. Do not modify any
code outside this asterisk zone.
6. Test out the sample collision responses included in the starter code by uncommenting the
indicated lines, such as xw = O; under whiteRight. The code listed here executes when a collision
occurs between the white ball and the right wall. These examples are meant to help you
understand how the simulation behaves.
7. Program the eight wall collisions first. Enforce the correct coefficient of restitution, and check
to make sure these impacts behave as you expect before you move on. Make sure your code is
8. Once you're done with the walls, program the response to the collision between the two balls;
this is the whiteBlack case in the switch statement. Enforce conservation of momentum and the
correct coefficient of restitution, and check to make sure the system behaves as you expect by
testing a variety of cases. Note that the MATLAB function dot may be useful for computing the
components of ball velocity parallel and perpendicular to the line of impact. Make sure your
code is well commented. It is likely that you will need to draw diagrams and write equations to
figure out the correct way to handle this collision. You do not need to turn in your written
analysis, but you should still do it.
9. If you cannot get the full ball-ball collisions to work correctly, program the simplified case
where the two balls have the same mass, and state that you did so in the submission email. This
simplification will earn you a maximum of 1.5 points out of the three possible on this task.
10. Use your simulation to investigate the following question, which you will answer in the
body of your submission email. Set up a situation where the two balls have different radii,
different masses, and a ball-on-ball coefficient of restitution equal to one. Given these
parameters, is it possible to have a moving ball A hit a stationary ball B so that A stops right
after the impact? If so, explain the initial conditions that yield this result. If not, briefly explain
why this result cannot be achieved. If you are not able to get your simulation working perfectly,
you may still answer this question based on your understanding of particle impact dynamics,
though you won't have the simulation available to test out your hypotheses.
This material may consist of step-by-step explanations on how to solve a problem or examples of proper writing, including the use of citations, references, bibliographies, and formatting. This material is made available for the sole purpose of studying and learning - misuse is strictly forbidden.
function [value, isterminal, direction] = check_billiards_collisions(~,state)
% This function determines whether two billiard balls have hit the walls or
% each other. There are nine possible collisions that can occur: each ball
% with each wall, plus the two balls together.
% The input t is time in seconds.
% The input state is our state vector. It has eight elements, arranged in a column.
% The state vector's first element, state(1), is the horizontal position of the white ball's center, in meters.
% The state vector's second element, state(2), is the vertical position of the white ball's center, in meters.
% The state vector's third element, state(3), is the horizontal velocity of the white ball's center, in m/s.
% The state vector's fourth element, state(4), is the vertical velocity of the white ball's center, in m/s.
% The state vector's fifth element, state(5), is the horizontal position of the black ball's center, in meters.
% The state vector's sixth element, state(6), is the vertical position of the black ball's center, in meters.
% The state vector's seventh element,state(7), is the horizontal velocity of the black ball's center, in m/s.
% The state vector's eigth element, state(8), is the vertical velocity of the black ball's center, in m/s.
% Declare wall locations as global.
global leftWallX rightWallX lowerWallY upperWallY
global whiteBallRadius blackBallRadius
global whiteLeft whiteRight whiteLower whiteUpper blackLeft blackRight blackLower blackUpper whiteBlack
% Pull our positions out of the state vector. They are all we need.
xw = state(1);
yw = state(2);
xb = state(5);
yb = state(6);
% The value of each test is the distance from the ball in question to the
% object with which we are checking, making sure to subtract out the ball
% radius. These are all positive values that become zero when the
% collision occurs.
value(whiteLeft,1) = xw - leftWallX - whiteBallRadius;
value(whiteRight,1) = rightWallX - xw - whiteBallRadius;
value(whiteLower,1) = yw - lowerWallY - whiteBallRadius;
value(whiteUpper,1) = upperWallY - yw - whiteBallRadius;
value(blackLeft,1) = xb - leftWallX - blackBallRadius;
value(blackRight,1) = rightWallX - xb - blackBallRadius;...