# How To Program a fully controllable spring (Damped) motion

i'm working on a project that needs realistic 2D Animation, none linear motion you might say, i have tried many things such as Hermite Curve and a few other motions and in the end Spring Motion.
because my spring is based on a simple distance, velocity, friction base i can not control the time that spring finishes and the number of times it pases through it's target, after some study i found out that i need to program damping motion and i found an article
http://www.ncsu.edu/per/Articles/MarchewkaAbbott&Beichner.pdf
which uses Total energy to solve the problem, but even after reading it many times i couldn't program it !
which is the reason i'm here , please i really need help
my goals : programming a fully controllable damping motion that is
1- has a progress based on an interval between 0 and 1 (meaning whatever the motion is, you can play it like a movie, you give the percent and it returns the position)
2- you can determine the number of times that object passes through the target before it settles down here is the final function that i have in mind :

 /*Returns the position object must be at the specified percent */
springAtPercent(startPosition, targetPosition, cycleNum /*the Number Of Times That Object Passes Through The Target*/, percent /*between 0.0 and 1.0*/);


i think the function makes it prettly clear what i need, i have searched a lot and found lots of result ! but none of them is good enough or has enough maneuver capability to achieve complicated animations that is why i need something like this.
in the end here is an image i made for a question in StackExchange physics section which resulted in the article i mentioned, hope it helps to further clarify my question. (i'm programming in pure C++, i also use ogre, in the end programming language doesn't matter that much, so if you don't know c++ feel free to use any language you like)
Thanks.

Edit 1:
i have tried to implant MickLH answer
but i have a few problem!
i can not get the friction right ! i know what Newton's Method is and how it works but how many times do i need to calculate the new friction ?
and when you say friction = 0 do you mean the x(n+1) = x + f(x) / f'(x) with n = 0 , should have an x = 0 ? because i tried that and the result of friction is very small! (as you can see the code)
i think the stiffness formula is right ! because i tried with test variable(which is 7.164850514373732 from your code) and with that constant i can get a right stiffness, which means the only problem i have is with the friction ! can you explain what should happen so that i get a right friction ? (note i can not use any math library because i'm using an company engine), thx

the lines i have written for friction and stiffness :

// as you can see the very first friction is calculated using friction = 0 which means e ^ 0 = 1
ATVector2 friction = csc(Ogre::Math::PI * cycleNum) * (((2 * Ogre::Math::PI * cycleNum * Ogre::Math::Cos(Ogre::Math::PI * cycleNum)) / decayTime) - ((2 * Ogre::Math::PI * cycleNum * 1 * pull) / (decayTime * stopped)));
ATVector2 aftePow =  ATVector2(pow(e, -(decayTime * friction.x) / 2), pow(e, -(decayTime * friction.y) / 2));
ATVector2 frictionplus1Part1 = Ogre::Math::PI * cycleNum * (2 * Ogre::Math::Cos(Ogre::Math::PI *  cycleNum) * aftePow * stopped - decayTime * friction * pull - 2 * pull);
ATVector2 frictionplus1Part2 = decayTime * (Ogre::Math::Sin(Ogre::Math::PI * cycleNum) * aftePow * stopped - (Ogre::Math::PI * cycleNum * pull));
ATVector2 frictionplus1 = frictionplus1Part1 / frictionplus1Part2;
ATVector2 test = ATVector2(1.842068074395237, 1.842068074395237);
ATVector2 stiffness = ((decayTime * decayTime) * (test * test) + 4 * (Ogre::Math::PI * Ogre::Math::PI) * (cycleNum * cycleNum))
/ (4 * decayTime * decayTime);


there is also csc which is 1 / sin(pi *4) which is divide on 0 , which if we fix, can be a large number about 672, is this okay as well ?
here is another test which even after 100 time calculating friction i can only get 1.#INF0000 which is invalid!
here is the test :

//here, the initial value for friction is set to zero
mFriction = ATVector2::ZERO;
ATVector2 aftePow = ATVector2::ZERO;
/*
decayTime = 5;
cycleNum = 4;
stopped = 0.01;
pull = ATVector(1, 1);
*/

for (int i = 0; i < 100; i++)
{

// here i calculate e^ ... to make things simpler
aftePow =  ATVector2(pow(e, -(decayTime * mFriction.x) / 2), pow(e, -(decayTime * mFriction.y) / 2));
mFriction = csc(Ogre::Math::PI * cycleNum) * (((2 * Ogre::Math::PI * cycleNum * Ogre::Math::Cos(Ogre::Math::PI * cycleNum)) / decayTime) - ((2 * Ogre::Math::PI * cycleNum * aftePow * pull) / (decayTime * stopped)));
}

//after the loop mFriction = ATVector(1.#INF0000, 1.#INF0000)


Edit 2 :
for now the friction is not going to be an issue! because i found an even bigger problem!
in MickLH answer, the friction is calculated on the onUpdate() and it uses the previous friction and that makes it useless inside a game with delta time on each update, why?
for example when decayTime is 5 it means if the animation is played for example 10 time in every second, it will play a smooth and good spring! but if some kind of lag occurs that means the friction was not updated correctly! and the animation will be played on decayTime + lagTime !
it also take away control of the element for example :
you can not play an specific part of the animation, for example in a 5 sec animation you only need the last 2 second ! using MickLH answer you can only do that if you write some kind of loop and do the calculation 10 * 3 = 30 time ! so the friction is on the right amount !
in the end The Answer does not meet my need for getting the position on a specific percent ! , which is also something i mentioned at the beginning of my question
so once again if possible, anyone, please i need a function that can return a position based on Percent, based on DeltaTime, based on the time that passes, something accurate , i say this once again, it must become something like
this :

springAtPercent(startPosition, targetPosition, cycleNum /*the Number Of Times That Object Passes Through The Target*/, percent /*between 0.0 and 1.0*/);


MickLH is it possible to write your answer based on DeltaTime, percent ?
Thx

• What level maths do you know? – MickLH Aug 16 '15 at 15:30
• @MickLH, i'm a University student at Software Engineering major, i'm not that good so if its not something super advance i should be able to understand it. – ali ahmadi Aug 16 '15 at 15:43

I'm using Hooke's Law here as the definition of a spring. ($$\F = -k * X\$$)

Given the derivatives of position and velocity, are velocity and force respectively, we can construct a differential equation for the stretching of the spring.

$$\frac{d^2}{dt^2}y(t)=-\textit{Stiffness }y(t) - \mathit{Friction} \bigg{(}\frac{d}{dt}y(t)\bigg{)}$$

Which is just a damped harmonic oscillator, and since we already know that only the under-damped case need analysis, we can obtain a nice solution:

$$\small{ y(t) = e^{-\frac{\textit{Friction }t}{2}}\bigg{(}\%k_1 \text{ sin}\big(\frac{t\sqrt{4\textit{ Stiffness} - Friction^2}}{2}\big) + \%k_2 \text{ cos} \big( \frac{t\sqrt{4 \textit{ Stiffness} - \textit{Friction}^2}}{2} \big) \bigg)}$$

We can eliminate the unknowns by knowing initial position $$\ y(0) = pull \$$ and velocity $$\ \frac{d}{dt}y(0)=0 \$$

$$\small{ y(t) = \frac { pull \sqrt{4\:\textit{Stiffness} - \textit{Friction}^2}} { \scriptsize{ \sqrt{4\:\text{Stiffness} - \textit{Friction}^2}\:e^\frac{Friction\:t}{2} \text{cos}\big(\frac{\sqrt{4\:\text{Stiffness} - \text{Friction}^2} t}{2}\big) - Friction\;e^\frac{Friction\;t}{2} \text{sin} \big( \frac{\sqrt{4\:\text{Stiffness} - \text{Friction}^2} t}{2} \big) }}}$$

Ugly. But we have a closed form solution in time :)

Before we can solve for the parameters you're interested in, we have to address one small ambiguity: The damped harmonic oscillator never stops, only decays. I will use a threshold where we consider motion "stopped", and solve for the peak which attains this amplitude.

Now, from the solution above, I have obtained that the set of peaks are generated by: Which conveniently brings the n-th peak into the relationship.

Since you've specified both the decay time, and the number of crossings, we can derive a relationship between the two constants we are trying to compute: And now by substituting the time of the peak into the amplitude function, we can express the magnitude of the peak after a target number of equilibrium crossings: Setting this equal to our "Stopped" threshold, we can then solve for the desired Friction coefficient: I have obtained an implicit solution, and although it looks like it could maybe be solved by the W function in closed form, we would then need to implement the W function in C++, so I have not investigated that. Instead, I have taken the approach of reformulating the problem as root-finding, and numerically approximating the soluiton using Newton's Method. This strategy has yielded the following fixed point iteration: I have not analysed the convergence properties of this formula with respect to the initial guess. With that said, seeding with Friction = 0 appears to place the attractor into the correct basin.

As a quick test, I've plugged in these values:

CycleNum = 4
DecayTime = 5
Stopped = 0.01
Pull = 1


Stiffness = 7.164850514373732
Friction = 1.842068074395237


Which I have then plugged into the most basic Euler integrator, with a step size of 1/10, and it produced this graph: Which appears to be right on target, crossing 4 times to achieve a peak value of 0.01, after 5 seconds of simulated time.

EDIT

Here's my prototype code, I let maxima do newton's method:

params:[DecayTime = 5.0, Stopped = 0.01, Pull = 1.0, CycleNum = 4.0]$load("mnewton")$
Friction = csc(%pi*CycleNum)*(2*%pi*CycleNum*cos(%pi*CycleNum)/DecayTime
-2*%pi*CycleNum*%e^-(DecayTime*Friction/2)*Pull
/(DecayTime*Stopped))$subst(params, %)$
mnewton(%, Friction, 0);


[[Friction=1.842068074395236]]

I also continued with maxima while completing stiffness:

subst(flatten([%,params]), Stiffness = (DecayTime^2*Friction^2+4*%pi^2*CycleNum^2)/(4*DecayTime^2)), numer;


Stiffness=7.164850514373732

So I plugged these constants into a quick little C++ program:

#include <stdio.h>

class SpringAndMassEuler {
public:
float Stiffness, Friction;
float MassOffset, MassVelocity;

SpringAndMassEuler() {
Stiffness = 7.164850514373732;
Friction = 1.842068074395237;
MassOffset = 1;
MassVelocity = 0;
}

float Update(float dt) {
MassVelocity -= dt * (Stiffness * MassOffset + Friction * MassVelocity);
MassOffset += dt * MassVelocity;
return MassOffset;
}
};

int main() {
float t = 0;

SpringAndMassEuler sim;
printf("[%f, %f],\n", t, sim.MassOffset);

// Choosing this timestep to play it safe with lame Euler
float dt = 0.1;

int c = 50;
while (c--) {
sim.Update(dt);
t+=dt;
printf("[%f, %f],\n", t, sim.MassOffset);
}
}


Which produced the data pairs that I copy pasted back into maxima:

plot2d([discrete, [ ..... ]]);

• i will try to implant this by tomorrow , hope it works, for now it looks like something that can work :) so Thank You. – ali ahmadi Aug 16 '15 at 18:37
• i have a few question 1- where should i use friction [n+1] ? 2- in friction equation , you have used friction ! when u say pow(e, (-decayTime * fiction) / 2) 3- for me to get the current position , i should get y(t) right ? i mean i give t , the passedTime , right ? – ali ahmadi Aug 17 '15 at 13:05
• if possible, can you copy ur code here ? i cant get the friction right, if possible i'd like to analyise your code to see how its implanted, thx – ali ahmadi Aug 17 '15 at 13:53
• please check my "Edit" in the main post, i still have some problem – ali ahmadi Aug 18 '15 at 14:16
• is it possible that because you are using maxima , somehow it fixes the divide on zeros and ignores them ? and now that i'm implanting them on c++ its causing problem ?, i know you have spend quite some time on your answer and that is exactly why i m asking you to help me one more time so i can get some results, please check my main post, i just cant get the friction right, Thank you for your help. – ali ahmadi Aug 20 '15 at 9:15