%	<pre>
%	Damped Free Vibration Demo (General solution)
%	Developed by Dr. Ho Sung Lee on 8/30/98 for ME240
%	
%	--------How to use this M-file for Matlab------------------------
%	Make a copy this whole file and paste this on a new file 
%	under Matlab folder or directory
%	saving as a name of 'vib.m' and then open Matlab. 
%	Type vib to run this vibration demo with a variety of inputs.
%	-----------------------------------------------------------------
%
clear all;clf;
%
%  Get the User inputs.
%

zeta=input('Enter Damping Ratio (Zeta), try 0.2    ');
wn=input('Enter Natural Circular Frequency (wn) in rad/s, try 30   ');
x0=input('Enter Initial Condition (Xo) in meter, try 0   ');
v0=input('Enter Initial Condition (Vo) in meter/sec, try 1   ');
t_max=input('Enter the max. time (sec) for plotting, try 1    ');

wd=wn*sqrt(abs(zeta.^2-1));
lamda1=-zeta*wn+wd;
lamda2=-zeta*wn-wd;
%
N=1000;	%No. of data point
%
t=0:t_max/(N-1):t_max;	%build an array for time
%
if zeta<1	%Subcritical Damping
    A=(v0+zeta*wn*x0)/wd;
    x=exp(-zeta*wn*t).*(A*sin(wd*t)+x0*cos(wd*t)); %displacement
    xdot=-zeta*wn*exp(-zeta*wn*t).*(A*sin(wd*t)+x0*cos(wd*t))...
         +exp(-zeta*wn*t).*(wd*A*cos(wd*t)-wd*x0.*sin(wd*t)); %velocity
elseif zeta==1	%Critical Damping
    
    x=x0*exp(-wn*t)+(v0+wn*x0)*t.*exp(-wn*t);; %displacement
    
    xdot=-wn*x0*exp(-wn*t)+B*exp(-wn*t).*(1-wn*t);%velocity
elseif zeta>1	%Supercritical Damping
    A=(v0-x0*lamda2)/(lamda1-lamda2);
    B=(x0*lamda1-v0)/(lamda1-lamda2);
    x=A*exp(lamda1*t)+B*exp(lamda2*t);  %displacement
    xdot=A*lamda1*exp(lamda1*t)+B*lamda2*exp(lamda2*t); %velocity
end
%    

subplot(211),plot(t,x);xlabel('time(sec)');ylabel('Displacement (m)');
title(['Damped Free System; Inputs:' ' zeta=' num2str(zeta) ... 
       ', wn=' num2str(wn) ', x0=' num2str(x0) ', v0='num2str(v0) ]);
subplot(212),plot(t,xdot),xlabel('time(sec)');ylabel('Velocity (m/s)');
	
  </pre>