Spring Simulation
Contents
About this Project
Simulating springs is generally pretty straight-forward in terms of its mathematics and programming. The displacement can be used to calculate the force, and the velocity can be used to calculate the next position as well as the damping force. This is very useful for simulating springs in a complex setting or where very small time steps (in comparison to its natural frequency) are possible - not so much, however, in scenarios where no guarantees about the timesteps can be made.
The motivation for this project, for instance, was simulating gun recoil for a game. When faced with this challenge, I remembered how the recoil in Half-Life 2 felt very intuitive and smooth to me, and I decided to mimic this behavior. The best way I knew how was to put both the camera and the gun on a set of 'springs' that are given some velocity every time the gun fires. Perhaps the most obvious example of this is the .357 magnum:
Important to notice in that clip is that the camera does not 'jump' to a new position as many indie games, and especially shooter games on Roblox, would have done. Instead, the camera tilts upwards and pans slightly to the side, and then returns to close to its original position. Likewise, the gun has similar movement relative to the camera. Ultimately, we'd implement 5 springs; two for the camera (tilt and pan) and three for the gun (tilt, pan and backwards movement). The final product looks like this:
Now that we've seen what the goal is and what the final result looks like, we can explore how we got there.
Basic Simulation
You may be familiar with the basic implementation of a spring. For the sake of completeness, I'll put an example here:
-- @Author: AtAltitude <[email protected]> --
-- @Description: Basic mass-on-a-spring physics simulation --
-------------------------------------------------------------
--Create new table with __index meta-field to allow access to methods from objects
local Spring = {}
Spring.__index = Spring
--Constructor that takes spring constant (K) and damping (D) as inputs
function Spring.new(K, D)
--Initialize a new Spring object, with the fields:
local new = {
K = K; --Spring constant
D = D; --Damping constant
Position = 0; --Position of the spring; 0 is no compression/tension
Velocity = 0; --Velocity of the end of the spring; 0 means no movement
}
--Set the metatable and return the new object
return setmetatable(new, Spring)
end
--Method to update the spring simulation
function Spring:step(deltaT)
--Add position
self.Position = self.Position + self.Velocity * deltaT
--Apply force - for this, we define that the mass at the end of the spring is 1 kg
--This calculates the force resulting from both spring compression/tension and
--damping on the same line
local force = -(self.Position * self.K + self.Velocity * self.D)
self.Velocity = self.Velocity + force * deltaT
--Return position and velocity
return self.Position, self.Velocity
end
--Method to add velocity to the spring
function Spring:addVelocity(v)
self.Velocity = self.Velocity + v
end
--Return contents of the module
return Spring
This works well if you're able to call Spring:step() with very small values of deltaT compared to its natural frequency. However, this has some obvious downsides: For one, you have to keep calling Spring:step() to run the simulation, and the simulation quality quickly deteriorates with larger time steps. My use-case required relatively stiff springs with a high natural frequency - it worked well if the game was running at 60 frames per second, but became unstable to the point where the velocity kept increasing without any further nudges.
I mean, it makes sense. If the time step is big enough, the force being applied will not only cancel the existing velocity, but actually accelerate the mass in the opposite direction at a higher speed than it was originally at. The problem actually makes itself worse, because both the spring force and the damping force are now going to be higher, which causes the same issue on the next physics step. This quickly diverges to speeds so high that they can no longer be represented by a Lua number. This happens both when the spring constant or damping are too high, and when the framerate drops too low to accurately simulate the spring.
There are three possible solutions to this issue:
- Limit the position and velocity so that these cases cannot occur
- Run the simulation at a constant time-step, if needed run multiple physics steps per frame
- Implement a more stable spring
With the clearly jarring behavior and unfair advantage/disadvantage (depending on the case) that the first two solutions would give players on slower computers, as well as the obvious performance implications of running more physics iterations per frame when the player's computer is already struggling, the only real solution seems to be implementing a more stable spring simulation. The question is: How?
A Physical Approach
The nice thing about a mass on a spring is that its movement is essentially deterministic. The physics behind this sort of thing are very well-understood and fall under the term "simple harmonic motion". The equation for the position of a spring in terms of time looks like this:
This is already a step in the right direction - we immediately remove the need for a constant simulation update because we can immediately calculate the position of the spring at any given time. Only trouble is, anything more complex than a single nudge can't be simulated with this equation alone, and we only get the position and no velocity.
Well, that's not entirely true. We can, for example, get the velocity at any point by calculating the derivative of that function. Because we have both a position and a velocity (or, mathematically speaking, a point and a slope), we should in principle be able to calculate v0 and an offset to t in a way that nudges the spring exactly the way we want it to.
First, let's modify the equation so that we have these offsets available to us:
I've gone ahead and replaced v0/ω with a here, because we will no longer be calculating the initial velocity, but rather the amplitude of the oscillation. I've also added a time offset b, because our simulation doesn't necessarily start at the center position. At any rate, with that we have all our variables accounted for, and we can now get the derivative of that function:
This is important, because it allows us to calculate the velocity of the mass on the spring. The amplitude of the spring a is dependent on the total energy in the system, and knowing the velocity, spring constant, and position of the mass allows us to calculate it:
We can then re-arrange the equation for the potential energy in the spring to give us displacement. Giving it the total energy in the system will give us the amplitude a:
From there, we can calculate the offset b. The damping curve the sine is multiplied by is equal to 1 when t is equal to 0, so we can safely ignore it when working backwards:
And we're almost done! At this point, we only need to calculate ω. For this we can once again turn to physics, which has an answer ready for us:
At this point, we can define the mass of the spring - I've gone ahead and just specified the mass as 1 kg, so we can just remove any terms of m from the equations:
With that, we've now reduced all relevant values to a set of equations that rigidly define the position of our spring at any moment in time after we last touch the simulation parameters. This is exactly what we want - the simulation is entirely independent of any time steps. This allows for a stable simulation even in cases where the time step is significantly longer than a single oscillation of the spring.
Result
Having reduced the entire simulation down to a set of pretty straight-forward equations, we can now implement a new spring simulation:
-- @Author: AtAltitude <[email protected]> --
-- @Description: Improved mass-on-a-spring physics simulation --
----------------------------------------------------------------
--Create new table with __index meta-field to allow access to methods from objects
local Spring = {}
Spring.__index = Spring
--Constructor to initialize
function Spring.new(K, D)
--The spring follows simple harmonic motion modelled by the function
-- f(x) = a * sin(w * (x + b)) * e^(-D*x)
--As this function is re-calculated once every physics step, we have to persist all values
local new = {
K = K;
D = D;
a = 0;
b = 0;
LastTime = 0;
}
return setmetatable(new, Spring)
end
function Spring:step()
--This calculates and returns the spring's current position and velocity
local x = os.clock() - self.LastTime
local K, D, a, b, sin, cos, exp, sqrt, asin = self.K, self.D, self.a, self.b, math.sin, math.cos, math.exp, math.sqrt, math.asin
local w = sqrt(K)
return a * sin(w * (x + b)) * exp(-D*x),
a * (w * cos(w * (x + b)) * exp(-D*x) - D * sin(w * (x + b)) * exp(-D*x))
end
function Spring:addVelocity(deltaV)
--Calculate current time since simulation was last updated
local x = os.clock() - self.LastTime
--This nudges the spring by adding deltaV to its current velocity
--First, we calculate what the current velocity actually is by computing the value of:
-- f'(x) = a * (w * cos(w * (x + b)) * e^(-D*x) - D * sin(w * (x + b)) * e^(-D*x))
--Then we add the nudge to the current velocity to calculate how fast we should actually be going
local K, D, a, b, sin, cos, exp, sqrt, asin = self.K, self.D, self.a, self.b, math.sin, math.cos, math.exp, math.sqrt, math.asin
local w = sqrt(K)
local newVelocity = deltaV + a * (w * cos(w * (x + b)) * exp(-D*x) - D * sin(w * (x + b)) * exp(-D*x))
--Next, we get the current position of the spring
local currentPos = a * sin(w * (x + b)) * exp(-D*x)
--Now that we have all information, we update the timestamp of the last update to the simulation
--We do this because adding velocity to the spring essentially gives us a new set of initial conditions, and we use those
--to start a new simulation
self.LastTime = os.clock()
--We can also calculate the total energy contained within the system, as well as the undamped amplitude
--We assumed a mass of 1 in the system
local energy = 0.5 * (K*currentPos^2 + newVelocity^2)
self.a = sqrt(2*energy/K)
--At this point, if the energy is 0, that means we have no movement in the system. The below equations would
--divide by zero in this case, so we should capture it here
if (energy == 0) then
self.a = 0
self.b = 0
return
end
--With this new value for a, we can also calculate b
local baseAngle = asin(currentPos/self.a)/w
--The actual value of b depends on velocity as well, though, not just position. If the velocity is negative,
--we have to move it outside of the [-90, 90] degree area covered by asin.
if (newVelocity >= 0) then
self.b = baseAngle
else
self.b = math.pi/w-baseAngle
end
end
return Spring