clc;
clearvars;
close all;

% dados do sistema primário
m1=300; % massa (kg)
lv=1.8; % comprimento da viga (m)
Iv=1.8e-5; % momento de inércia de área da viga (m^4)
E=200e9; % módulo de elasticidade longitudinal (Pa)

% dados operacionais do sistema primário
rpm=750; % velocidade angular (rpm)
om=rpm*((2*pi)/60); % frequência de operação (rad/s)
Xor=8e-4; % amplitude permanente de vibração harmônica (m)

% rigidez da viga (N/m)
k1=(3*E*Iv)/(lv^3);

% frequência natural do sistema primário (rad/s)
om1=sqrt(k1/m1);

% projeto do neutralizador
  % massa do neutralizador (kg)
m2=30; 
  % razão entre massas
mi=m2/m1;
  % frequência natural do neutralizador
om2=om1/(1+mi);
  % constante de rigidez do neutralizador
k2=m2*om2^2;
  % coeficiente de amortecimento viscoso do neutralizador
c2=2*m2*om2*sqrt((3*mi)/(8*(1+mi)));
  % visualização dos resultados
disp(' ')
disp(['A frequência do neutralizador é ',num2str(om2),' rad/s.'])
disp(' ')
disp(['A massa do neutralizador é ',num2str(m2),' kg.'])
disp(' ')
disp(['A constante de rigidez do neutralizador é ',num2str(k2),' N/m.'])
disp(' ')
disp(['O coeficiente de amortecimento do neutralizador é ',num2str(c2),' kg/s.'])
disp(' ')

% frequências naturais do sistema composto
  % coeficientes da equação biquadrada
c(1)=1;
c(2)=0;
c(3)=-((k1/m1)+(1+(m2/m1))*(k2/m2));
c(4)=0;
c(5)=(k1*k2)/(m1*m2);
  % raízes da equação biquadrada
rs=roots(c); 
  % raízes ordenadas
rs=sort(rs); 
  % frequências naturais
omn=rs(3:4);
  % visualização das frequências naturais
disp(['As frequências naturais do sistema composto são ',num2str(omn(1)),...
      ' rad/s e ',num2str(omn(2)),' rad/s.'])
disp(' ')

% funções resposta em frequência
  % coef. de amortecimento do sistema primário (kg/s)
c1=0; % hipótese de amortecimento desprezável
  % faixa de frequência (rad/s)
omr=logspace(1,3,1000);
  % FRF do sistema primário
H=1./((-omr.^2*m1+k1)+(1i*c1*omr)); % (m/N)
  % FRFs do sistema composto
ab=-omr.^2*m1+1i*omr*(c1+c2)+(k1+k2);
bb=1i*omr*c2+k2;
cb=-omr.^2*m2+1i*omr*c2+k2;
H1=cb./(ab.*cb-bb.^2); % (m/N)
H2=bb./(ab.*cb-bb.^2); % (m/N)
  % visualização
figure(1)
semilogx(omr,20*log10(abs(H)),'r-',omr,20*log10(abs(H1)),'b-',...
         omr,20*log10(abs(H2)),'g-')
xlabel('Frequência (rad/s)');
ylabel('Magnitude de FRF (dB ref. 1m/N)');
axis([30 300 -170 -60]);
grid on
legend('SP','SC (máq.)','SC (ndv)')

% amplitude de vibração do sistema primário
  % função resposta em frequência original
Hsp=1./((-om^2*m1+k1)+(1i*c1*om)); % (m/N)
  % amplitude de excitação (N)
F=Xor/abs(Hsp);
  % função resposta em frequência modificada
ab=-om^2*m1+1i*om*(c1+c2)+(k1+k2);
bb=1i*om*c2+k2;
cb=-om^2*m2+1i*om*c2+k2;
H1sc=cb./(ab.*cb-bb.^2); % (m/N)
  % nova amplitude permanente
Xmd=F*abs(H1sc);
  % visualização da nova amplitude permanente
disp(['A nova amplitude de vibração permanente é ',num2str(Xmd),' m, ',...
      'enquanto que a amplitude original era ',num2str(Xor),' m.'])
disp(' ')