% Viga engastada-livre em translação lateral com motor dc
% (modelo discretizado simples, sgdl)

clc
clearvars
close all

% primeira etapa: Abordagem teórica

% dados dimensionais e materiais da viga
ep=6.32e-3;  % espessura (m)
lg=49.8e-3;  % largura (m)
cp=602e-3;   % comprimento (m) (até o motor)
E=200e9;     % módulo de elasticidade (Pa)
ro=7850;     % densidade volumétrica (kg/m^3)

% constante de rigidez da viga
I=lg*ep^3/12;    % momento de inércia de área (m^4)
kv=3*E*I/(cp^3); % constante de rigidez (N/m)

% constante de rigidez do sistema
k=kv; % N/m

% massa da viga
mv=ep*lg*cp*ro; % kg

% massa de motor e complementos
mm=0.245; % kg

% massa do sistema
if mv > 0.1*mm
  m=mm+0.23*mv; % kg
else
  m=mm; % kg
end

% frequência natural
omn=sqrt(kv/m); % rad/s
fn=omn/(2*pi); % Hz

% razão de amortecimento
zt=0.001; % estimada por experiência

% visualização de parâmetros modais
disp(' ')
disp(['A frequência natural circular é ',...
      'igual a ',num2str(omn),' rad/s.'])
disp(' ')
disp(['A frequência natural é igual a ',...
       num2str(fn),' Hz.'])
disp(' ')
disp(['A razão de amortecimento é igual a ',...
       num2str(zt),'.'])
disp(' ')

% receptância
  % cálculo de FRF
ft=linspace(5,20,1024); % Hz
omt=2*pi*ft; % rad/s
Hbt=(1/m)*(1./((omn^2-omt.^2)+1i*(2*zt*omn*omt))); % m/N
  % geração de espectros
figure(1)
subplot(2,1,1)
plot(ft,20*log10(abs(Hbt)),'b-');
axis([5 20 -80 -10])
grid on;
title('Espectro de amplitude')
xlabel('frequência (rad/s)')
ylabel('amplitude (dB ref.1 m/N)')
subplot(2,1,2)
plot(ft,atan2(imag(Hbt),real(Hbt)),'b-');
axis([5 20 -3*pi/2 pi/2])
grid on;
title('Espectro de fase')
xlabel('frequência (rad/s)')
ylabel('fase (rad)')

% % segunda etapa: Abordagem experimental
% 
% % nomeação de arquivos de dados
% arqdt='impacto_t2.xlsx'; % dados no tempo
% arqdf='impacto_f2.xlsx'; % dados na frequência
% 
% % leitura de arquivos de dados
% dadost=xlsread(arqdt,'A3:C5566');
% dadosf=xlsread(arqdf,'C1:BZV10');
% 
% % atribuição e visualização de dados no tempo
% % (força e aceleração)
%   % atribuição
% t=dadost(:,1); % tempo (s)
% fo=dadost(:,2); % força (N)
% ac=dadost(:,3); % aceleração (m/s^2)
%   % visualização
% figure(2)
% subplot(2,1,1)
% plot(t,fo,'r-')
% axis([0 60 -20 20])
% grid on
% title('Força x Tempo')
% xlabel('tempo(s)')
% ylabel('força (N)')
% subplot(2,1,2)
% plot(t,ac,'g-')
% axis([0 60 -30 30])
% grid on
% title('Aceleração x Tempo')
% xlabel('tempo(s)')
% ylabel('aceleração (m/s^{2})')
%  
% % atribuição e visualização de dados na frequência
% % (inertância, em partes real e imaginária e em amplitude e fase)
%   % atribuição (uso de estimadores H1 e H2)
% fe=dadosf(1,:); % frequência (Hz)
% % cor1=dadosf(2,:);  % coerência 
% I1_r=dadosf(3,:);  % I1, parte real 
% I1_i=dadosf(4,:);  % I1, parte imaginária 
% I1_a=dadosf(5,:);  % I1, amplitude 
% % I1_f=dadosf(6,:);  % I1, fase 
% % I2_r=dadosf(7,:);  % I2, parte real 
% % I2_i=dadosf(8,:);  % I2, parte imaginária 
% I2_a=dadosf(9,:);  % I2, amplitude 
% % I2_f=dadosf(10,:); % I2, fase
%   % cálculo de receptância via inertância
% Ibe=I1_r+1i*I1_i; % inertância complexa, estimador H1
% ome=2*pi*fe; % frequência (rad/s)
% Hbe=-1*(Ibe./(ome.^2)); % receptância complexa
%   % visualização
% figure(3)
% subplot(2,1,1)
% plot(fe,20*log10(I1_a),'b-',fe,20*log10(I2_a),'c-')
% axis([5 20 -10 60])
% grid on
% title('Inertância - Espectro de amplitude')
% xlabel('frequência (Hz)')
% ylabel('amplitude (dB ref 1 m/(s^{2}N)')
% legend('estimador H1','estimador H2')
% subplot(2,1,2)
% plot(fe,20*log10(abs(Hbe)),'b-')
% axis([5 20 -80 -10])
% grid on
% title('Receptância - Espectro de amplitude')
% xlabel('frequência (Hz)')
% ylabel('amplitude (dB ref 1 m/N)')

% % terceira etapa: Comparação entre abordagens
% 
% % comparação entre receptâncias: teórica x experimental 
% figure(4)
% plot(ft,20*log10(abs(Hbt)),'b-',fe,20*log10(abs(Hbe)),'c-')
% axis([5 20 -80 -10])
% grid on
% title('Receptância - Espectro de amplitude')
% xlabel('frequência (Hz)')
% ylabel('amplitude (dB ref 1 m/N)')
% legend('teórica','experimental')

% % quarta etapa: Atualização de modelo
% 
% % regeneração de receptância (manual)
% fr=ft; % Hz
% omr=omt; % rad/s
% mr=m; % kg
% fnr=10.23; % Hz
% omnr=2*pi*fnr; % rad/s
% ztr=0.0044;
% Hbr=(1/mr)*(1./((omnr^2-omr.^2)+1i*(2*ztr*omn*omr))); % m/N
% 
% % comparação entre receptâncias: regenerada x experimental
% figure(5)
% plot(fr,20*log10(abs(Hbr)),'b-',fe,20*log10(abs(Hbe)),'c-')
% axis([5 20 -80 -10])
% grid on
% title('Receptância - Espectro de amplitude')
% xlabel('frequência (Hz)')
% ylabel('amplitude (dB ref 1 m/N)')
% legend('regenerada','experimental')