% Frequências naturais, modos de vibrar, vibrações amortecidas
% e receptâncias de ponte rolante modelada com 2 gdl


clc
clearvars
close all


% dados do sistema
ngdl=2; % número de graus de liberdade
m1=5000; % massa da carga
m2=1000; % massa do carrinho
k1=3e5; % constante de rigidez associada ao cabo
k2=3.09e6; % constante de rigidez associada à viga
zeta=[0.02; 0.001]; % razões de amortecimento

% matrizes de inércia e rigidez
M=[m1 0; 0 m2];
K=[k1 -k1; -k1 k1+k2];


% 1a. parte - Cálculo de frequências naturais e modos de vibrar

% passo 1 - decomposição de Choleski, transposição e inversas
L=chol(M);
Lt=L';
Li=inv(L);
Lti=inv(Lt);

% passo 2 - matriz de rigidez normalizada pela massa
Ktil=Lti*K*Li;

% passo 3 - problema de autovalores
[V,D]=eig(Ktil);
  % frequências naturais em rad/s
omn=sqrt(diag(D));
  % ordenamento das frequências naturais
[omn,iv]=sort(omn);
  % apresentação das frequências naturais
disp('As frequências naturais, em rad/s, são')
disp(' ')
for im=1:length(omn)
  disp(['ômega_',num2str(im),' = ',num2str(omn(im)),' rad/s'])
end
disp(' ')

% passo 4 - autovetores normalizados e ordenados
P=zeros(size(V));
for il=1:length(iv)
  P(:,il)=V(:,iv(il));
end

% passo 5 - matriz modal, transposições e inversa
%          (modos de vibrar em colunas da matriz modal)
T=Li*P; 
Tt=T';
Pt=P';
Ti=Pt*L;
  % visualização dos modos
figure(1)
plot(T,'o--')
vcf=1:1:size(T,2); % vetor de índices de coordenadas físicas
set(gca,'xtick',vcf); % restrição correspondente nos "ticks"
grid on
xlabel('índice de coordenada física')
ylabel('amplitude modal')
legend('modo 1','modo 2')
title('Modos de vibrar')


% % 2a. parte - Cálculo de vibrações amortecidas
% 
% % seleção do caso de interesse
% caso=1;
% 
% % passo 6 - matriz espectral
% Lb=Pt*Ktil*P;
% 
% % passo 7 - condições iniciais em coordenadas modais e  
% % passo 8 - resolução das equações pertinentes
%   % vetor de tempo
% t=linspace(0,50,1024); 
%   % casos considerados
% switch caso
%   case 0 % sem determinação de vibrações
% disp('Frequências naturais e modos de vibrar determinados.')
% disp(' ')
% return
% 
%   case 1 % vibrações livres amortecidas 
%   % condições iniciais
% x0=[0.0705; 0.0065];
% xp0=[0; 0];
%   % condições iniciais em coordenadas modais
% r0=Ti*x0;
% rp0=Ti*xp0;
%   % respostas em coordenadas modais (amort. subcrítico)
% ni=omn.*sqrt(1-zeta.^2);
% delta=zeta.*omn;
% C1=r0;
% C2=(r0.*delta+rp0)./ni;
% C=sqrt(C1.^2+C2.^2);
% fi=atan2(C1,C2);
% rcom=zeros(length(omn),length(t));
% for ir=1:length(omn)
%   rcom(ir,:)=C(ir)*exp(-delta(ir)*t).*sin(ni(ir)*t+fi(ir));
% end
%   % visualização das respostas em coordenadas modais
% figure(3)
% subplot(2,1,1)
% title('Vibrações modais livres amortecidas')
% plot(t,rcom(1,:),'r-')
% grid on
% xlabel('tempo (s)')
% ylabel('coordenada modal 1 (m)')
% axis([min(t) max(t) min(min(rcom')) max(max(rcom'))]) 
% subplot(2,1,2)
% plot(t,rcom(2,:),'r-')
% grid on
% xlabel('tempo (s)')
% ylabel('coordenada modal 2 (m)')
% axis([min(t) max(t) min(min(rcom')) max(max(rcom'))])
% 
%   case 2 % vibrações forçadas sob excitação harmônica
%   % condições iniciais
% x0=[0; 0];
% xp0=[0; 0];
%   % condições iniciais em coordenadas modais
% r0=Ti*x0;
% rp0=Ti*xp0;
%   % forças atuantes
% F=[0; 10]; % amplitudes
% om=3; % frequência
%   % vetor de amplitudes de forças modais
% Fm=Tt*F;
%   % forças modais harmônicas
% fm=Fm*sin(om*t);
%   % visualização das forças modais
% figure(2)
% plot(t,fm(1,:),'r-',t,fm(2,:),'b-')
% grid on
% xlabel('tempo (s)')
% ylabel('força modal (N)')
% legend('força 1','força 2')
% title('Forças modais')
%   % respostas em coordenadas modais
%     % respostas permanentes
% k=omn.^2;
% e=om./omn;
% H=(1./sqrt((1-e.^2).^2+(2*zeta.*e).^2));
% R=(Fm./k).*H;
% teta=atan2((2*zeta.*e),(1-e.^2));
% tetar=teta;
% rper=zeros(length(omn),length(t));
% for ir=1:length(omn)
%   rper(ir,:)=R(ir)*sin(om*t-tetar(ir));
% end
%     % respostas transientes (amort. subcrítico, ci's nulas)
% ni=omn.*sqrt(1-zeta.^2);
% delta=zeta.*omn;
% psi=atan2((2*zeta.*sqrt(1-zeta.^2)),(2*zeta.^2-(1-e.^2)));
% rtran=zeros(length(omn),length(t));
% for ir=1:length(omn)
%   rtran(ir,:)=R(ir)*(e(ir)*exp(-delta(ir)*t)/sqrt(1-zeta(ir)^2)).*...
%               sin(ni(ir)*t+psi(ir));
% end
%     % respostas completas
% rcom=rper+rtran;
%   % visualização das respostas em coordenadas modais
% figure(3)
% subplot(2,1,1)
% title('Vibrações modais sob excitação harmônica')
% plot(t,rper(1,:),'r-',t,rtran(1,:),'g-',t,rcom(1,:),'b-')
% grid on
% xlabel('tempo (s)')
% ylabel('coordenada modal 1 (m)')
% legend('per.','tran.','com.')
% axis([min(t) max(t) min(min(rcom')) max(max(rcom'))]) 
% subplot(2,1,2)
% plot(t,rper(2,:),'r-',t,rtran(2,:),'g-',t,rcom(2,:),'b-')
% grid on
% xlabel('tempo (s)')
% ylabel('coordenada modal 2 (m)')
% legend('rper','tran.','com.')
% axis([min(t) max(t) min(min(rcom')) max(max(rcom'))]) 
% 
%   case 3 % vibrações forçadas sob excitação impulsiva
%   % condições iniciais
% x0=[0; 0];
% xp0=[0; 0];
%   % condições iniciais em coordenadas modais
% r0=Ti*x0;
% rp0=Ti*xp0;
%   % forças atuantes
% I=[10; 0]; % intensidades
%   % vetor de intensidades de forças modais
% Im=Tt*I;
%   % respostas em coordenadas modais (amort. subcrítico)
% ni=omn.*sqrt(1-zeta.^2);
% delta=zeta.*omn;
% rcom=zeros(length(omn),length(t));
% for ir=1:length(omn)
%   rcom(ir,:)=(Im(ir)./ni(ir)).*exp(-delta(ir)*t).*sin(ni(ir)*t);
% end
%   % visualização das respostas em coordenadas modais
% figure(3)
% subplot(2,1,1)
% title('Vibrações modais sob excitação impulsiva')
% plot(t,rcom(1,:),'r-')
% grid on
% xlabel('tempo (s)')
% ylabel('coordenada modal 1 (m)')
% axis([min(t) max(t) min(min(rcom')) max(max(rcom'))]) 
% subplot(2,1,2)
% plot(t,rcom(2,:),'r-')
% grid on
% xlabel('tempo (s)')
% ylabel('coordenada modal 2 (m)')
% axis([min(t) max(t) min(min(rcom')) max(max(rcom'))]) 
% 
%   otherwise
% disp('Caso não especificado.')
% disp(' ')
% return
% end
% 
% % passo 9 - respostas globais em coordenadas físicas
% x=T*rcom;
%   % contribuições dos modos
% xm1=T(:,1)*rcom(1,:);
% xm2=T(:,2)*rcom(2,:);
%   % visualização das respostas em coordenadas físicas
% figure(4)
% subplot(2,1,1)
% title('Vibrações em coordenadas físicas')
% plot(t,xm1(1,:),'r-',t,xm2(1,:),'g-',t,x(1,:),'b-')
% grid on
% xlabel('tempo (s)')
% ylabel('vibração da massa 1 (m)')
% legend('modo 1','modo 2','global')
% axis([min(t) max(t) min(min(x')) max(max(x'))]) 
% subplot(2,1,2)
% plot(t,xm1(2,:),'r-',t,xm2(2,:),'g-',t,x(2,:),'b-')
% grid on
% xlabel('tempo (s)')
% ylabel('vibração da massa 2 (m)')
% legend('modo 1','modo 2','global')
% axis([min(t) max(t) min(min(x')) max(max(x'))]) 


% % 3a. parte - Representação global de receptâncias
% % (aplicável a ngdl = 2; comentar a 2a. parte acima)
% 
%   % vetor de frequências
% omx=linspace(0,70,2048);
%   % alocação de espaço
% H11=zeros(1,length(omx));
% H21=zeros(1,length(omx));
% H12=zeros(1,length(omx));
% H22=zeros(1,length(omx));  
%   % matriz de receptâncias
%     % montagem da matriz de amortecimento
% A(:,1)=[1 1]'; A(:,2)=omn.^2;
% b=2*omn.*zeta;
% z=A\b; % ver Eq. (7), Módulo 18
% af=z(1); % constante alfa
% bt=z(2); % constante beta
% C=af*M+bt*K; % matriz de amortecimento
%     % montagem da matriz de receptâncias
% for ix=1:length(omx)
%   H=inv(-omx(ix)^2*M+1i*omx(ix)*C+K); % ver Eq. (12), Módulo 18
%   H11(ix)=H(1,1);
%   H12(ix)=H(1,2);
%   H21(ix)=H(2,1);
%   H22(ix)=H(2,2);
% end
%   % visualização da matriz de receptâncias
%     % espectro de amplitude
% figure(5)
% subplot(2,2,1)
% plot(omx,20*log10(abs(H11)),'b-')
% grid on
% xlabel('Frequência (rad/s)')
% ylabel('H_{11} - Amplitude (dB re 1m/N)')
% subplot(2,2,2)
% plot(omx,20*log10(abs(H12)),'b-')
% grid on
% xlabel('Frequência (rad/s)')
% ylabel('H_{12} - Amplitude (dB re 1m/N)')
% subplot(2,2,3)
% plot(omx,20*log10(abs(H21)),'b-')
% grid on
% xlabel('Frequência (rad/s)')
% ylabel('H_{21} - Amplitude (dB re 1m/N)')
% subplot(2,2,4)
% plot(omx,20*log10(abs(H22)),'b-')
% grid on
% xlabel('Frequência (rad/s)')
% ylabel('H_{22} - Amplitude (dB re 1m/N)')
%     % espectro de fase
% figure(6)
% subplot(2,2,1)
% plot(omx,unwrap(atan2(imag(H11),real(H11))),'r-')
% grid on
% xlabel('Frequência (rad/s)')
% ylabel('H_{11} - Fase (rad)')
% axis([0 70 -8 2])
% subplot(2,2,2)
% plot(omx,unwrap(atan2(imag(H12),real(H12))),'r-')
% grid on
% xlabel('Frequência (rad/s)')
% ylabel('H_{12} - Fase (rad)')
% axis([0 70 -8 2])
% subplot(2,2,3)
% plot(omx,unwrap(atan2(imag(H21),real(H21))),'r-')
% grid on
% xlabel('Frequência (rad/s)')
% ylabel('H_{21} - Fase (rad)')
% axis([0 70 -8 2])
% subplot(2,2,4)
% plot(omx,unwrap(atan2(imag(H22),real(H22))),'r-')
% grid on
% xlabel('Frequência (rad/s)')
% ylabel('H_{22} - Fase (rad)')
% axis([0 70 -8 2])


% % % 4a. parte - Representação detalhada de receptâncias
% % % (aplicável a ngdl = 2; comentar 2a. e 3a. partes acima)
% 
%   % vetor de frequências
% omx=linspace(0,70,2048);
%   % matriz de receptâncias via parâmetros modais
% H11=zeros(1,length(omx));
% H21=zeros(1,length(omx));
% H12=zeros(1,length(omx));
% H22=zeros(1,length(omx));
% for ix=1:length(omx)
%   daux1=1/((omn(1)^2-omx(ix)^2)+1i*(2*zeta(1)*omn(1)*omx(ix)));
%   daux2=1/((omn(2)^2-omx(ix)^2)+1i*(2*zeta(2)*omn(2)*omx(ix)));
%   Daux=diag([daux1 daux2]);
%   Haux=T*Daux*Tt;
%   H11(ix)=Haux(1,1);
%   H21(ix)=Haux(2,1);
%   H12(ix)=Haux(1,2);
%   H22(ix)=Haux(2,2);
% end
%   % visualização da matriz de receptâncias
% figure(7)
% subplot(2,2,1)
% plot(omx,20*log10(abs(H11)),'b-')
% grid on
% xlabel('frequência (rad/s)')
% ylabel('Módulo de H_{11} (dB re 1m/N)')
% subplot(2,2,2)
% plot(omx,20*log10(abs(H12)),'b-')
% grid on
% xlabel('frequência (rad/s)')
% ylabel('Módulo de H_{12} (dB re 1m/N)')
% subplot(2,2,3)
% plot(omx,20*log10(abs(H21)),'b-')
% grid on
% xlabel('frequência (rad/s)')
% ylabel('Módulo de H_{21} (dB re 1m/N)')
% subplot(2,2,4)
% plot(omx,20*log10(abs(H22)),'b-')
% grid on
% xlabel('frequência (rad/s)')
% ylabel('Módulo de H_{22} (dB re 1m/N)')


% % 5a. parte - Representação de receptância selecionada
% % (aplicável a ngdl = 2; comentar 2a. e 3a. partes acima)
% 
%   % seleção de receptância
% es=1;
% er=2;
%   % montagem de receptância selecionada via parâmetros modais
% Hsrx=zeros(length(omn),length(omx));
% Hsr=zeros(1,length(omx));
% for ix=1:length(omx)
%   for im=1:length(omn)
%     Hsrx(im,ix)=(T(es,im)*T(er,im))/...
%                 ((omn(im)^2-omx(ix)^2)+1i*(2*zeta(im)*omn(im)*omx(ix)));
%   end
%   Hsr(ix)=sum(Hsrx(:,ix));
% end
%   % visualização da função receptância selecionada
% figure(8)
% ifrf=[num2str(es),num2str(er)];
%     % espectro de amplitude
% mfrf=['H_{',ifrf,'} - Amplitude (dB re 1m/N)'];
% tfrf=['Receptância H_{',ifrf,'}'];
% subplot(2,1,1)
% plot(omx,20*log10(abs(Hsrx(1,:))),'r-',...
%      omx,20*log10(abs(Hsrx(2,:))),'g-',...
%      omx,20*log10(abs(Hsr)),'b-')
% grid on
% xlabel('Frequência (rad/s)')
% ylabel(mfrf)
% title(tfrf)
% legend('modo 1','modo 2','total')
%     % espectro de fase
% afrf=['H_{',ifrf,'} - Fase (rad)'];
% subplot(2,1,2)
% plot(omx,atan2(imag(Hsrx(1,:)),real(Hsrx(1,:))),'r-',...
%      omx,atan2(imag(Hsrx(2,:)),real(Hsrx(2,:))),'g-',...
%      omx,atan2(imag(Hsr),real(Hsr)),'b-')
% grid on
% xlabel('Frequência (rad/s)')
% ylabel(afrf)
% legend('modo 1','modo 2','total')


% % 6a. parte - Estimação de parâmetros modais
% % (aplicável a ngdl = 2; comentar 2a. e 3a. partes acima)
% 
%   % frequências naturais e razões de amortecimento
%   % via receptância individual
%   % (rodar linhas 423 a 428, com linhas 429 a 451 comentadas,)
%   % (para leitura de valores; na sequência, rodar sem comentários nessas)
%   % (linhas, para inserção de valores e cálculo de parâmetros)
% figure(9) % uso de H21
% plot(omx,20*log10(abs(H21)),'b-') 
% grid on
% xlabel('Frequência (rad/s)')
% ylabel('H_{21} - Amplitude(dB re 1m/N)')
% title('Receptância H_{21}')
% for im=1:length(omn)
%   disp(' ')
%   omr=input([num2str(im),'a. frequência de ressonância? ']);
%   disp(' ')
%   om1=input(['Primeira frequência da ',num2str(im),'a. banda de -3 dB? ']);
%   disp(' ')
%   om2=input(['Segunda frequência da ',num2str(im),'a. banda de -3 dB? ']);
%   disp(' ')
%     % frequência natural estimada
%   omn2(im)=omr;
%     % razão de amortecimento estimada
%   zeta2(im)=(om2-om1)/(2*omr);
% end
% for im=1:length(omn)
%   % visualização dos parâmetros:
%   disp(' ')
%   disp(['A ',num2str(im),'a. frequência natural estimada é ',...
%         num2str(omn2(im)),' rad/s.'])
%   disp(' ')
%   disp(['A ',num2str(im),'a. razão de amortecimento estimada é ',...
%         num2str(zeta2(im)),'.'])
%   disp(' ')
% end