% Analytic solution of cantilevered beam
clear all;
clc
format short e

% Dimensions and material
a = 0.01;
b = 0.01;
L = 1;
x = 0:0.01:L;
S = a*b;
J = (1/12)*a*b^3;
E = 2.1e11;
ro = 7830;

% Load
F = -50; %[N]   

%% Static analysis

% Cantilevered beam deformation
y = (F/(6*E*J))*(x.^3 + 3*L*x.^2);
ymax = F*L^3/(3*E*J);
y_static = {num2str(ymax)};

% Static deformation of cantilevered beam
figure(1), clf
plot(x,y,'LineWidth',2), hold on
title('Static deformation of cantilevered beam')
xlabel('x [m]'), ylabel('y [m]')
axis([-0.1 1.1 -0.5 0.5])
grid on

% Maximal stress
sigmamax = F*L/(J/(b/2));

display('STATIC ANALISIS')
display({'Max. def. [m]' 'Max. Stress [Pa]'; 
        num2str(ymax)      num2str(sigmamax)})

%% Modal analysis

% Coeficients
b1L = 1.88; 
b2L = 4.69; 
b3L = 7.85; 

b1 = b1L/L;
b2 = b2L/L;
b3 = b3L/L;

% Natural frequencies [Hz]
f0_1 = (b1L^2)*sqrt(E*J/(ro*S*L^4))/(2*pi);
f0_2 = (b2L)^2*sqrt(E*J/(ro*S*L^4))/(2*pi);
f0_3 = (b3L)^2*sqrt(E*J/(ro*S*L^4))/(2*pi);

% 1. Mode shape
ksi1 = (cosh(b1L) + cos(b1L))/(sinh(b1L) + sin(b1L));
mode_1 = cosh(b1*x) - cos(b1*x) - ksi1*(sinh(b1*x) - sin(b1*x));
% 2. Mode shape
ksi2 = (cosh(b2L) + cos(b2L))/(sinh(b2L) + sin(b2L));
mode_2 = cosh(b2*x) - cos(b2*x) - ksi2*(sinh(b2*x) - sin(b2*x));
% 3. Mode shape
ksi3 = (cosh(b3L) + cos(b3L))/(sinh(b3L) + sin(b3L));
mode_3 = cosh(b3*x) - cos(b3*x) - ksi3*(sinh(b3*x) - sin(b3*x));

% Mode shapes of cantilevered beam
figure(2), clf
plot(x,mode_1,'LineWidth',2), hold on
plot(x,mode_2,'g','LineWidth',2), hold on
plot(x,mode_3,'r','LineWidth',2), hold on
title('Mode shapes of cantilevered beam')
legend('Mode 1','Mode 2','Mode 3')
xlabel('x [m]'), ylabel('y [m]')
axis([-0.1 1.1 -5 5])
grid on

display('MODAL ANALISIS')
display({'Mode 1 [Hz]' 'Mode 2 [Hz]' 'Mode 3 [Hz]'; 
        num2str(f0_1) num2str(f0_2) num2str(f0_3)})

%% Harmonic analisis

% Frquency vektor [rad/s]
w = 2*pi*(0:0.1:100);
% POI
x = L;

r = ((ro*S*w.^2)/(E*J)).^(1/4);
num_w = (cosh(r*L) + cos(r*L)).*(sin(r*x) - sinh(r*x)) - (sinh(r*L) + sin(r*L)).*(cos(r*x) - sinh(r*x));
den_w = 2*E*J*r.^3.*(1 + cos(r*L).*cosh(r*L));
y_w = F*num_w./den_w;

% Frquency vektor [Hz]
freq = w/2/pi;  

% Frequency response of the tip of cantilevered beam
figure(3), clf
semilogy(freq,abs(y_w),'- .'), hold on
title('Frequency response of the tip of cantilevered beam')
xlabel('freq [Hz]'), ylabel('y [m]')
grid on

%% Transient analysis
% TOO MUCH ...


