% Sistema com 1 gdl sob excitação harmônica
% Exercício de controle de vibrações

clc;
clearvars;
close all;

% dados do sistema
m=95; % kg
c=6.2e2; % kg/s (original: 6.2e2; localizada: 1.5e4)
k=4.0e5; % N/m (original: 4.0e5; alteração: 1.5e6)

% vibração limite
Xl=100; % (1e-6)m

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

% frequência natural (em rad/s)
omn=sqrt(k/m); 
% frequência natural (em Hz)
fn=omn/(2*pi); 
% coeficiente de amortecimento crítico
cc=2*m*omn; 
% razão de amortecimento
zeta=c/cc;
% 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(' ')

% cálculo e visualização da excitação e da vibração permanente

% vetor de tempo
t=linspace(0,2,500); % s
% excitação harmônica
F=80; % amplitude, N
fe1=1; % frequência 1, Hz
fe2=10; % frequência 2, Hz
f1=F*sin(2*pi*fe1*t); % N
f2=F*sin(2*pi*fe2*t); % N
% vibração permanente
  % razões entre frequências
e1=fe1/fn;
e2=fe2/fn;
  % fatores de amplitude
He1=(1/sqrt((1-e1^2)^2+(2*zeta*e1)^2));
He2=(1/sqrt((1-e2^2)^2+(2*zeta*e2)^2));
  % amplitudes da vibração permanente
X1=(F/k)*He1; % m
X2=(F/k)*He2; % m
  % fatores de fase
tetae1=atan2((2*zeta*e1),(1-e1^2)); % rad
tetae2=atan2((2*zeta*e2),(1-e2^2)); % rad
  % vibração permanente
xp1=X1*sin(2*pi*fe1*t-tetae1);
xp2=X2*sin(2*pi*fe2*t-tetae2); 
% visualização de curvas
figure(1);
subplot(2,1,1)
plot(t,f1,'r-',t,xp1*1e6,'b-');
grid on;
xlabel('t (s)');
ylabel('f(t)(N), xp(t)(\mum)');
legend('f(t)','xp(t)')
subplot(2,1,2)
plot(t,f2,'r-',t,xp2*1e6,'b-');
grid on;
xlabel('t (s)');
ylabel('f(t)(N), xp(t)(\mum)');
legend('f(t)','xp(t)')
