% Sistema com 1 gdl sob ensaio de impacto

clc
clearvars
close all

% parâmetros do sistema
m=1; % kg
c=0.2; % kg/s
k=4; % N/m

% dados da força
F0=4; % N
t0=0.05; % s
I=F0*t0; % N.s

% parâmetros intermediários
omn=sqrt(k/m); % frequência natural (rad/s)
fn=omn/(2*pi); % frequência natural (Hz)
cc=2*m*omn; % coef. de amortecimento crítico (kg/s)
zeta=c/cc; % razão de amortecimento
delta=zeta*omn;
ni=omn*sqrt(1-zeta^2);
  % visualização dos parâmetros intermediários
disp(' ')
disp(['A frequência circular natural é ',num2str(omn),' rad/s.'])
disp(' ')
disp(['A razão de amortecimento é ',num2str(zeta),'.'])
disp(' ')

% % 1o. passo - Cômputo e visualização de receptância e mobilidade
%   % cômputo de receptância e mobilidade
%     % vetor de frequências
% om=logspace(0,1,512);
%     % receptância
% H=1./(((-om.^2)*m+k)+1i*om*c);
%     % mobilidade
% Hv=1i*om.*H;
%   % visualização de receptância e mobilidade
%     % receptância (módulo e argumento)
% figure(1)
% subplot(1,2,1)
% loglog(om,abs(H),'r-')
% axis([1e0 1e1 1e-2 1e1])
% xlabel('Frequência (rad/s)')
% ylabel('Módulo da Receptância (m/N)')
% grid on
% subplot(1,2,2)
% semilogx(om,angle(H),'r-')
% xlabel('Frequência (rad/s)')
% ylabel('Argumento da Receptância (rad)')
% grid on
%     % mobilidade (módulo e argumento)
% figure(2)
% subplot(1,2,1)
% loglog(om,abs(Hv),'b-')
% axis([1e0 1e1 1e-1 1e1])
% xlabel('Frequência (rad/s)')
% ylabel('Módulo da Mobilidade (m/sN)')
% grid on
% subplot(1,2,2)
% semilogx(om,angle(Hv),'b-')
% xlabel('Frequência (rad/s)')
% ylabel('Argumento da Mobilidade (rad)')
% grid on

% % 2o. passo - Vibração (deslocamento) na frequência
%   % determinação da vibração na frequência
% X=H*I;
%   % visualização da vibração na frequência (módulo e argumento)
% figure(3)
% subplot(1,2,1)
% loglog(om,abs(X),'m-')
% title('Espectro de amplitude')
% xlabel('Frequência (rad/s)')
% ylabel('Módulo da Vibração (m.s)')
% grid on
% subplot(1,2,2)
% semilogx(om,angle(X),'m-')
% title('Espectro de fase')
% xlabel('Frequência (rad/s)')
% ylabel('Argumento da Vibração (rad)')
% grid on

% % 3o. passo - Força e vibração no tempo
% tf=linspace(0,0.5,1024);
% tx=linspace(0,50,1024);
%   % determinação da força no tempo
% f=F0*((tf>=0)-(tf>=t0));
%   % determinação da vibração no tempo
% x=(I/(m*ni))*exp(-delta*tx).*sin(ni*tx);
%   % visualização de força e vibração no tempo
% figure(4)
% subplot(1,2,1)
% plot(tf,f,'r-')
% xlabel('Tempo (s)')
% ylabel('Força (N)')
% axis([0 0.5 0 5])
% grid on
% subplot(1,2,2)
% plot(tx,x,'r-')
% xlabel('Tempo (s)')
% ylabel('Vibração (m)')
% axis([0 50 -0.1 0.1])
% grid on

% % 4o. passo - Obtenção de sinais "experimentais"
%   % definições preliminares (experimento)
% fmax=15; 
%     % frequência de amostragem (Hz)
% fs=2.56*fmax;
%     % número de pontos
% N=2048;
%   % definições associadas
% dt=1/fs; % resolução no tempo (s)
% T=N*dt; % duração dos sinais experimentais (s)
% df=1/T; % resolução em frequência (Hz)
%   % vetores de frequência e tempo
% te=0:dt:(T-dt); % vetor de tempo (s)
% fe=(0:N-1)*df;  % vetor de frequência (Hz)
% ome=2*pi*fe; % (rad/s)
%   % sinais experimentais sem ruído
% fe=F0*((te>=0)-(te>=t0));
% xe=(I/(m*ni))*exp(-delta*te).*sin(ni*te);
%   % visualização
% figure(5)
% subplot(2,1,1)
% plot(te,fe,'b.')
% xlabel('Tempo (s)')
% ylabel('Força (N)')
% axis([0 50 0 5])
% grid on
% subplot(2,1,2)
% plot(te,xe,'b.')
% xlabel('Tempo (s)')
% ylabel('Vibração (m)')
% axis([0 50 -0.1 0.1])
% grid on

% % 5o. passo - Estimação de receptância e mobilidade "experimentais"
%   % estimação via transformada rápida de Fourier
% Xe=fft(xe)*(1/fs); % com correção para sinal transiente
% Fe=fft(fe)*(1/fs); % com correção para sinal transiente
% He=Xe./Fe;
% Hve=1i*ome.*He;
%   % visualização
%     % receptância (módulo e argumento)
% figure(6)
% subplot(1,2,1)
% loglog(ome,abs(He),'b.')
% axis([1e-1 1e1 1e-2 1e1])
% xlabel('Frequência (rad/s)')
% ylabel('Módulo da Receptância (m/N)')
% grid on
% hold on
% loglog(om,abs(H),'r-')
% hold off
% legend('experimental','teórico')
% subplot(1,2,2)
% semilogx(ome,unwrap(angle(He)),'b.')
% axis([1e-1 1e1 -1.2*pi +1.2*pi])
% xlabel('Frequência (rad/s)')
% ylabel('Argumento da Receptância (rad)')
% grid on
% hold on
% loglog(om,angle(H),'r-')
% hold off
% legend('experimental','teórico')
%     % mobilidade (módulo e argumento)
% figure(7)
% subplot(1,2,1)
% loglog(ome,abs(Hve),'b.')
% axis([1e-1 1e1 1e-2 1e1])
% xlabel('Frequência (rad/s)')
% ylabel('Módulo da Mobilidade (m/sN)')
% grid on
% hold on
% loglog(om,abs(Hv),'r-')
% hold off
% legend('experimental','teórico')
% subplot(1,2,2)
% semilogx(ome,unwrap(angle(Hve)),'b.')
% axis([1e-1 1e1 -1.2*pi +1.2*pi])
% xlabel('Frequência (rad/s)')
% ylabel('Argumento da Mobilidade (rad)')
% grid on
% hold on
% loglog(om,angle(Hv),'r-')
% hold off
% legend('experimental','teórico')

% % 6o. passo - estimação de frequência natural e razão de amortecimento
% % (1) rodar as linhas 199 a 206, mantendo as linhas 207 a 224 comentadas;
% % (2) registrar aproximações da frequência de ressonância e das frequências
% %     da banda de - 3 dB, a partir do gráfico;
% % (3) retirar comentários das linhas 207 a 224 e rodar, para inserir
% %     os valores obtidos em (2) na tela de comandos;
% % (4) buscar estimativas calculadas na tela de comandos.
%   % geração do gráfico
% omt=logspace(0,1,2048);
% Ht=1./(((-omt.^2)*m+k)+1i*omt*c);
% figure(8)
% plot(omt,20*log10(abs(Ht)),'r-');
% xlabel('Frequência (rad/s)')
% ylabel('Módulo da Receptância (dB)')
% grid on
%   % cômputo numérico
% omr=input('Frequência de ressonância? ');
% disp(' ')
% om1=input('Primeira frequência da banda de -3 dB? ');
% disp(' ')
% om2=input('Segunda frequência da banda de -3 dB? ');
% disp(' ')
%   % frequência natural estimada
% omn2=omr;
%   % razão de amortecimento estimada
% zeta2=(om2-om1)/(2*omr);
%   % visualização dos parâmetros
% disp(['A frequência natural estimada é igual a ',...
%       num2str(omn2),' rad/s.'])
% disp(' ')
% disp(['A razão de amortecimento estimada é igual a ',...
%       num2str(zeta2),'.'])
% disp(' ')

% % 7o. passo - Gráfico de Nyquist
% % (1) comentar as linhas 192 a 224;
% % (2) rodar as linhas 234 a 248, mantendo as linhas 249 a 264 comentadas;
% % (3) obter aproximações da frequência natural e das frequências
% %     correspondentes ao ângulo teta igual a 90 graus;
% % (4) retirar comentários das linhas 249 a 264 e rodar;
% % (5) buscar estimativas calculadas na tela de comandos.
%   % geração de gráfico
% figure(9)
% plot(real(Hv),imag(Hv),'r.')
% xlabel('Parte real da Mobilidade (m/sN)')
% ylabel('Parte imaginária da Mobilidade (m/sN)')
% axis equal
% grid on
%   % indicação de frequências
% for k=1:length(Hv)
%   if (abs(Hv(k)) >= 0.95*max(abs(Hv))) 
%     text(real(Hv(k)),imag(Hv(k)),['  ',num2str(om(k),4)],'FontSize',9)
%   elseif max(imag(Hv))-abs(imag(Hv(k))) <= 2e-2
%     text(real(Hv(k)),imag(Hv(k)),['  ',num2str(om(k),4)],'FontSize',9,...
%          'Rotation',90)
%   end
% end
%   % cômputo numérico
% omn3=input('Frequência natural? ');
% disp(' ')
% oma=input('Frequência superior (90 graus)? ');
% disp(' ')
% omb=input('Frequência inferior (90 graus)? ');
% disp(' ')
%   % razão de amortecimento estimada
% zeta3=(omb-oma)/(2*omn3);
%   % visualização dos parâmetros
% disp(['A frequência natural estimada é igual a ',...
%       num2str(omn3),' rad/s.'])
% disp(' ')
% disp(['A razão de amortecimento estimada é igual a ',...
%       num2str(zeta3),'.'])
% disp(' ')