% Sistema com 1 gdl sob excitação harmônica pela base

clc;
clearvars;
close all;

% dados fixos
c=2e3; % coeficiente de amortecimento (kg/s)
k=4e4; % constante de rigidez (N/m)
Y=0.01; % amplitude do deslocamento de base (m)

% dados variáveis
m=1007; % massa do automóvel (kg) [1007 1585]
v=20; % velocidade do automóvel (km/h) [20 40 60 80 110]

% passo 1 - cálculo e visualização
%           de parâmetros do sistema

% parâmetros do sistema
omn=sqrt(k/m); % frequência natural (rad/s)
fn=omn/(2*pi); % frequência natural (Hz)
cc=2*m*omn; % coeficiente de amortecimento crítico (kg/s)
zeta=c/cc; % razão de amortecimento
% visualização da frequência natural
% e da razão de amortecimento
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(zeta),'.'])
disp(' ')

% % passo 2 - razão entre amplitude de vibração
% %           e amplitude da base (deslocamentos)
% 
% ev=logspace(-1,1,500);
% % razão entre amplitudes
% Tv=(sqrt(1+(2*zeta*ev).^2)./...
%     sqrt((1-ev.^2).^2+(2*zeta*ev).^2)); 
% % visualização da razão entre amplitudes
% figure(1)
% semilogx(ev,20*log10(Tv),'-b');
% grid on;
% xlabel('razão entre frequências \epsilon')
% ylabel('razão entre amplitudes T(\epsilon) (dB ref.1m/m)')

% % passo 3 - cálculo e visualização do deslocamento de base 
% %           e da vibração permanente
% 
% t=linspace(0,2,500);
% % deslocamento harmônico de base
% L=6; % distância entre picos (m)
% omb=(2*pi)*(1/L)*v*(1000/3600);
% y=Y*sin(omb*t);
% % vibração permanente
%   % razão entre frequências
% e=omb/omn;
%   % amplitude de vibração permanente
% X=Y*(sqrt(1+(2*zeta*e)^2)/sqrt((1-e^2)^2+(2*zeta*e)^2));
%   % fase da vibração permanente
% tetax=atan2(2*zeta*e,(1-e^2))-atan2(2*zeta*e,1);
%   % expressão com amplitude e fase
% xp=X*sin(omb*t-tetax);
% % visualização da amplitude de vibração permanente
% % e da frequência de excitação 
% disp(['A amplitude da vibração permanente é igual a ',...
%        num2str(X),' m.'])
% disp(' ')
% disp(['A frequência circular da vibração permanente é igual a ',...
%        num2str(omb),' rad/s.'])
% disp(' ')
% disp(['A razão entre frequências é igual a ',...
%       num2str(e),' .'])
% disp(' ')
% % visualização do movimento da base 
% % e da vibração permanente
% figure(2);
% plot(t,y,'r-',t,xp,'b-');
% axis([0 2 -0.04 0.04])
% grid on;
% xlabel('t (s)');
% ylabel('y(t) (m), xp(t) (m)');
% legend('y(t)','xp(t)');