% Exercise 2_1: Frequency Response Function fo the cantilever beam
% use international unit system
clear all
clc

% Data:

Href=1e-6;                  % Reference used in dB representation
Omax=3000;                  % Maximum value of circular frequency
N=10000;                    % Discretization points in the frequeny domain

L=0.4;                      % length of the beam
E=20e10;                    % Elasticity modulus of the steel material
ro=7850;                    % density
b=40e-3;                    % Width
h=10e-3;                     % Height
%k2=1e8;
%k1=k2;
%k3=1e5;
%keq1=(k1+k2)*k3/((k1+k2)+k3);
% L=0.6;                      % length of the beam
% E=7.1e10;                    % Elasticity modulus of the steel material
% ro=2700;                    % density
% b=40e-3;                    % Width
% h=4e-3;                     % Height
mb=b*h*L*ro;                 % beam mass        
I=b*h^3/12;                  % second order moment of the area
kb=3*E*I/L^3;                % Equivalent Stiffness of the beam [N/m]
k=kb%+keq1

m=0.203*mb						% 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]')
 xlabel('frequency [Hz]')