%% FREQUENZE PROPRIE - modello analitico-numerico
% Meccanica dei continui per frequenze proprie energy harvester
% Ricerca delle frequenze proprie e modi di vibrare di una trave
% vincolo come cantilever, frequenze al variare della sporgenza
% CONFIGURAZIONE LIBERA
% definire param geometrici-meccanici e CC!

close all
clear all
clc

%% Definizione dei paramtri della trave
%geometria [m]
b=38.10/1000;
h=0.61/1000;

L_clamp=46.36/1000;

A=b*h;
J=b*h^3/12;

%materiale [N kg m]  QUESTI DATI NON SONO FORNITI!
E=24.0e9;
rho=3500;

EJ=E*J;
m=rho*A;
L=L_clamp;


%% Soluzione analitica
syms x A B C D g L

% soluzione generica eq della trave
w=sym(A*cos(g*x)+B*sin(g*x)+C*cosh(g*x)+D*sinh(g*x));

dw=diff(w,x);
d2w=diff(dw,x);
d3w=diff(d2w,x);

%% condizioni al contorno
cc1=subs(w,x,0);
cc2=subs(dw,x,0);
cc3=subs(d2w,x,L);
cc4=subs(d3w,x,L);

cc=[cc1; cc2; cc3; cc4];
F=[A B C D];

for i=1:4
    for j=1:4
        H(i,j)=diff(cc(i),F(j));
    end
end

determinante=det(H);

determ=simplify(determinante);
display('Determinante di H:')
pretty(determ)

Hs=subs(subs(H,'L*g','gammaL'),'g',1);

H_fun=matlabFunction(Hs);


%% Soluzione numerica

xx=[-20:.1:20];
figure
plot(xx,subs(subs(determ,'L*g',xx),'g',1),'m')
xlim([-20 20])
ylim([-100 100])
hold on, plot(xx,0)

%primo valore di gammaL
% gammaL1=fzero(ww,[0 8])

%% Chiamata programma per soluzione numerica
% in input matrice H e parametri meccanici della trave


display('Calcolo frequenze e modi...')
pause
nmax=10/pi;

gammaL_i=calcolo_frequenze(H_fun,EJ,m,L_clamp,nmax);


%% Frequenze al variare della sporgenza

d=[40:.1:60]/1000;

n=1; %fondamentale
% omega=n^2*pi^2./L.^2*sqrt(E*J/rho/A); %trave appoggio appoggio
 
gamma_i=gammaL_i(n)./d;
omega=gamma_i.^2.*sqrt(EJ/m);

freq=omega./(2*pi); %freq in Hz

figure
plot(d,freq,'r')
hold on
plot(L_clamp,(gammaL_i(n)./L_clamp).^2.*sqrt(EJ/m)/(2*pi),'gp','linewidth',2)
grid on
title('Frequenze proprie cantilever al variare della sporgenza');
xlabel('L [m]')
ylabel('f [Hz]')