Angry Birds Space: A Custom Physics Engine
Demonstration of the game for several, different initial conditions
Motivation
After I took my first physics class in high school, I was deeply unsatisfied. Studying the motion of free-falling cannonballs was interesting, but the simplistic models we were constructing were far from the more complex ones that make up our reality. Determining the motion of a water droplet falling from leaf to leaf or a cardboard box thrown at a wall, at the time, seemed like an impossible task. Fast-forward to my first year at Northwestern, and I was studying topics like force balance in a mass-spring-damper system, which was simple at first but became incredibly tedious as the more masses, springs, dampers were added. I was astounded. Surely, there must be another approach to analyzing complex systems that doesn’t require spending dozens of hours drawing diagrams and manipulating equations.
There was. Enter the fall quarter of my second year, and I enrolled in the MECH_ENG 314 course, titled “Machine Dynamics.” Here, I learned about an elegant approach to modeling physical systems: Lagrangian mechanics. Albeit an incomplete theory still being researched today, Lagrangian mechanics lends itself to understanding dynamic systems, including ones that are more conveniently represented in non-Cartesian coordinates (like a double pendulum or quad-jointed robot), with multiple constraints and rigid bodies impacting one another.
Being one of my favorite class I’ve taken thus far, I decided, for my final project, to tackle something ambitious and fun: a replication of an Angry Birds Space level. Angry Birds Space is a physics-based game developed by Rovio Entertainment that is similar to Angry Birds but different because instead of having one constant gravitational force that acts in one direction, it has multiple pulling objects towards the center of planets.
The Inspiration: Angry Birds Space Game by Rovio Entertainment
Design Specifications
The level I created has one bird, two pigs, and two planets. The bird and pigs were modeled as squares with rotational inertia. The bird and the pigs can collide with one another and also with the planets’ surfaces. The planets, assumed to have infinite mass relative to the bird and pigs, are fixed and will not be modeled as dynamic bodies. Each of the planets will have a gravitational field that acts towards the planet’s center. The planets’ gravitational fields only act on bodies within a distance from the planet’s center. Limiting the effect of each planet’s gravitational field leads to more interesting game dynamics. For instance, if the bird is not aimed properly and doesn’t enter the gravitational field, the direction of the bird’s velocity remains unchanged. Also, the bird, if its velocity is large enough, can enter, curve, and escape a planet’s gravitational field.
Unlike the actual Angry Birds Space game, the pigs do not disappear after being hit by the bird. I want to model the impacts between the bodies. If the pigs can be “defeated,” then there are only impacts between the bodies and the planets. In addition, I modeled the loss of energy due to drag forces that only act on bodies inside a planet’s gravitational field and plastic impacts with a coefficient of restitution, alpha. I made the assumption that space is a perfect vacuum and no drag force acts on objects in space. In other words, drag force only acts on bodies within one of the planets’ gravitational fields.
There are nine degrees of freedom in my system. They are the x coordinates, y coordinates, and angles from vertical of the three bodies. I defined four frames: the world frame and three translated and then rotated frames, one for each body. I labelled these frames numerically: {0}, {1}, {2}, and {3}. The transformation matrices between the frames are labelled as g01, g02, and g03 (a transformation from the world frame, {0}, to body frame {N}.
Rigid Body Frames and Configuration Variables
Calculations
Because I want each planet’s gravitational field to be localized, each body’s potential energy changes depending on its location. Typically, the Lagrangian is calculated by subtracting the total kinetic energy by the total potential energy. However, this approach would necessitate that I write one Euler-Lagrange (E-L) equation for every combination of the bodies’ locations. There is a faster approach to determining the equations of motion: writing down the Lagrangian for each body independent of the other bodies. Then, I can calculate the equations of motion for (1) the body in space, (2) the body in planet 1’s gravitational field, and (3) the body in planet 2’s gravitational field. One body at a time, I can determine the position of the body and use the corresponding equations of motion to integrate that body’s acceleration. If I did this using the traditional method, then I would need to calculate 27 E-L equations and evaluate, potentially, 27 conditionals in the simulate function. Since I chose the same dimensions and density for each body, I only need to calculate a total of 3 E-L equations using my non-traditional approach.
To calculate the E-L equations, I first found the body velocities by unhatting the product of the inverse of the transformation matrices and the time derivative of the transformation matrices. Then, I calculated the inertia matrices by taking a triple integral over the volume of a cube. Using the body velocities and inertia matrices, I calculated the kinetic energy (see image to the left). The gravitational potential energy for an object in a planet’s gravitational field is equal to V = -km⁄r where k is a parameter relating to the strength of the planet’s gravitational pull, m is the mass of the body, and r is the distance between the center of mass of the body and the center of the planet. The Lagrangian is equal to the kinetic energy minus the potential energy.
General Mechanics Equations
Within each planet’s gravitational field, the bodies experience a drag force. Friction can be defined as the negative velocity gradient of the general dissipation function (see image on the left). Drag force is proportional to the square of a body’s velocity. Hence, the dissipation function for drag is D = 1⁄6·ρ·Cd·A·v3 where ρ is the density of the air, Cd is the drag coefficient (~0.8 for cubes), A is the reference area (s2 for cubes), and v is the velocity of the body. See the image on the right for how the dissipation function is factored into the forced E-L equations.
To reduce the amount of work I had to do, I made the assumption that impacts between two bodies can only occur when they both occupy the same planet’s gravitational field or space. In nearly all cases, this assumption holds. With this assumption, I only need to write down three sets of impact equations: (1) impacts in space, (2) impacts in planet 1’s gravitational field, and (3) impacts in planet 2’s gravitational field. I can choose which Lagrangian (with the correct potential energy for the impacting bodies) to use, since I know where the impact is occurring.
There are three bodies (all of which are squares) and two planets. Therefore, there are 72 impact conditions and ϕ(q) equations. I wrote an algorithm to calculate these constraints for me to minimize the risk of human error. For the impact conditions between a corner of body 1 and a side of body 2, ϕ(q) equals (p1-p2)2 + (p2-p3)2 - (p1-p3)2 where p1 and p3 are the coordinates of body 2’s corners that, when connected, form a side and p2 is the coordinate of body 1’s corner. For the impact conditions between the corners of a body and a planet’s surface, ϕ(q) = (x-xcen)2 + (y-ycen)2 - R2 where x and y are coordinates of the body’s corner, xcen and ycen are coordinates of the planet’s center, and R is the radius of the planet. The impact equations (see image), along with the corresponding impact condition ϕ(q), give the resulting velocities of any two bodies colliding with each other. Due to the complexity of the impact equations, the solutions for q̇(τ+) were solved for numerically at the time of the impact. To model plastic impacts between the bodies and the planet, I used the following heuristic: I solved for the elastic case, replaced lambda with α·λ, and then solved again for q̇(τ+).
General Dissipation Function, Modified E-L Equations, and Impacts
Implementation
The SymPy package for Python was used to perform most of the symbolic math. Numerically solving the impact equations is computationally intensive. The purpose of this project was to apply my knowledge of Lagrangian mechanics, not to write fast, optimized code. That being said, it may take a few minutes (approximately 3) to run one simulation. Below is a snippet of code that constructs the system of Euler-Lagrange equations. L is the system’s Lagrangian. funcs are a vector of length 3N, where the first N elements are the configuration variables, which the system has a total of 9, the next N elements are their derivatives (velocity), and the last N elements are the second derivatives (acceleration). D is the dissipation function and t is the symbolic variable for time. The code is using SymPy’s built-in functions to take derivatives with respect to different variables. Then, the terms are aggregated back into a symbolic equation, forming the forced E-L equations.
def euler_equations(L, funcs, D, t):
N = len(funcs)//3
q = Matrix([funcs[0:N]])
qdot = Matrix([funcs[N:2*N]])
J = Matrix([L])
dJdq = J.jacobian(q)
dJdqdot = J.jacobian(qdot)
d_dJdqdot_dt = dJdqdot.diff(t)
D_mat = Matrix([D])
dDdqdot = D_mat.jacobian(qdot).T
lhs = (d_dJdqdot_dt-dJdq).T
EL_eqn = sym.Eq(lhs,dDdqdot)
return EL_eqn
A significant portion of the code was written using Python’s eval() function, which dynamically executes strings inputs. This was done due to the sheer number of equations of motion and impact conditions that I needed to keep track of. Instead of defining over hundreds of variables, I devised a variable naming scheme that would allow me to generate these equations without having to verbosely write every one of them out. Take the code below that formulates the 72 impact conditions between the three squares. We are evaluating ϕ(q), the version for impacts between rigid bodies, for every combination of body 1’s corner and body 2’s side. The naming scheme is written as “phi_”, followed by the ID of body 1 (an integer 1 to 3), the ID of body 2 (also, an integer 1 to 3), the ID of body 1’s corner (an integer 1 to 4), and lastly the ID of body 2’s side (also, an integer 1 to 4).
exec("body_pairs=combinations(list(range(1,{N1})),2)".format(N1=i+1),globals(),_locals)
for pair in _locals['body_pairs']:
for corner in range(1,5):
for side in [(1,2),(2,3),(3,4),(4,1)]:
exec("phi_{b1}{b2}{c}{p1}=(((r{b2}_c{p1}-r{b1}_c{c}).T*(r{b2}_c{p1}-r{b1}_c{c}))[0])**0.5+(((r{b1}_c{c}-r{b2}_c{p2}).T*(r{b1}_c{c}-r{b2}_c{p2}))[0])**0.5-(((r{b2}_c{p1}-r{b2}_c{p2}).T*(r{b2}_c{p1}-r{b2}_c{p2}))[0])**0.5\nphi_eqns.append(phi_{b1}{b2}{c}{p1})".format(b1=pair[0],b2=pair[1],c=corner,p1=side[0],p2=side[1]),globals(),_locals)
Simulation Outcomes
I ran multiple simulations, each with different initial conditions, to check that the impacts between all the bodies were working properly. I will only describe and rationalize the simulations titled “1 Pig Impact,” “Test W/ Alpha=1.0,” and “Test W/ Alpha=0.85” in the video. The other simulations appear reasonable for the same reasons I will give for the simulations I choose to describe below.
In the “1 Pig Impact” simulation, the pigs are initially at rest in space. Both are outside both planet’s gravitational field. At first, they remain stationary, since no gravitational forces or drag forces exist outside the planets’ gravitational field. The bird is also initially in space but is traveling with a positive horizontal velocity. As expected, the bird continues to travel in the same direction, at the same speed, since, in space, there are no gravitational forces or drag forces acting on the bird. When the bird enters the gravitational field of planet 1, it begins to curve as excepted due to the gravitational force acting inwards towards the center of the planet. I also observe that the bird slows down slightly due to the drag force now acting on the bird. The initial horizontal velocity of the bird is great enough that the bird is able to escape planet 1’s gravitational field after curving a bit. Once it enters space again, the bird travels at the same speed and direction it had when it left planet 1’s gravitational field. Shortly after, the bird collides with one of the pigs, which sends the pig flying towards planet 2. When the pig enters planet 2’s gravitational field, it begins accelerating towards the surface and impacts the surface with considerable speed. The pig bounces off the surface, escapes planet 2’s gravitational field, and flies off into space. This makes sense because most of the kinetic energy from the pig’s impact with the bird was transferred to the pig, providing it enough energy to impact with planet 2’s surface and escape the planet’s gravitational pull. Meanwhile, the bird, after colliding with the pig, is redirected towards planet 1. Having lost most of its kinetic energy to the pig, the bird bounces along the surface of planet 1, pulled towards the center of planet 1 by gravity and unable to escape planet 1’s gravitational field. The second pig does not experience impacts with the other pig or the bird, so it remains stationary.
To check that the plastic impacts were working, I compared two simulations (“Test W/ Alpha=1.0” and “Test W/ Alpha=0.85”) with the same initial conditions but different values of α, the coefficient of restitution. In the simulation with α = 0.85, the bird loses energy upon contact with the planet’s surface and bounces to a lower height compared to the simulation with α = 1.0. Overall, all the simulations I ran behaved in the way I expected and support that my code works.
