GraphingCalculator 4;
Window 46 13 855 1400;
PaneDivider 443;
FontSizes 14;
StackPanes 1;
Slider 0 100;
SliderSteps 500;
SliderControlValue 0;
2D.Scale 0.1 0.25 5 1;
2D.BottomLeft -5.703125 -2.86875;
2D.Axes 0;
2D.GraphPaper 0;
Text "Phase space trace of (under-)damped oscillator.
Version 0.71 4-24-14

With k=9 and m=1, the natural frequency W=3, so critical damping (W=d) occurs when b=6.   Change range of b to [0,1] to see very light damping. Otherwise b=[0,10]

To do:
(i) 
  ";
Color 4;
Expr k=9,m=1;
Color 5;
MathPaneSlider 1;
Expr b=slider([0,10,100]);
Color 8;
Expr d=b/(2*m);
Expr W=sqrt(k/m);
Color 7;
Expr a=sqrt(W^2-d^2);
Color 4;
Expr c=sqrt(d^2-W^2);
Color 6;
Expr A_0=1;
Text "A – undamped amplitude; k – spring constant; m – mass, b – damping constant; d – damping parameter; a – angular frequency; c – overdamped angular ""frequency""; W – natural frequency; X – position; V – velocity

Three sets of initial conditions such that all three solutions match at t=0.   

The first set has x(0)=A and phase constant p=0 but does not try to match initial velocities: A=A0, B=0,C=A,D=A,G=0,p=0
The second set has x(0)=A and phase constant p=0 while matching initial velocities, v_x(0)=-dA: A=A0, B=A/2,C=A/2,D=A,G=0,p=0
The third set matches the initial position and velocity of the undamped case, x(0)=A0 and v_x(0)=0: A=A0/cosp, B=(c-d)/2c*A0, C=(c+d)/2c*A0,D=A0,G=d*A0,p=-atan(d/a)
";
Color 2;
Expr A=A_0,B=0,C=A_0,D=A_0,G=0,p=0;
Color 3;
Expr function(X,s)=branch(if(A*e^(-(d*s))*cos([a*s+p]),d<W),if([D+G*s]*e^(-(d*s)),d=W),if([B*e^(-(c*s))+C*e^(c*s)]*e^(-(d*s)),d>W));
Color 17;
Expr function(V,s)=branch(if(-(a*A*e^(-(d*s))*sin([a*s+p]))-(d*A*e^(-(d*s))*cos([a*s+p])),d<W),if([G-(d*D)-(d*G*s)]*e^(-(d*s)),d=W),if([-(B*[c+d]*e^(-(c*s)))+C*[c-d]*e^(c*s)]*e^(-(d*s)),d>W));
Color 6;
Expr function(k,x)=X*cos([f*x-P]);
Color 7;
Expr function(l,x)=-(f*X*sin([f*x-P]));
Color 17;
Grain 1;
Expr vector(x,y)=vector(function(X,n*t)+function(k,n*t),function(V,n*t)+function(l,n*t));
Color 17;
Expr vector(x,y)=vector(function(X,n)+function(k,n),function(V,n)+function(l,n));
Text "
F0 – magnitude of driving force; f – angular frequency of driving force; X – amplitude of driven oscillation; P – phase of driven oscillation";
Color 3;
Expr F_0=slider([0,10,20]);
Color 4;
MathPaneSlider 5;
Expr f=slider([0,10,50]);
Color 2;
Expr function(F,x)=F_0*cos([f*x]);
Color 5;
Expr X=[F_0/m]/sqrt([W^2-f^2]^2+4*d^2*f^2);
Color 6;
Expr P=atan([2*d*f/(W^2-f^2)]);
Text "Energy (approximate and exact) and potential energy superimposed:";
Color 2;
Expr function(E,n)=1/2*k*function(X,n)^2+1/2*m*function(V,n)^2;
Color 8;
Expr function(U,x)=1/2*k*x^2;
Color 17;
Expr y=1/2*k*A^2;
Color 17;
Expr y=1/2*k*A^2*e^(-(2*d*n));
Color 17;
Expr y=function(E,0);
Color 17;
Expr y=function(E,n);
Color 17;
Expr function(U,x);
Color 17;
Expr vector(x,y)=vector(function(X,n)+function(k,n),1/2*k*function(X,n)^2);
Text "Turning points:";
Color 17;
Expr abs(x)=A;
Color 17;
Expr abs(x)=A,y>0;
Color 17;
Expr abs(x)=A*e^(-(d*n));
Color 17;
Expr abs(x)=A*e^(-(d*n)),y>0;
Text "Oscillating mass superimposed";
Color 5;
Expr vector(x,y)=vector(function(X,n)+function(k,n),0);
MathPaneSlider 200;
Expr h=slider([0,0.05]);
Color 7;
Expr abs(y)<h,0<x<function(X,n)+function(k,n);
Color 7;
Expr abs(y)<h,function(X,n)+function(k,n)<x<0;
Text "Driving force superimposed:";
Color 17;
Expr vector(x,y)=vector(function(F,n),0);
Color 17;
Expr abs(y)<h,0<x<function(F,n);
Color 17;
Expr abs(y)<h,function(F,n)<x<0;
Text "
Author: David A. Craig <<http://web.lemoyne.edu/~craigda/>
";