% Problems 15-21
%script for graphing the solution to a mass-spring-dashpot
% system and its envelope; this version treats all cases 
function msd_system(m,c,k,t0,t1,x0,v0)
%m mass; %k springconstant; c damping; t0/t1 start/end time
%x0 initial displacement; v0 initial velocity
t = linspace(t0,t1,100);
p = c/(2*m);
omega0 = (k/m)^.5;
disp(sprintf(' omega0 = %5.3e  p = %5.3e ',omega0,p))
% undamped analog
      U1 = x0;
      U2 = v0/omega0;
      AU = (U1^2+U2^2)^.5;
      deltau = atan(U2/U1);
      xu = AU*cos(omega0*t-deltau);
if omega0 >= p % decide if it is underdamped 
   omega = ((k/m)-p^2)^.5;
   disp(sprintf(' omega = %5.3e ',omega))
   if omega0 == p then % this means it is critically damped
      disp('system is critically damped')
      % solution is x(t) = exp(-p*t)*(C1 + C2*t)
      C1 = x0;
      C2 = v0+p*x0;
      x = exp(-p*t).*(C1 + C2*t);
      plot(t,x,t,xu,'k')
      xlabel(' time t')
      ylabel('x(t)')
   else
      disp('system is underdamped')
      % Now x(t) = exp(-pt)*(C1*cos(omega*t)+C2*sin(omega*t))
      T = 2*pi/omega;
      C1 = x0;
      C2 = (v0+C1*p)/omega;
      A=(C1^2+C2^2)^.5;
      delta = atan(C2/C1);
      disp(sprintf('pseudoperiod T = %5.3e Amplitude = %5.3e  delta = %5.3e ',...
               T,A,delta))
      y1 = A*exp(-p*t);
      y2 = -A*exp(-p*t);
      x =y1.*(cos(omega*t-delta));
      plot(t,x,t,y1,'.',t,y2,'.',t,xu,'k')
      xlabel(' time t')
      ylabel('x(t)')
   end
else
   disp('system is overdamped')
   % x(t) = C1*exp(p1*t) + C2*exp(p2*t)
   pp = (p^2-omega0^2)^.5;
   p1 = -p + pp; p2 = -p-pp;
   det = -2*pp;
   C1 = (x0*p2-v0)/det;
   C2 = (v0-x0*p1)/det;
   disp(sprintf('p1 = %5.3e  p2 = %5.3e  C1 = %5.3e  C2 = %5.3e',p1,p2,C1,C2))
   x = C1*exp(p1*t)+C2*exp(p2*t);
   plot(t,x,t,xu,'k')
   xlabel(' time t')
   ylabel('x(t)')
   title('overdamped case')
end
   