Simple Harmonic Motion

An Alternative To Lerp

What is a lerp?

The easy, boring bit

Linear Interpolate (Lerp) is a simple mathematical function that is used to find a value somewhere between two inputs. Specifically, we’re looking at the use of Lerp in relation to a transform position.

There are 3 parameters in a Lerp. The first two, a & b, could also be called Start & End, Low & High or Min & Max, depending on context.

The final parameter, t, is sometimes called Time, Alpha, Delta, Blend, or Mix.

Most math libraries support a & b being floating point numbers or vectors (and therefore also colours). In all cases, t is a floating point number.

Some functions clamp t into the range of 0 -> 1, others allow any value for t.

a + (ba) * t;

It’s that simple. Whilst it’s good to understand the maths behind the function, bare in mind most engines have a math library which has the function built in.

How do we use it?

The classic method

There are many uses for Lerp. The one which we will be focusing on in this article is;

Each frame, interpolate a vector position from point a to point b by a fraction of delta time.

Specifically we’re looking at the following case:

lerp(currentPosition, targetPosition, deltaTime * strength)

By moving some fraction of the remaining distance between where we are currently, and where we would like to be, we move in progressively smaller steps towards our target, giving the appearance of smooth motion.

As we can see above, the movement eases in. Over the course of the movement, the velocity is largest on the first frame and smallest on last frame. This is because the distance between a and b is largest on the first frame, and lerp moves a fraction of the distance between a and b per frame.

We don’t have any concept of persistent velocity here, so the moment the target moves we get a sharp change in velocity followed by a smooth ease in.

Simple Harmonic Motion

The cool method.

An alternative way to move from point a to point b which takes current velocity into account, is called Simple Harmonic Motion. Since velocity it taken into account, an object in motion will not suddenly change path when a new target appears, it will alter course smoothly before attempting to reach a new equilibrium.

The Simple Harmonic Class class uses two parameters: Angular Frequency, and Damping Ratio. To oversimplify:

The angular frequency parameter affects how fast the current state will move towards the equilibrium state (A bit like the t value in a lerp).

The damping ratio parameter affects how ‘springy’ the motion will be.

We can under-damp (DR < 1) to allow the movement to overshoot and then attempt to recover. Under-damped harmonic motion is a really great way to add juice to UI, character movement, and physics simulations.

We can critically damp (DR == 1) to perfectly move towards the point as fast as the AF allows, without overshooting. A critically dampened spring is very similar to a lerp, but current velocity is taken into account leading to a smoother move.

We can over-damp (DR > 1), which ensures we never oscillate, but we might not reach equilibrium as quickly as possible. We haven’t found much use for over-damped springs in our game, but there are probably some cool ones out there!

For a more in depth explanation from someone who knows far more than us, check out Ryan Juckett

3D Position Interpolation

How many dimensions?

The maths being spring motion in a single dimension can easily be transferred across to multidimensional equations by simply isolating each dimension as its own spring, calculating all the springs, then recombining into a multidimensional vector.

The code included on this page has methods for 1D, 2D, and 3D calculations. If anyone was feeling especially clever it wouldn’t be difficult to add methods for 4D/Colour interpolation.

The original code comes from Ryan Juckett. We have translated his code into C# and made it Unity friendly. There are two stages to its use – The first is to call CalcDampedSpringMotionParams in order to convert from AF/DR into a set of coefficients which can be used in the Calculate method. This can be done every frame if you’re tweaking values, but should really be cached at the start of your application to avoid the expensive function call in the hot path.

We’ve opted to use ref parameters for current state and velocity.

using UnityEngine;

namespace DigitalSalmon {
	/// Harmonic motion allows spring dampened calculations in ND space.
	/// The Calculate methods can be used along with transform assignments to create very juicy motion.
	/// http://www.ryanjuckett.com/programming/damped-springs/
	public static class HarmonicMotion {
		//-----------------------------------------------------------------------------------------
		// Type Definitions:
		//-----------------------------------------------------------------------------------------

		public struct DampenedSpringMotionParams {
			// newPos = posPosCoef*oldPos + posVelCoef*oldVel
			public float PosPosCoef, PosVelCoef;
			// newVel = velPosCoef*oldPos + velVelCoef*oldVel
			public float VelPosCoef, VelVelCoef;
		}

		//-----------------------------------------------------------------------------------------
		// Public Methods:
		//-----------------------------------------------------------------------------------------

		public static void Calculate(ref Vector3 state, ref Vector3 velocity, Vector3 targetState, DampenedSpringMotionParams springMotionParams) {
			float velocityX = velocity.x;
			float stateX = state.x;
			float targetX = targetState.x;
			Calculate(ref stateX, ref velocityX, targetX, springMotionParams);

			float velocityY = velocity.y;
			float stateY = state.y;
			float targetY = targetState.y;
			Calculate(ref stateY, ref velocityY, targetY, springMotionParams);

			float velocityZ = velocity.z;
			float stateZ = state.z;
			float targetZ = targetState.z;
			Calculate(ref stateZ, ref velocityZ, targetZ, springMotionParams);

			velocity = new Vector3(velocityX, velocityY, velocityZ);
			state = new Vector3(stateX, stateY, stateZ);
		}

		public static void Calculate(ref Vector2 state, ref Vector2 velocity, Vector2 targetState, DampenedSpringMotionParams springMotionParams) {
			float velocityX = velocity.x;
			float stateX = state.x;
			float targetX = targetState.x;
			Calculate(ref stateX, ref velocityX, targetX, springMotionParams);

			float velocityY = velocity.y;
			float stateY = state.y;
			float targetY = targetState.y;

			Calculate(ref stateY, ref velocityY, targetY, springMotionParams);

			velocity = new Vector2(velocityX, velocityY);
			state = new Vector2(stateX, stateY);
		}

		public static void Calculate(ref float state, ref float velocity, float targetState, DampenedSpringMotionParams springMotionParams) {
			float oldPos = state - targetState; // update in equilibrium relative space
			float oldVel = velocity;

			state = oldPos * springMotionParams.PosPosCoef + oldVel * springMotionParams.PosVelCoef + targetState;
			velocity = oldPos * springMotionParams.VelPosCoef + oldVel * springMotionParams.VelVelCoef;
		}

		public static DampenedSpringMotionParams CalcDampedSpringMotionParams(float dampingRatio, float angularFrequency) {
			const float epsilon = 0.0001f;

			// force values into legal range
			if (dampingRatio < 0.0f) dampingRatio = 0.0f;
			if (angularFrequency < 0.0f) angularFrequency = 0.0f;

			// if there is no angular frequency, the spring will not move and we can
			// return identity
			if (angularFrequency < epsilon) { return new DampenedSpringMotionParams {PosPosCoef = 1.0f, PosVelCoef = 0.0f, VelPosCoef = 0.0f, VelVelCoef = 1.0f}; } if (dampingRatio > 1.0f + epsilon) {
				// over-damped
				float za = -angularFrequency * dampingRatio;
				float zb = angularFrequency * Mathf.Sqrt(dampingRatio * dampingRatio - 1.0f);
				float z1 = za - zb;
				float z2 = za + zb;

				float e1 = Mathf.Exp(z1 * Time.deltaTime);
				float e2 = Mathf.Exp(z2 * Time.deltaTime);

				float invTwoZb = 1.0f / (2.0f * zb); // = 1 / (z2 - z1)

				float e1_Over_TwoZb = e1 * invTwoZb;
				float e2_Over_TwoZb = e2 * invTwoZb;

				float z1e1_Over_TwoZb = z1 * e1_Over_TwoZb;
				float z2e2_Over_TwoZb = z2 * e2_Over_TwoZb;

				return new DampenedSpringMotionParams {
					PosPosCoef = e1_Over_TwoZb * z2 - z2e2_Over_TwoZb + e2,
					PosVelCoef = -e1_Over_TwoZb + e2_Over_TwoZb,
					VelPosCoef = (z1e1_Over_TwoZb - z2e2_Over_TwoZb + e2) * z2,
					VelVelCoef = -z1e1_Over_TwoZb + z2e2_Over_TwoZb
				};
			}
			if (dampingRatio < 1.0f - epsilon) {
				// under-damped
				float omegaZeta = angularFrequency * dampingRatio;
				float alpha = angularFrequency * Mathf.Sqrt(1.0f - dampingRatio * dampingRatio);

				float expTerm = Mathf.Exp(-omegaZeta * Time.deltaTime);
				float cosTerm = Mathf.Cos(alpha * Time.deltaTime);
				float sinTerm = Mathf.Sin(alpha * Time.deltaTime);

				float invAlpha = 1.0f / alpha;

				float expSin = expTerm * sinTerm;
				float expCos = expTerm * cosTerm;
				float expOmegaZetaSin_Over_Alpha = expTerm * omegaZeta * sinTerm * invAlpha;

				return new DampenedSpringMotionParams {
					PosPosCoef = expCos + expOmegaZetaSin_Over_Alpha,
					PosVelCoef = expSin * invAlpha,
					VelPosCoef = -expSin * alpha - omegaZeta * expOmegaZetaSin_Over_Alpha,
					VelVelCoef = expCos - expOmegaZetaSin_Over_Alpha
				};
			}
			else {
				// critically damped
				float expTerm = Mathf.Exp(-angularFrequency * Time.deltaTime);
				float timeExp = Time.deltaTime * expTerm;
				float timeExpFreq = timeExp * angularFrequency;

				return new DampenedSpringMotionParams {
					PosPosCoef = timeExpFreq + expTerm,
					PosVelCoef = timeExp,
					VelPosCoef = -angularFrequency * timeExpFreq,
					VelVelCoef = -timeExpFreq + expTerm
				};
			}
		}
	}
}