%XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX
%X    Exercise 3: Response of the one degree of freedom system when       X
%X                    it is hamonically solicited                         X
%XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX

close all
clear all
clc

% System data--------------------------------------------------------------

omega=[2 10 50];           % excitation frequency in [rad/s]

massa=10;                   % equivalent mass of the system
ksi=0.05;                 % damping ratio    
k=omega(2)^2*massa;         % equivalent stiffness
omegan=sqrt(k/massa);		% no damped natual frequency [rad/s]
c=2*massa*ksi*omegan;       % equivalent viscous damping
fn=omegan/(2*pi)            % no damped natural frequency [Hz]

% Excitation data----------------------------------------------------------

A=0;
B=2;

% Time domain data---------------------------------------------------------
N=2000;                     % number of points
detaT=0.00431;              % increment time
T=detaT*N;                  % final time

t=0:detaT:T;                % time vector

F=(A-1i*B)/2;               % Complex amplitude
H=1./(-massa*omega.^2+1i*c*omega+k);    % FRF - receptance                                         
X=H*F;                                  % Response for diferent excitation frequencies 
 
for n=1:3
    for it=1:length(t)
    f(n,it)=F*exp(1i*omega(n)*t(it))+conj(F)*exp(-1i*omega(n)*t(it));
    x(n,it)=X(n)*exp(1i*omega(n)*t(it))+conj(X(n))*exp(-1i*omega(n)*t(it));
    end
end
%--------------------------------------------------------------------------
%                              Results
%--------------------------------------------------------------------------

scrsz = get(0,'ScreenSize');
figure('OuterPosition',[0 scrsz(4)/20 scrsz(3) scrsz(4)*18/20])
figure(1)

% Ploting the Excitations

for n=1:3
subplot(5,3,n)
plot(t,f(n,:))
ylabel('Excitation [EU]')
axis([0 5 -2.5 2.5])
    if n==1, title('Low Frequency')
    elseif n==2, title('Natural Frequency')
                 xlabel('Time [s]')
    elseif n==3, title('High Frequency')
    end
end

% Ploting the Responses

for n=1:3
subplot(5,3,n+3)
plot(t,x(n,:))
    if n==1, axis([0 5 -6e-3 6e-3])
             ylabel('Response [EU]')
    elseif n==2, axis([0 5 -2.5e-2 2.5e-2])
                 xlabel('Time [s]')
    elseif n==3, axis([0 5 -8e-4 8e-4])
    end
end

%--------------------------------------------------------------------------
%  Response in the frequency domain for a given excitation frequency
%--------------------------------------------------------------------------

Dome=0.013;
ome=0:Dome:1.2*max(omega);          % vector frequency (variable)

HH=1./(-massa*ome.^2+1i*c*ome+k);   % FRF - receptance 


HHH=1./(-massa*omega.^2+1i*c*omega+k);% Value of the FRF for the 
                                                   % excitation frequency 
XXX=HHH*F;                  % Response for the excitation frequency
FF=zeros(3,length(ome));    % Bulding a vector force
XX=zeros(3,length(ome));

for n=1:3
for ii=1:length(ome)
    if abs(omega(n)-ome(ii))<=Dome 
           FF(n,ii)=F;
           XX(n,ii)=XXX(n); 
   end    
end
end

for n=1:3
 subplot(5,3,n+6)
 plot(ome,abs(HH))
 axis([0 60 0 1.5e-2])
    if n==1, ylabel('|Receptance| [m/N]')
    elseif n==2,  xlabel('frequency [rad/s]')
    end
end

for n=1:3
 subplot(5,3,n+9)
 plot(ome,abs(FF(n,:)))
  axis([0 60 0 2])
    if n==1, ylabel('|Excitation| [N]')
    elseif n==2, xlabel('frequency [rad/s]')
    end
end
 
for n=1:3 
subplot(5,3,n+12)
plot(ome,abs(XX(n,:)))
    if n==1,axis([0 60 0 3e-3])
            ylabel('|Response| [m]')
    elseif n==2,axis([0 60 0 1.5e-2])
                xlabel('frequency [rad/s]')
    elseif n==3,axis([0 60 0 2e-4])
    end
end



