Jarlold's Physics Cook Book¶
These notes are designed to hold the simplest (easiest to understand?) implementations I could think of for physics simulators in "Game Physics Engine Development". They are not meant to be computationally efficient, and are thus written in Python. It is my intent to implement these examples using no special libraries at all. That rule does not apply to rendering the simulations, for which I will primarly use MatPlotLib.
The idea behinds this "Cook Book" is that you can easily read and convert the contents into whatever language you are using, in order to implement a functional physics simulator with ease of effort. Thus the explanations will be kept brief, and notes will only really be used for when I deviate from the book, or as a quick reminder of the contents of the book. If you are seeking to actually understand how a physics engine works, then maybe consult the original explanations in the book (it's a real page turner!).
First Some Math¶
Basic vector arithmetic. You probably won't have to implement this yourself unless you're doing something really funny. Like that guy who ran mario on a dishwasher.
class Vector3:
def __init__(self, x=0, y=0, z=0):
self.x, self.y, self.z = x, y, z
# Get how long the vector is
def magnitude(self):
return (self.x**2 + self.y**2 + self.z**2)**0.5
# Make the length = 1 without changing the direction
def normalize(self):
if (self.x, self.y, self.z) == (0, 0, 0):
return Vector3()
if (length := self.magnitude()) > 0:
return self.scalar_multiply(1/length)
# Rotate the vector 180 degrees (point backwards)
def invert(self):
a = Vector3(self.x, self.y, self.z)
a.x *= -1
a.y *= -1
a.z *= -1
return a
# Makes the vector `scalar` times longer, in the same direction
def scalar_multiply(self, scalar):
a = Vector3(self.x, self.y, self.z)
a.x *= scalar
a.y *= scalar
a.z *= scalar
return a
def euclidean_distance(other):
return (
(self.x - other.x)**2 +
(self.y - other.y)**2 +
(self.z - other.z)**2)**0.5
# Add the components of each vector
def __add__(self, vect):
a = Vector3(self.x, self.y, self.z)
a.x += vect.x
a.y += vect.y
a.z += vect.z
return a
# Like above but negative
def __sub__(self, vect):
a = Vector3(self.x, self.y, self.z)
return a + (vect.invert())
# Returns the scalar product (sometimes also called the dot product
def __mul__(self, vect):
if type(vect) == Vector3:
return self.x*vect.x + self.y*vect.y + self.z*vect.z
else:
return self.scalar_multiply(vect)
# Returns the vector product (also called the cross product)
def __mod__(self, vect):
return Vector3(
self.y * vect.z - self.z*vect.y,
self.z * vect.x - self.x*vect.z,
self.x * vect.y - self.y*vect.x)
# Makes debugging a bit easier
def __str__(self):
return "({}, {}, {})".format(self.x, self.y, self.z)
def __eq__(self, vect):
return self.x == vect.x and self.y == vect.y and self.z == vect.z
def __repr__(self):
return "Vector3"+self.__str__()
And here are some test cases, for if you do decide to implement it on your own.
a = Vector3(1, 2, 3)
b = Vector3(3, 2, 1)
print("Subtraction:", (a-b) == Vector3(-2, 0, 2))
b *= 2
print("Assigned multiplication: ", b == Vector3(6, 4 ,2))
print("Multiplication assignment: " + str(b == Vector3(6, 4, 2)))
print("Scalar Product", a * b == 20)
print("Vector Product: ", (Vector3(1, 2, 3) % Vector3(1, 5, 7)) == Vector3(-1, -4, 3))
Subtraction: True Assigned multiplication: True Multiplication assignment: True Scalar Product True Vector Product: True
Point Mass Physics¶
I'm using "Point Mass" because I think the title "Particle Physics" would be misleading. But I'll switch back to using the word particle in a moment.
Basic Equations¶
We'll implement a particle physics system that uses Newtonian physics and a force accumulator. Essentially we're just making a Vector3, and applying the following mathematical rules to it: $$\vec{F}_{net} = \Sigma_{i=0...n}\vec{F}_n$$ $$\vec{a} = \frac{\vec{F}_{net}}{m}$$ $$\vec{v} = \int_0^t{\vec{a}}$$ $$\vec{p} = \int_0^t{\vec{v}}$$
Where:
- $\vec{F}$ is the force generated by some constraint (like a spring, or gravity).
- $\vec{F}_{net}$ is the sum of the above foces.
- $\vec{a}$ is of course, acceleration (in the book denoted $\ddot{p}$.
- $\vec{v}$ is velocity, (in the book denoted $\dot{p}$.<br.
- $\vec{p}$ is the position of the particle.
Integration Hacks¶
Here, when we write $\int$, we're actually referring to a discrete sum of boxes on the function. It's essentially the first method we're taught for estimating the integral in calculus 1.
So every tick, we'll multiply the y value (our acceleration or velocity) by $\delta t$, which is the "framerate" of our simulation (here it's set to be 60).
Floating Point Problems¶
Note that in our code, since we always divide by mass and never use it in it's original form, we'll store it as the inverse (so we can just multiply it). This helps with floating point numbers being generally shit.
Python Force Registry¶
We'll also register our forces in the particle object itself, instead of having an external force generator registry. If I were to do this in c++, I'd use a hashmap and an external class or a static variable for this (for performance reasons mentioned in the book). However, since I've written this in Python, there's actually very little performance bonus from switching to such a method.
class Particle:
def __init__(self, v=Vector3(), p=Vector3(), mass=1):
if type(p) != Vector3:
raise TypeError("Passed a non-vector as position.")
if type(v) != Vector3:
raise TypeError("Passed a non-vector as velocity.")
self.velocity = v
self.position = p
self.inverse_mass = 1.0/mass
self.force_registry = [] # A list of callables that generate forces.
# In languages without call overloading, use a class.
def register_force(self, force_generator):
self.force_registry.append(force_generator)
def tick(self, delta_t):
# Compute F_net as the sum of all the forces
f_net = Vector3()
for force_generator in self.force_registry:
f_net += force_generator(self)
# We're doing our integration by calling realy small time frames, and making a box
# estimate. Similar to how we're first taught in calculus 1.
# Then update the acceleration according to the force
acceleration = f_net * self.inverse_mass
# Update the velocity according to the acceleration
self.velocity += acceleration * delta_t
# Update the position in accordance to the velocity
self.position += self.velocity * delta_t
Now let's do a little test by simulating one particle in 3D space with a few forces applied to it. We'll graph it and make sure it looks like it's working.
import matplotlib.pyplot as plt
import numpy as np
import math
%matplotlib ipympl
ax = plt.figure().add_subplot(projection='3d')
# Create a particle with an initial velocity of 10 in every direciton
p = Particle(Vector3(10, 10, 10))
# Register gravity as a force that always returns (0, -9.81, 0)
p.register_force( lambda x: Vector3(0, -9.81, 0) )
# For shits and giggles, let's add some wobbly forces for no reason
#p.register_force( lambda particle: Vector3(0, -100.81*math.sin(particle.position.x), 0) )
#p.register_force( lambda particle: Vector3(-2, 0, 2) )
# Run the simulation and graph it in matplotlib
x, y, z = [], [], []
for i in range(60*5): # We'll run the simulation at 60tps
p.tick(1.0/60)
x.append(p.position.x)
y.append(p.position.y)
z.append(p.position.z)
ax.scatter(x, z, y) # Matplotlib has a different opinion about the Z axis
<mpl_toolkits.mplot3d.art3d.Path3DCollection at 0x755dfa01dd90>
Some Force Generators¶
Recipes for some force generators that are generally useful (and were detailed in the book :P)
Gravity¶
This is really better as a simple function, or even better as a built-in to the physics particle class (for performance reasons). I've done it here as a callable class as the most basic example of a force generator.
class GravityGenerator:
def __init__(self, strength):
self.force = Vector3(0, -strength, 0)
def __call__(self, particle):
return self.force*(1.0/particle.inverse_mass)
g = GravityGenerator(9.81)
str(g(p))
'(0.0, -9.81, 0.0)'
Drag¶
We often estimate the force of air resitance as a quadratic function of the speed of the object. We can implement that like below. This is another force generator that would be more efficient built into the particle class, but for consitency case, I've done it as a callable.
class DragGenerator:
def __init__(self, linear_coefficient, quadratic_coefficient):
self.l = linear_coefficient
self.q = quadratic_coefficient
def __call__(self, particle):
# The magnitude of the force of air resitance is dependent
# on the speed and the square of the speed
speed = particle.velocity.magnitude()
drag_coefficient = self.l * speed + self.q * (speed**2)
# Then we'll apply that to all axis
return particle.velocity.normalize() * -drag_coefficient
d = DragGenerator(0.2, 0.001)
p = Particle(v=Vector3(10, 30, 100))
str(d(p).magnitude())
'31.976176963403034'
Springs¶
Now we'll add in springs. This is the first case of a force generator that actually benefits from being a force generator. The force of a spring increases linearly as we move away from the springs resting point (Hook's law). Note that the way we implement this, each particle bound by the spring needs it's own spring generator!
class SpringGenerator:
def __init__(self, other_particle, spring_constant, resting_length):
self.other_particle = other_particle
self.spring_constant, self.resting_length = spring_constant, resting_length
def __call__(self, particle):
# Compute the distance between the two particles
delta_p = particle.position - self.other_particle.position
length = delta_p.magnitude()
# Compute the force of the spring.
# Difference between current distance and resting length, times the spring constant
force_magnitude = abs(length - self.resting_length) * self.spring_constant
r = delta_p.normalize()
# Now we'll calculate the direction of the spring force
return r * -force_magnitude
p1 = Particle(p=Vector3(2, 2, 2))
p2 = Particle(p=Vector3(2, 2, 2))
sg = SpringGenerator(p1, 0.1, 1)
a = sg(p2)
print(a)
(-0.0, -0.0, -0.0)
Anchored Springs¶
Sometimes we would like to pretend the spring is stuck to a point of infinite mass. Rather than compute the infinite mass back and forth, we can just assume it doesn't move and create something like the following.
class AnchoredSpringGenerator:
def __init__(self, anchor_point, spring_constant, resting_length):
self.anchor_point = anchor_point
self.spring_constant = spring_constant
self.resting_length = resting_length
def __call__(self, particle):
# Distance between particle and anchor point
distance = particle.position - self.anchor_point
displacement = distance.magnitude() # displacement is vectorless distance
# Magnitude of the pulling force generated is porportional to how far away
# the particle is from the resting length
magnitude = abs(displacement - self.resting_length) * self.spring_constant
# Return that magnitude in the opposite direction as the distance
return distance.normalize() * -magnitude
# Make an anchor that tries to pull to 0, 2, 0, with a spring constant 1
# and a resting length of 1
a = AnchoredSpringGenerator(Vector3(0, 2, 0), 1, 1)
# Then make a particle at 0, 0, 0 (2 units away) and see if the
# spring pulls it by (2 - 1)*1 = 1 units
a(Particle())
Vector3(-0.0, 1.0, -0.0)
Putting it All Together¶
Let's try simulating two particles suck together by a very strong spring. One particle will be launched rather firmly in an arc. Both objects will experience normal gravity.
p1 = Particle(p=Vector3(2, 2, 2), v=Vector3(0, 10, 30))
p2 = Particle(p=Vector3(2, 2, 2), mass=20)
p1.register_force(SpringGenerator(p2, 100, 2))
p2.register_force(SpringGenerator(p1, 100, 2))
p1.register_force(GravityGenerator(9.81))
p2.register_force(GravityGenerator(9.81))
p1_hist = []
p2_hist = []
for _ in range(60*1):
p1.tick(1.0/60)
p1_hist.append( (p1.position.x, p1.position.y, p1.position.z) )
p2.tick(1.0/60)
p2_hist.append( (p2.position.x, p2.position.y, p2.position.z) )
import matplotlib.pyplot as plt
import numpy as np
import math
%matplotlib ipympl
ax = plt.figure().add_subplot(projection='3d')
p1_xs, p1_ys, p1_zs = zip(*p1_hist)
p2_xs, p2_ys, p2_zs = zip(*p2_hist)
ax.scatter(p1_xs, p1_zs, p1_ys) # P1 in blue
ax.scatter(p2_xs, p2_zs, p2_ys) # P2 in orange
<mpl_toolkits.mplot3d.art3d.Path3DCollection at 0x755df9b1cc20>
It's a bit hard to get a "feel" for whether things are working or not. But we can go through a checklist of things that we know should be true:
- The particle with greater mass is less effected by forces.
- The lighted object is more affected by forces.
- Both particles are affected by gravity, and gravity accelerates them equally.
- The particles seem bound together by the spring, and curve to try and be closer. The lighter particle does the majority of the curving
Maybe for a more useful example, let's try simulating the terminal velocity of a point mass falling from a table after being pushed by a cat.
p1 = Particle(p=Vector3(0, 100, 0), v=Vector3(0.1, 0, 0))
p1.register_force(GravityGenerator(9.81))
p1.register_force(DragGenerator(0.1, 0.1))
p1_hist = []
p1_velocity_hist = []
for _ in range(60*5):
p1.tick(1.0/60)
p1_hist.append( (p1.position.x, p1.position.y, p1.position.z) )
p1_velocity_hist.append(p1.velocity.y)
x, y, z = zip(*p1_hist)
ax = plt.figure().add_subplot(projection='3d')
ax.scatter(x, z, y)
<mpl_toolkits.mplot3d.art3d.Path3DCollection at 0x755df9b1e090>
ax = plt.figure().add_subplot()
#ts = [i/100.0 for i in range(0, 100*60*5, 100/60.0)]
ts = [ i * 1/60.0 for i in range(len(p1_velocity_hist)) ]
plt.xlabel("Time (minutes)")
plt.ylabel("Velocity in Y Direction")
plt.scatter(ts, p1_velocity_hist)
plt.axhline(y=p1_velocity_hist[-1], color='r', linestyle='-')
plt.text(-0.2, p1_velocity_hist[-1] + 0.25, "Terminal Vel. :" + str(p1_velocity_hist[-1]))
plt.show()
So it looks like the terminal velocity for our point mass is approximately -9.42 unkown units that we didn't specify. Isn't that useful to know?
Finally let's try a particle stuck to an anchored spring. This should oscilate, like those things they put behind doors. Those springy things, you know what I'm talking about right?
p1 = Particle(v=Vector3(20, 0, 0)) # Give that boy some momentum!
p1.register_force(AnchoredSpringGenerator(Vector3(-2,0,0), 1, 2))
#p1.register_force(GravityGenerator(9.81))
p1_hist = []
for _ in range(60*1):
p1.tick(1.0/60*20)
p1_hist.append( (p1.position.x, p1.position.y, p1.position.z) )
xs, ys, zs = zip(*p1_hist)
ax = plt.figure().add_subplot()
plt.xlabel("Time Frame")
plt.ylabel("X Pos. of Particle")
plt.plot(xs)
[<matplotlib.lines.Line2D at 0x755df9a747d0>]
ax = plt.figure().add_subplot(projection='3d')
ax.scatter(xs, zs, ys)
<mpl_toolkits.mplot3d.art3d.Path3DCollection at 0x755df9adcbf0>
And as I remember from the physics class I had to drop, that is accurate. If we set it to a non-harmonic frequency, it even has the small and large peaks! (Guess I didn't need to enter the complex plane to simulate it after all.)
Hard Constraints¶
Our current force-generator simulation has a lot of limitations. In real life collisions act like springs, even very hard objects bend a little and result in some springy-ness. This unfrotunately, happens very quickly, and with very large, or very small, numbers. That makes it unrealistic to do it in our physics simulator. Instead, we'll develop a new system comprised of the following main components:
- Collision resolution: Which determines the velocities of objects after they bump into each other.
- Interpenetration resolution: Which moves objects outside of each other in the event that they wind up going inside each other.
- Resting contact detection: Which prevents resting objects from constantly generating collision checks.
Collision Detection, Closing Velocity, and Contact Normals¶
If the distance between two sphere is less than the sum of their radii, it should reason that those two spheres are colliding.
But simply saying "there is a collision" isn't enough, we also need to know how quickly, and how they should separate.
How quickly they collide is called the "closing velocity" and it is computed as such
$$
v_s = \dot{p}_a \cdot \widehat{(p_b - p_a)} + \dot{p}_b \cdot \widehat{(p_a - p_b)}
$$
where $v_s$ represents the rate at which the objects are moving away from each other.
How much of that velocity should we assign to each object? That's called the restitution factor, and we'll simply declare it as a coefficient on the separating velocity $v_s$ $$ v\prime_s = -cv_s $$
Now we have the task of computing in what direction we should be moving the objects to un-collide them. This is actually very complicated for non-trivial shapes, but for spheres the following equation suffices. $$ \hat{n} = \widehat{(p_a - p_b)} $$
Putting it all together, we get the final equation for the exit velocity as a vector, which is $$ v_s = (\dot{p}_a - \dot{p}+b) \cdot \hat{n} $$ Conveniently, this is a much simpler and faster to computer equation.
class RigidSphere:
def __init__(position: Vector3, velocity: Vector3, scale=1, mass=1):
self.position = position
self.velocity = velocity
self.inverse_mass = 1.0/mass
self.radius = 1
self.force_registry = []
def register_force(self, force_function):
self.force_register.append(force_function)
def get_acceleration(self):
# net force is sum of all other forces
f_net = sum([ fgen(self) for fgen in self.force_registry])
# f = ma
return f_net * self.inverse_mass
def is_colliding(self, other):
return self.position.euclidean_distance(other) < (self.radius + other.radius)
When two sphere collide, their velocities change. This function does that bit.
def resolve_sphere_velocity(sphere1, sphere2, dt, restitution=0.5): # Untested
# This part only works for spheres, the rest should be ok though
contact_normal = (sphere1.position - sphere2.position).normalize()
# Figure out how quickly the objects are moving apart
separating_velocity = contact_normal * (sphere1.velocity - sphere2.velocity);
# Make sure they're actually colliding... and not drifting apart on their own
if separating_velocity > 0:
return
# Once they bounce off they'll leave with some amount of their original energy
# (but not all of it)
new_sep_velocity = -separating_velocity * restitution
# This is part of resting object detection. Instead of just checking if the objects
# are moving towards eachother, we'll check if they're accelerating to eachother as well.
sep_velocity_from_acceleration = sphere1.get_acceleration() - sphere2.get_acceleration()
sep_velocity_from_acceleration *= contact_normal * dt;
# So based off the acceleration, we're going to try and look one frame ahead,
# and if the object would collide again (because of the acceleration), we'll
# factor that into our adjustment. This hack allows us to avoid implementing normal
# forces.
if sep_velocity_from_acceleration < 0:
new_sep_velocity += restitution * sep_velocity_from_acceleration
if new_sep_velocity < 0: new_sep_velocity = 0
total_imass = sphere1.inverse_mass + sphere2.inverse_mass
# If both particles have negative mass (our flag for infintie mass)
# then just like my real life, nothing we do here matters.
if (total_imass <= 0): return
# Ratio of how much to speed the objects up by based off their mass
# I'm not a physicist but like, mass*velocity is momentum right?
# and we're like, conserving that momentum?
# so it's like... conservation of like.. momentum?
impulse_per_imass = contact_normal * (new_sep_velocity - separating_velocity)/total_imass
sphere1.velocity += impulse_per_imass * sphere1.inverse_mass
sphere2.velocity += impulse_per_imass * sphere2.inverse_mass
When two sphere collide, more often than not, one will wind up inside the other before we get around to resolving the collision. So this moves one outside of the other.
def resolve_sphere_interpenetration(sphere1, sphere2, dt): # Untested
total_imass = sphere1.inverse_mass + sphere2.inverse_mass
# Your mama so fat that her inverse mass is a negative
# number (our flag for infinite inverse mass)
if total_imass <= 0: return
# Difference to their middles
penetration = sphere1.position.euclidean_distance(sphere2.position)
penetration -= sphere1.radius + sphere2.radius # This right, right?
# From each to according to his displacement, to each according to his inverse mass
movement_per_imass = contactNormal * (-penetration / total_imass)
sphere1.position += movement_per_imass * sphere1.inverse_mass
sphere2.position += movement_per_imass * sphere2.inverse_mass
def resolve_sphere_collision(sphere1, sphere2, dt):
resolve_sphere_velocity(sphere1, sphere2, dt)
resolve_sphere_interpenetration(sphere1, sphere2, dt)
Contact Resolver Algorithm¶
In order to resolve all the contacts in our simulation, we have to do the following:
- Adjust the velocities
- Fix any interpenetrations for all the collisions/contacts. This leaves us with the question of what order and how many times to do each step. If we resolve the velocity or interpenetration of one collision, it might affect the separating velocity or interpenetration of a different collision. Fortunately some nerd did a bit of math, and proved that in frictionless scenarios, if we just keep fixing the interpenetrations, eventually will fix them all.
A Better Visualization Method¶
We're gonna need to be able to generate animations of points, and later maybe animations of primitive shapes, so I'll briefly put together a library for doing that. I've written it in OpenGL, os it won't render properly in Jupyter. Maybe there's some way to force it to record a video and embed that, but for now just run it yourself!
Unfortunately this is an OpenGL implementation, and does not integrate nicely with Jupyter notebooks, so you'll have to run things for yourself to test. I'll expand on this side-library as needed, but will generally try and keep it simple.
This just renders a rotating cube and a translating sphere as a demonstration.