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

clc;
clearvars;
close all;

% dados do sistema
m=100; % kg
c=500; % kg/s
k=30000; % N/m

% passo 1 - 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(' ')

% % passo 2 - cálculo e visualização da
% % excitação e da vibração
% % (permanente, transiente e completa)
%
% % vetor de tempo
% t=linspace(0,3,1000);
% % excitação harmônica
% F=80; % amplitude, N
% om=10; % frequência, rad/s
% f=F*sin(om*t);
% % vibração permanente
% % razão entre frequências
% e=om/omn;
% % fator de amplitude
% He=(1/sqrt((1-e^2)^2+(2*zeta*e)^2));
% % amplitude da vibração permanente
% X=(F/k)*He;
% % fator de fase
% tetae=atan2((2*zeta*e),(1-e^2));
% % expressão da vibração permanente
% xp=X*sin(om*t-tetae);
% % vibração transiente
% % fase
% psie=atan2((2*zeta*sqrt(1-zeta^2)),(2*zeta^2-(1-e^2)));
% % expressão da vibração transiente
% xt=X*(e*exp(-zeta*omn*t)/sqrt(1-zeta^2)).*...
% sin(omn*sqrt(1-zeta^2)*t+psie);
% % vibração completa
% x=xp+xt;
% % visualização da razão entre frequências
% disp(['A razão entre frequências é igual a ',...
% num2str(e),'.'])
% disp(' ')
% % visualização de curvas
% figure(1);
% subplot(2,1,1)
% plot(t,f,'m-');
% axis([0 3 -100 100]);
% grid on;
% xlabel('t (s)');
% ylabel('f(t) (N)');
% subplot(2,1,2)
% plot(t,xp*1000,'c-',t,xt*1000,'g-',t,x*1000,'b-');
% axis([0 3 -20 20]);
% grid on;
% xlabel('t (s)');
% ylabel('xp(t), xt(t), x(t) (mm)');
% legend('xp(t)','xt(t)','x(t)');

% % passo 3 - cálculo e visualização dos fatores
% % de amplitude e fase (deslocamento)
%
% % vetor de razão entre frequências
% ev=logspace(-1,1,500);
% % fator de amplitude
% Hev=(1./sqrt((1-ev.^2).^2+(2*zeta*ev).^2));
% % fator de fase
% tetaev=atan2((2*zeta*ev),(1-ev.^2));
% % visualização:
% figure(2)
% subplot(2,1,1)
% semilogx(ev,20*log10(Hev),'r-');
% axis([0.1 10 -40 20])
% grid on;
% xlabel('razão entre frequências \epsilon')
% ylabel('fator de amplitude H(\epsilon) (dB ref.1m/N)')
% subplot(2,1,2)
% semilogx(ev,tetaev,'r-')
% grid on;
% xlabel('razão entre frequências \epsilon')
% ylabel('fator de fase \phi(\epsilon)')

% passo 4 - verificação de comportamento com alterações
% em frequência e amortecimento

% ---> alterar: om = 17 rad/s (linha 44)
% (perto da frequência natural)

% ---> alterar: om = 30 rad/s (linha 44)
% (acima da frequência natural)

% ---> alterar: om = 17 rad/s (linha 44)
% c = 1000 kg/s (linha 9)
% (perto da frequência natural)
% (inserção de amortecimento: amplitude e transiente)

% ---> alterar: c = 250 kg/s (linha 9)
% (perto da frequência natural)
% (redução de amortecimento: amplitude e transiente)

% ---> alterar: c = 100 kg/s (linha 9)
% (perto da frequência natural)
% (amortecimento baixo: variação abrupta de fase,
% ressonância quase na freq. natural,
% fase como indicador de freq. natural)

% ---> alterar: c = 500 kg/s (linha 9)
% k = 60000 N/m (linha 10)
% (perto da frequência natural)
% (dessintonização: alternativa de controle)