% Exercise 3: Response of the one degree of freedom system when 
%             it is hamonically solicited
clear all
clc
%
% System data
%
massa=10;                    % equivalent mass of the system
ksi=0.005;                  % damping ratio    
k=1e5;                      % 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
omega=[omegan/3  omegan  3*omegan];           % excitation frequencies in [rad/s]  -- only first 3 are in figure 1
%
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 it=1:length(t)
    f1(it)=F*exp(1i*omega(1)*t(it))+conj(F)*exp(-1i*omega(1)*t(it));
    f2(it)=F*exp(1i*omega(2)*t(it))+conj(F)*exp(-1i*omega(2)*t(it));
    f3(it)=F*exp(1i*omega(3)*t(it))+conj(F)*exp(-1i*omega(3)*t(it));     
    x1(it)=X(1)*exp(1i*omega(1)*t(it))+conj(X(1))*exp(-1i*omega(1)*t(it));
    x2(it)=X(2)*exp(1i*omega(2)*t(it))+conj(X(2))*exp(-1i*omega(2)*t(it));
    x3(it)=X(3)*exp(1i*omega(3)*t(it))+conj(X(3))*exp(-1i*omega(3)*t(it));
end
%
%
% Results
%
figure(1)
subplot(2,3,1)
plot(t,f1)
axis([0 0.5 -3 3])
ylabel('Excitation [EU]')


subplot(2,3,4)
plot(t,x1)
axis([0 0.5 -5e-5 5e-5])
ylabel('Response [EU]')

subplot(2,3,2)
plot(t,f2)
axis([0 0.5 -3 3])

subplot(2,3,5)
plot(t,x2)
axis([0 0.5 -5e-5 5e-5])

subplot(2,3,3)
plot(t,f3)
axis([0 0.5 -3 3])
xlabel('Time [s]')

subplot(2,3,6)
plot(t,x3)
axis([0 0.5 -5e-5 5e-5])
xlabel('Time [s]')

% 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 
%
figure(2)
subplot(3,1,1)
plot(ome,20*log10(abs(HH)))
ylabel('|Receptance|')
title(' ref=1 [m/N]')
%
FF=zeros(1,length(ome));            % bulding a vector force 
XX=zeros(1,length(ome));            % bulding a vector displacement

 for n=1:length(omega)                         % Excitation and response for each frequency                
omegaext=omega(n);                      
                                   
FFF=(A-1i*B)/2;                                % Complex amplitude of the excitation

HHH=1./(-massa*omegaext.^2+1i*c*omegaext+k);   % Value of the FRF for the 
                                               % excitation frequency    
XXX=HHH*FFF;                                   % Response

for ii=1:length(ome)
    if abs(omegaext-ome(ii))<=Dome 
           FF(ii)=FFF;
           XX(ii)=XXX; 
   end    
end

end

subplot(3,1,2)
plot(ome,abs(FF))
ylabel('|Excitation| [N]')
axis([0 400 0 1.1])

subplot(3,1,3)
plot(ome,abs(XX))               
ylabel('|Response| [m]')
xlabel('frequency [rad/s]')
%alteramos os eixos na figura 1 e o eixo y da excitação na figura 2

%f3 é uma função harmônica?
