%% 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 CON MASSA ALL'ESTREMITA'
% 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;

% massa concentrata (parallelepipedo)
Mc=15.6e-3; %[kg]
l_Mc=16e-3; %lunghezza massa [m];
s_Mc=1.5e-3; % spessore massa [m]
Jc=1/12*Mc*(l_Mc^2+4*s_Mc^2) ;%[kg]*[m^2]

L1=35.55e-3; %[m] %posizione massa
L2=L-L1; 

%% Soluzione analitica
syms x1 x2 g om t F11 F12 F13 F14 F21 F22 F23 F24 l1 l2 ej ra mc jc

% soluzione generica eq della trave, 2 domini
w1=sym((F11*cos(g*x1)+F12*sin(g*x1)+F13*cosh(g*x1)+F14*sinh(g*x1))*exp(1j*om*t));
w2=sym((F21*cos(g*x2)+F22*sin(g*x2)+F23*cosh(g*x2)+F24*sinh(g*x2))*exp(1j*om*t));

dw1_dx=diff(w1,x1);
d2w1_dx=diff(dw1_dx,x1);
d3w1_dx=diff(d2w1_dx,x1);

dw2_dx=diff(w2,x2);
d2w2_dx=diff(dw2_dx,x2);
d3w2_dx=diff(d2w2_dx,x2);

dw2_dt=diff(w2,t);
d2w2_dt=diff(dw2_dt,t);
d3w2_dtdx=diff(d2w2_dt,x2);
%% condizioni al contorno
cc1=subs(w1,x1,0);
cc2=subs(dw1_dx,x1,0);
cc3=subs(w1,x1,l1)-subs(w2,x2,0);
cc4=subs(dw1_dx,x1,l1)-subs(dw2_dx,x2,0);
cc5=subs(d2w2_dx,x2,l2);
cc6=subs(d3w2_dx,x2,l2);
cc7=subs(d2w2_dx,x2,0)-subs(d2w1_dx,x1,l1)-jc/ej*subs(d3w2_dtdx,x2,0); %equilibrio alla rotazione
cc8=subs(d3w1_dx,x1,l1)-subs(d3w2_dx,x2,0)-mc/ej*subs(d2w2_dt,x2,0); %equilibrio alla traslazione

cc=[cc1; cc2; cc3; cc4; cc5; cc6; cc7; cc8];
cc=cc/exp(om*t*1i);
cc=simplify(cc);

F=[F11 F12 F13 F14 F21 F22 F23 F24];

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

determinante=det(H);

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


H_g=subs(H,'om^2','g^4*ej/ra');
H_g=subs(H_g,'ej',EJ);
H_g=subs(H_g,'ra',m);
H_g=subs(H_g,'l1',L1);
H_g=subs(H_g,'l2',L2);
H_g=subs(H_g,'mc',Mc);
H_g=subs(H_g,'jc',Jc);

determ_g=det(H_g);

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

H_fun=matlabFunction(H_g);
det_fun=matlabFunction(determ_g);

%% Soluzione numerica

% xx=[-200:.1:200];
% figure
% plot(xx,det_fun(xx),'m')
% % xlim([-20 20])
% ylim([-1000 1000])
% hold on, plot(xx,0)

% % primo valore di gammaL
% gamma1=fzero(det_fun,[1 100])

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


display('Calcolo frequenze e modi...')

nmax=200/pi;

gamma_i=calcolo_frequenze_massa(H_fun,EJ,m,nmax,[L1 L2]);


omega=gamma_i.^2.*sqrt(EJ/m);
freq=omega./(2*pi); %freq in Hz
display(' ')
display('Frequenze proprie:');
for nn=1:3
    display(['Modo ', num2str(nn), ': ', num2str(freq(nn)), ' Hz']);
end

%% Prima frequenza propria al variare della posizione della massa
display(' ')
display('Calcolo al variare della posizione della massa...');

d=[34:.5:38]/1000; %posizioni massa dall'incastro
nmax=50/pi;

% matrice solo funzione di gamma, L1 e L2
H_d=subs(H,'om^2','g^4*ej/ra');
H_d=subs(H_d,'ej',EJ);
H_d=subs(H_d,'ra',m);
H_d=subs(H_d,'mc',Mc);
H_d=subs(H_d,'jc',Jc);

for dd=1:length(d)
    
    L1_d=d(dd);
    L2_d=L_clamp-L1_d;
    
    H_dd=subs(H_d,'l1',L1_d);
    H_dd=subs(H_dd,'l2',L2_d);
    
    H_fun_d=matlabFunction(H_dd);
    
    %variante della funzione per il calcolo senza interruzioni o plottaggio di
    %grafici (no verbose mode)
    gamma_i=calcolo_frequenze_massa_silent(H_fun_d,EJ,m,nmax,[L1_d L2_d]);
    
    freq_d(dd)=1/2/pi*(gamma_i(1))^2*sqrt(EJ/m);
    
end

figure
plot(d,freq_d,'r')
grid on
title('Prima frequenza cantilever al variare della posizione della massa');
xlabel('L [m]')
ylabel('f [Hz]')