%XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX
%X              Exercise 2: Frequency Response Function                   X
%XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX
 
% use international unit system

close all
clear all
clc

scrsz = get(0,'ScreenSize');

% Data---------------------------------------------------------------------

Href=1e-6;                  % Reference used in dB representation
Omax=3000;                  % 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.1;                    % 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
  
%-------------------------------------------------------------------------- 
%                        Linear Bode diagram 
%--------------------------------------------------------------------------
figure('OuterPosition',[0 scrsz(4)*12/20 scrsz(3)/4 2*scrsz(4)/5])
figure(1)
 subplot(2,1,1)              
 plot(omega/(2*pi),abs(H))
 xlabel('frequency [Hz]')
 ylabel('Receptance [m/N]')
 title('Spectrum')
 axis([0 10*freqn 0 2e-5])
 
 subplot(2,1,2)
 plot(omega/(2*pi),angle(H))
 xlabel('frequency [Hz]')
 ylabel('Phase angle [rad]')
 axis([0 10*freqn -pi 0]) 

%-------------------------------------------------------------------------- 
%                             dB Bode diagram 
%--------------------------------------------------------------------------
figure('OuterPosition',[scrsz(3)/4 scrsz(4)*12/20 scrsz(3)/4 2*scrsz(4)/5])
figure(2)                   
 subplot(2,1,1)
  plot(omega/(2*pi),20*log10(abs(H)/Href))
  xlabel('frequency [Hz]')
  ylabel('Receptance [dB]-ref=1m/N ')
  title('Spectrum')
  axis([0 10*freqn -30 30])
 
 subplot(2,1,2)
  plot(omega/(2*pi),angle(H))
  xlabel('frequency [Hz]')
  ylabel('Phase angle [N]')
  axis([0 10*freqn -pi 0])

%--------------------------------------------------------------------------
%                          Complete Bode diagram
%--------------------------------------------------------------------------

H1=1i*omega.*H;                         % mobilidade
 
H2=-omega.^2.*H;                        % inertância


figure('OuterPosition',[scrsz(3)*0 scrsz(4)/20 scrsz(3)/2 scrsz(4)/2])
figure(3)
 subplot(2,3,1)
  plot(omega/(2*pi),abs(H)) 
  ylabel('FRF Spectrum')
  axis([0 10*freqn 0 2e-5])
 
 subplot(2,3,2)
  plot(omega/(2*pi),abs(H1))
  title('Spectrum')
  axis([0 10*freqn 0 4e-3])
 
 subplot(2,3,3)
  plot(omega/(2*pi),abs(H2))
  xlabel('frequency [Hz]')
  axis([0 10*freqn 0 0.7])
 
 subplot(2,3,4)
  plot(omega/(2*pi),angle(H))
  ylabel('Phase angle [rad]')
  axis([0 10*freqn -pi pi])
  
 subplot(2,3,5)
  plot(omega/(2*pi),angle(H1))
  axis([0 10*freqn -pi pi])
 
 subplot(2,3,6)
  plot(omega/(2*pi),angle(H2))
  xlabel('frequency [Hz]')
  axis([0 10*freqn -pi pi])

%--------------------------------------------------------------------------
%                              Diagrama de Nyquist
%--------------------------------------------------------------------------

figure('OuterPosition',[scrsz(3)/2 scrsz(4)*12/20 scrsz(3)/2 scrsz(4)*4/10])
figure(4)
 subplot(1,3,1)
  plot(real(H),imag(H))
  xlabel('Real(H)')
  ylabel('Imag(H)')
  axis equal
 
 subplot(1,3,2)
  plot(real(H1),imag(H1))
  xlabel('Real(Y)')
  ylabel('Imag(Y)')
  axis equal
 
 subplot(1,3,3)
  plot(real(H2),imag(H2))
  xlabel('Real(A)')
  ylabel('Imag(A)')
  axis equal
 
%-------------------------------------------------------------------------- 
%                 Diagrama Real vs omega e Imag vs omega
%--------------------------------------------------------------------------

figure('OuterPosition',[scrsz(3)/2 scrsz(4)/20 scrsz(3)/2 scrsz(4)/2])
figure(5)
 subplot(2,3,1)
  plot(omega/(2*pi),real(H))
  ylabel('Real(FRF)')
  axis([0 10*freqn -1e-5 1e-5])
 
 subplot(2,3,2)
  plot(omega/(2*pi),real(H1))
  xlabel('frequency [Hz]')
  axis([0 10*freqn 0 4e-3])
 
 subplot(2,3,3)
  plot(omega/(2*pi),real(H2))
  axis([0 10*freqn -0.3 0.3])
  
 subplot(2,3,4)
  plot(omega/(2*pi),imag(H))
  ylabel('Imag(FRF)')
  axis([0 10*freqn -2e-5 0])
 
 subplot(2,3,5)
  plot(omega/(2*pi),imag(H1))
  xlabel('frequency [Hz]')
  axis([0 10*freqn -2e-3 2e-3])
 
 subplot(2,3,6)
  plot(omega/(2*pi),imag(H2))
  axis([0 10*freqn 0 0.6])
 
 
 
 
 
 
 