% Exercise 2: Frequency Response Function
% use international unit system
clear all
clc

% Data:

Href=1e-6;                  % Reference used in dB representation
Omax=1000;                  % Maximum value of circular frequency
N=10000;                    % Discretization points in the frequeny domain 
k=3.16e5;					% Stiffness parameter [N/m]
m=10;						% Mass of the system [kg]
ksi=0.01;                   % damping ratio = C/Cc 

omegan=sqrt(k/m);           % undamped natural frequency in [rad/s]
freqn=omegan/2/pi;          % undamped natural frequency in [Hz]
C=2*m*ksi*omegan;           % damping parameters [N.s/m]

Domega=Omax/N;
omega=0:Domega:Omax;        % frequency in [rad/s]

 H=1./(-m*omega.^2+1i*omega*C+k);        % receptance

figure(1)
 subplot(2,1,1)              % Linear Bode diagram
 plot(omega/(2*pi),abs(H))
 xlabel('frequency [Hz]')
 ylabel('Receptance [m/N]')
 title('Spectrum')
 subplot(2,1,2)
 plot(omega/(2*pi),angle(H))
 xlabel('frequency [Hz]')
 ylabel('Phase angle [rad]')

figure(2)                   % dB Bode diagram
 subplot(2,1,1)
 plot(omega/(2*pi),20*log10(abs(H)/Href))
 xlabel('frequency [Hz]')
 ylabel('Receptance [dB]-ref=1m/N ')
 title('Spectrum')
 subplot(2,1,2)
 plot(omega/(2*pi),angle(H))
 xlabel('frequency [Hz]')
 ylabel('Phase angle [N]')


 H1=1i*omega.*H;                         % mobilidade
% 
 H2=-omega.^2.*H;                        % inertância
% 
% Complete Bode diagram
% 
figure(3)
 subplot(2,3,1)
 plot(omega,abs(H)) 
 ylabel('FRF Spectrum')
 subplot(2,3,4)
 plot(omega,angle(H))
 ylabel('Phase angle [rad]')
 
 subplot(2,3,2)
 plot(omega,abs(H1))
 title('Spectrum')
 
 subplot(2,3,5)
 plot(omega,angle(H1))
 subplot(2,3,3)
 plot(omega,abs(H2))
 xlabel('frequency [Hz]')
 subplot(2,3,6)
 plot(omega,angle(H2))
 xlabel('frequency [Hz]')
% 
% Diagrama de Nyquist
% 
figure(4)
 subplot(3,1,1)
 plot(real(H),imag(H))
 xlabel('Real(H)')
 ylabel('Imag(H)')
 axis equal
 subplot(3,1,2)
 plot(real(H1),imag(H1))
 xlabel('Real(Y)')
 ylabel('Imag(Y)')
 axis equal
 subplot(3,1,3)
 plot(real(H2),imag(H2))
 xlabel('Real(A)')
 ylabel('Imag(A)')
 axis equal
% 
% Diagrama Real vs omega e Imag vs omega
% 
figure(5)
 subplot(2,3,1)
 plot(omega,real(H))
 ylabel('Real(FRF)')
 
 subplot(2,3,4)
 plot(omega,imag(H))
 ylabel('Imag(FRF)')
 subplot(2,3,2)
 plot(omega,real(H1))
 subplot(2,3,5)
 plot(omega,imag(H1))
 subplot(2,3,3)
 plot(omega,real(H2))
 xlabel('frequency [Hz]')
 subplot(2,3,6)
 plot(omega,imag(H2))
 xlabel('frequency [Hz]')