%
%	Forced Damped Vibrations (General solution)
%	Developed by Dr. Ho Sung Lee on 10/20/98 for ME240
%
%   This program solves the governing equation of the forced damped 
%   vibrations with the exciting force of F(t)=Fo*sin wt for the 
%   position as a function of time with the arbitrary initial conditions 
%   (x0 and v0).
%       
%       
clear all;clf;
%
%  Get the User inputs.
%
m=input('Enter Mass (m) in kg or slug, try 2       ');
k=input('Enter Spring Constant (k) in N/m or lb/ft, try 8      ');
zeta=input('Enter Damping Ratio (Zeta), try 0.125      ');
F0=input('Enter Excitation Force (Fo) in N or lb, try 20      ');
w=input('Enter Excitation Frequency (w) in rad/s, try 4      ');
x0=input('Enter Initial Displacement (Xo) in meter or ft, try 0      ');
v0=input('Enter Initial Velocity (Vo) in m/s or ft/s, try 0      ');
t_max=input('Enter the max. time (sec) for plotting, try 25      ');

wn=sqrt(k/m);%natural circular frequency
wd=wn*sqrt(abs(1-zeta^2));%damped natural circular frquency
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
%
D=F0/k/sqrt((2*zeta*w/wn)^2+(1-(w/wn)^2)^2);%Amplitude of particular sol.
A1=2*zeta*wn*w;
A2=wn^2-w^2;
xp=F0/m/(A1^2+A2^2)*(A2*sin(w*t)-A1*cos(w*t));

if zeta<1	%Subcritical Damping
    D1=F0/m/((2*zeta*wn*w)^2+(wn^2-w^2)^2);
    B1=x0+D1*A1;
    B2=v0/wd+zeta*wn*x0/wd+D1*w/wd*(wn^2*(2*zeta^2-1)+w^2);
    xh=exp(-zeta*wn*t).*(B1*cos(wd*t)+B2*sin(wd*t));%homogeneous solution   
    x=xh+xp;
elseif zeta==1	%Critical Damping

    B=v0+wn*x0;
    xh=x0*exp(-wn*t)+B*t.*exp(-wn*t);
    x=xh+xp;
elseif zeta>1	%Supercritical Damping

    A=(v0-x0*lamda2)/(lamda1-lamda2);
    B=(x0*lamda1-v0)/(lamda1-lamda2);
    xh=A*exp(lamda1*t)+B*exp(lamda2*t);   
    x=xh+xp;
end
%
subplot(311),plot(t,xh);xlabel('time(sec)');ylabel('Homogeneous xh(t)');    
title(['Forced Vibrations; Input:' ' zeta='num2str(zeta) ...
     ', wn=' num2str(wn) ', x0=' num2str(x0) ', v0='num2str(v0) ...
     ', F(t)=' num2str(F0) ' sin(' num2str(w) 't)' ', D='num2str(D) ]);
%
subplot(312),plot(t,xp),xlabel('time(sec)');ylabel('Particular xp(t)');
%
subplot(313),plot(t,x);xlabel('time(sec)');ylabel('Total Displacement x(t)');