Underdamped Harmonic Oscillations

underdampedoscillations.m

// NOTE: I've decided to use a much simpler approach to this as the complexity was starting to consume
// too much energy. Leaving this here for historical context.

// Building an Underdamped Harmonic Oscillator
//
// All equations derived from this wonderful paper by David Morin:
// http://www.people.fas.harvard.edu/~djmorin/waves/oscillations.pdf
//
// Note on mass: In order to simplify our calculations I'll assume that mass (m) is always 1. Where
//               m would otherwise be written I've folded it into the expression.
//
// To start, we pull out two equations from the paper:
// Oscillation  (37):  tween(t) = cos(ft) + (-v/f) * sin(ft)
// Underdamping (53):  tween(t) = exp(-Dc * t / 2)
// Note: An adjustment was made to (37) of negating v - without this the velocity behaves backward to
//       what you would expect.
//
// The underdamping is an exponential decay. The oscillation is a cosine wave. By themselves they're
// not particularly intriguing, but when we combine them, we get a delightful decaying oscillation.
//
// Underdamped oscillation: exp(-Dc * t / 2) * (cos(ft) + (-v/f) * sin(ft))
//
// Definitions:
// Dc = damping coefficient
// Sc = spring coefficient
// f = oscillation frequency
// t = animation tick (valid range: [0...d])
// d = animation duration
//
// How to calculate f (oscillation frequency) (Wu in the paper)
// ------------------------------------------------------------
//
// We simplify equation (50):  f = sqrt(Sc) * sqrt(1 - pow(Dc / (2 * sqrt(Sc)), 2))
//                        to:  f = sqrt(Sc - pow(Dc / 2, 2))
//
// Note that Sc nicely cancels itself out leaving us with a fairly simple equation.
//
//
// How to calculate Dc and Sc (spring) coefficients
// ------------------------------------------------
//
// A tween needs to predictably travel from 0 to 1, or 1 to 0 over a specific duration (d). In
// this case our tween should reach 0 at d, so let's solve the Dc and Sc variables by finding when
// our two equations (37 and 53) reach 0 with a given t value.
//
// Part 1: Underdamping
// The underdamping function (53) only has one variable in it, Dc, so let's start by solving for
// that. Because it is an exponential decay we can't solve for 0 (it never reaches it!), but we can
// solve for something close. Let's call it epsilon and define it as 0.001.
//
// exp(-Dc * t / 2) = epsilon
//      -Dc * t / 2 = log(epsilon)
//            Dc(t) = -2log(epsilon) / t
//
// We can now calculate a Dc for any given t duration.
//
// Part 2: Oscillation
// Calculating when the oscillation reaches zero is trivial when we don't have an initial velocity.
//
//                   cos(ft) = 0
//                        ft = acos(0)
//                         f = acos(0) / t
//                         f = sqrt(Sc - pow(Dc / 2, 2))
// sqrt(Sc - pow(Dc / 2, 2)) = acos(0) / t
//       Sc - pow(Dc / 2, 2) = pow(acos(0) / t), 2)
//                     Sc(t) = pow(acos(0) / t), 2) + pow(Dc / 2, 2)
//
// However, when we introduce initial velocity and attempt to solve for Sc again, we run in to an
// interesting problem:
//
// cos(ft) + (-v/f) * sin(ft) = 0
//           (-v/f) * sin(ft) = -cos(ft)
//            (v/f) * sin(ft) = cos(ft)
//            (v/f) * tan(ft) = 1
//                    tan(ft) = f/v
//                         ft = atan(f/v)
//
// At this point there is no way to solve for f (that I've been able to find) other than by using
// a root-finding algorithm (Newton's method, for example).
//
//

typedef CGFloat (^CalculateEquationBlock)(CGFloat x);

CGFloat NIFindRoot(CalculateEquationBlock fn, CalculateEquationBlock fndx, CGFloat guess) {
  CGFloat root = guess;
  while (1) {
    CGFloat newRoot = root - fn(root) / fndx(root);
    NSLog(@"%f", newRoot);
    if (fabsf(newRoot - root) < 0.001) { // or whatever epsilon you want
      return newRoot;
    }
    root = newRoot;
  }
  return root;
}

void NISpringTweenValuesForDuration(CGFloat duration, CGFloat initialVelocity, CGFloat* damping, CGFloat* frequency) {
  static const CGFloat epsilon = 0.001;

  // Calculates a damping and stiffness such that an underdamped harmonic oscillation settles to 0
  // at duration.

  // Calculate a damping such that our exponential decay reaches epsilon at duration.
  *damping = 2 * -logf(epsilon) / duration;

  // Calculate a stiffness such that our harmonic oscillation reaches zero at duration.

  CGFloat stiffness = NIFindRoot(^CGFloat(CGFloat x) {
    CGFloat freq = NICGFloatSqrt(x - NICGFloatPow((*damping) / 2, 2));
    if (isnan(freq)) freq = 0;
    return NICGFloatCos(freq * duration) + -initialVelocity / freq * NICGFloatSin(freq * duration);
  }, ^CGFloat(CGFloat x) {
    CGFloat freq = NICGFloatSqrt(x - NICGFloatPow((*damping) / 2, 2));
    if (isnan(freq)) freq = 0;
    // Grabbed from WolframAlpha cuz I'm lazy as balls.
    return (((initialVelocity - duration * freq * freq) * NICGFloatSin(freq * duration)
             - duration * initialVelocity * freq * cos(freq * duration))
            / (freq * freq));
  }, NICGFloatPow(*damping / 2, 2) + 1);

	*frequency = NICGFloatSqrt(stiffness - NICGFloatPow((*damping) / 2, 2));
	if (isnan(*frequency)) *frequency = 0;
}

CGFloat NISpringTween(CGFloat t, CGFloat initialVelocity, CGFloat damping, CGFloat frequency, NSInteger numberOfBounces) {
  CGFloat bounceScalar = ((CGFloat)numberOfBounces * 2) + 1;

  CGFloat dampener = NICGFloatExp(-damping * t / 2);
  CGFloat oscillator = (NICGFloatCos(frequency * t * bounceScalar)
                        + (-initialVelocity / frequency) * NICGFloatSin(frequency * t * (bounceScalar - 1)));

	return dampener * oscillator;
}

CGFloat NIScaledSpringTween(CGFloat from, CGFloat delta, CGFloat initialVelocity, CGFloat t, CGFloat damping, CGFloat frequency, NSInteger numberOfBounces) {
  if (delta < 0) {
    initialVelocity = -initialVelocity;
  }
	return from + delta - delta * NISpringTween(t, initialVelocity, damping, frequency, numberOfBounces);
}