%          MATLAB routine to calculate modal characteristics, 
%          namely eigenvalues, undamped natural frequencies for a 
%          beam with given boundary conditions
%
%
%
clear all
clc
%
%          Define material properties
%
rho = 7830 ;        % beam's mass density
E = 200e9 ;         % beam's Elasticity Modulus
%
%
%          Define beam's geometric properties
%
L = 1000e-3 ;       % Length
b = 25.41e-3 ;      % Width
h = 9.54e-3 ;       % Thickness
%
%          Calculate basic properties
%
A = b*h ;              % Cross section area
I = (b*h^3)/12 ;       % Polar moment of inertia
m = rho*(b*h*L) ;      % Beam' total mass 
md = m/L;              % mass per unit of length
%
k1 = 100 ;
k2 = 100 ;             % Define spring constants for soft suspension
%
i = sqrt(-1);
 %
 %      Characteristic equation for eigenvalues 
 %
 %
 EI = E*I;                        % Define which equivalent bending stiffness to use 
 wave = sqrt((EI)/(md*L^4)) ;     % dimension factor
 %
 Alfa = [0:0.0001:20];
 %
 for p = 1: length(Alfa)
     alfa = Alfa(p);
     %
     y11 = -EI*alfa^2 ;
     y12 =  0 ;
     y13 =  EI*alfa^2 ;
     y14 =  0 ;
     %
     %
     y21 =   k1;
     y22 =  -EI*alfa^3;
     y23 =   k1 ;
     y24 =   EI*alfa^3 ;
     %
     y31 = -EI*alfa^2*cos(alfa*L)  ;
     y32 = -EI*alfa^2*sin(alfa*L) ;
     y33 =  EI*alfa^2*cosh(alfa*L) ;
     y34 =  EI*alfa^2*sinh(alfa*L) ;
     %
     y41 =  EI*alfa^3*sin(alfa*L) -  k2*cos(alfa*L)  ;
     y42 = -EI*alfa^3*cos(alfa*L) -  k2*sin(alfa*L) ;
     y43 =  EI*alfa^3*sinh(alfa*L) - k2*cosh(alfa*L) ;
     y44 =  EI*alfa^3*cosh(alfa*L) - k2*sinh(alfa*L) ;
     %
     y = [y11 y12 y13 y14 ;
          y21 y22 y23 y24 ;
          y31 y32 y33 y34 ;
          y41 y42 y43 y44] ;
      %
     % roots of the characteristic equation for eigenvalues
     z(p) = det(y) ;    % my model
     %
 end
     z = z';
     %
    u = 1 ;
%
for j = 1:length(Alfa)-1
    if z(j)*z(j+1) < 0 
        an(u) = Alfa(j) ;   % search for eigenvalues model
        u = u + 1 ;
    else
    end
end
%
%
%    Calculate the natural frequencies in rad/s and Hz
%
wmodel = (an*L).^2*wave ;
fmodel = wmodel/(2*pi);
format short e
an                         % eigenvalues
fmodel                     % natural frequencies in Hz
%
%
%      %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%  END   %%%%%%%%%%%%%%%%%%%%%%%%%%%%