clc;
clearvars;
close all;

% dados do sistema primário
m1=7.3; % massa (kg)
om1=2*pi*3.45; % frequência natural (rad/s)
k1=m1*(om1)^2; % constante de rigidez (N/m)

% condições de projeto
F1=13;
omf=2*pi*3.45;
X2=0.025;

% frequência de sintonização
oms=omf;

% % passo 1 - projeto do neutralizador (NDV) não amortecido
%   % constante de rigidez do NDV
% k2=F1/X2; % (N/m)
%   % massa do NDV
% m2=k2/oms^2; % (kg)
%   % visualização dos resultados
% disp(' ')
% disp(['A frequência de sintonização inicial é ',num2str(oms),' rad/s.'])
% disp(' ')
% disp(['A massa inicial do neutralizador é ',num2str(m2),' kg.'])
% disp(' ')
% disp(['A constante de rigidez inicial do neutralizador é ',num2str(k2),' N/m.'])
% disp(' ')

% % passo 2 - 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 dos resultados
% disp(['As frequências naturais iniciais do sistema composto são ',...
%        num2str(omn(1)),' rad/s e ',num2str(omn(2)),' rad/s.'])
% disp(' ')

% % passo 3 - funções resposta em frequência individuais
%   % dados adicionais
% c1=1; % coef. de amortecimento do sistema primário
% c2=5; % coef. de amortecimento do neutralizador (0.001, 1, 2, 5)
%   % faixa de frequência
% om=logspace(1,2,1000);
%   % FRF do sistema primário
% H=1./((-om.^2*m1+k1)+(i*c1*om));
%   % FRFs do sistema composto
% ab=-om.^2*m1+i*om*(c1+c2)+(k1+k2);
% bb=i*om*c2+k2;
% cb=-om.^2*m2+i*om*c2+k2;
% H1=cb./(ab.*cb-bb.^2);
% H2=bb./(ab.*cb-bb.^2);
%   % visualização
% figure(1)
% semilogx(om,20*log10(abs(H)),'r-',om,20*log10(abs(H1)),'b-',...
%          om,20*log10(abs(H2)),'g-')
% xlabel('Frequência (rad/s)');
% ylabel('Magnitude de FRF (dB ref. 1m/N)');
% axis([10 100 -140 0]);
% grid on
% legend('sist. primário','sist. composto (máquina)','sist. composto (ndv)')

% % passo 4 - funções resposta em frequência em conjunto
%   % dados adicionais
% c1=1; % coef. de amortecimento do sistema primário
% c2=[0.001 1 2 5]; % coefs. de amortecimento do neutralizador (0.001, 1, 2, 5)
%   % faixa de frequência
% om=logspace(1,2,1000);
%   % FRF do sistema primário
% H=1./((-om.^2*m1+k1)+(1i*c1*om));
%   % FRFs do sistema composto
% for j=1:length(c2)
%   ab=-om.^2*m1+1i*om*(c1+c2(j))+(k1+k2);
%   bb=1i*om*c2(j)+k2;
%   cb=-om.^2*m2+1i*om*c2(j)+k2;
%   H1(j,:)=cb./(ab.*cb-bb.^2);
% end
%   % visualização
% figure(2)
% semilogx(om,20*log10(abs(H)),'r-',om,20*log10(abs(H1(1,:))),'b-',...
%          om,20*log10(abs(H1(2,:))),'b--',om,20*log10(abs(H1(3,:))),'c-',...
%          om,20*log10(abs(H1(4,:))),'c--')
% xlabel('Frequência (rad/s)');
% ylabel('Magnitude de FRF (dB ref. 1m/N)');
% axis([10 100 -120 -20]);
% grid on
% legend('SP (c1=1kg/s)','SC (c2=0.001kg/s)','SC (c2=1kg/s)','SC (c2=5kg/s)',...
%        'SC (c2=10kg/s)')

% % passo 5 - sintonização e amortecimento "ideais"
%   % dados adicionais
% c1=1; % coef. de amortecimento do sistema primário
% c2=[0.001 1 2 5]; % coefs. de amortecimento do neutralizador (0.001, 1, 2, 5)
%   % faixa de frequência
% om=logspace(1,2,1000);
%   % FRF do sistema primário
% H=1./((-om.^2*m1+k1)+(1i*c1*om));
%   % valores "ideais"
% mi=m2/m1; % razão de interesse entre massas
% m2id=mi*m1; % nova massa do neutralizador (kg)
% omid=om1*(1/(1+mi)); % frequência de sintonização "ideal" (rad/s)
% k2id=m2id*omid^2; % constante de rigidez "ideal" (N/m)
% c2id=(2*m2id*omid)*sqrt((3*mi)/(8*(1+mi))); % coef. de amort. "ideal" (kg/s)
%   % inclusão de amortecimento "ideal"
% c2=[c2 c2id];
%   % FRFs do sistema composto
% for j=1:length(c2)
%   ab=-om.^2*m1+1i*om*(c1+c2(j))+(k1+k2id);
%   bb=1i*om*c2(j)+k2id;
%   cb=-om.^2*m2id+1i*om*c2(j)+k2id;
%   H1(j,:)=cb./(ab.*cb-bb.^2);
% end
%   % visualização
% figure(2)
% semilogx(om,20*log10(abs(H)),'r-',om,20*log10(abs(H1(1,:))),'b-',...
%          om,20*log10(abs(H1(2,:))),'b--',om,20*log10(abs(H1(3,:))),'c-',...
%          om,20*log10(abs(H1(4,:))),'c--',om,20*log10(abs(H1(5,:))),'k-')
% xlabel('Frequência (rad/s)');
% ylabel('Magnitude de FRF (dB ref. 1m/N)');
% axis([10 100 -120 -20]);
% grid on
% legend('SP (c1=1kg/s)','SC (c2=0.001kg/s)','SC (c2=1kg/s)',...
%        'SC (c2=5kg/s)','SC (c2=10kg/s)','SC (c2 "ideal")')

% % passo 6 - novos parâmetros
%   % novas frequências naturais do sistema composto
%     % coeficientes da equação biquadrada
% c(1)=1;
% c(2)=0;
% c(3)=-((k1/m1)+(1+(m2id/m1))*(k2id/m2id));
% c(4)=0;
% c(5)=(k1*k2id)/(m1*m2id);
%   % raízes da equação biquadrada
% rs=roots(c); 
%   % raízes ordenadas
% rs=sort(rs); 
%   % frequências naturais
% omn=rs(3:4);
%   % visualização dos novos parâmetros
% disp(' ')
% disp(['A nova massa do neutralizador é ',num2str(m2id),' kg.'])
% disp(' ')
% disp(['A nova frequência de sintonização é ',num2str(omid),' rad/s.'])
% disp(' ')
% disp(['A nova constante de rigidez do neutralizador é ',num2str(k2id),' N/m.'])
% disp(' ')
% disp(['O coeficiente de amortecimento do neutralizador é ',num2str(c2id),...
%       ' kg/s'])
% disp(' ')
% disp(['As novas frequências naturais do sistema composto são ',num2str(omn(1)),...
%       ' rad/s e ',num2str(omn(2)),' rad/s.'])
% disp(' ')