% Sistema com 1 gdl sob excitação periódica
% (força em puncionadeira)

clc
clearvars
close all

% dados do sistema
m=500;
k=1.25e6;
zeta=0.1;

% frequência natural (rad/s)
omn=sqrt(k/m);

% período e frequência da força
T=0.5;
om1=2*pi/T;

% visualização de resultados preliminares
disp(' ')
disp(['A frequência natural circular é ',num2str(omn),' rad/s.'])
disp(' ')
disp(['A frequência fundamental da excitação é ',num2str(om1),' rad/s.'])
disp(' ')

% % 1o. passo
% % série de Fourier da força
% % (atribuir a nth, na linha 55, o número 
% %  de termos desejado, tal que nto >= 1;
% %  a princípio, nth = 20.)
% 
% % vetor de tempo
% tf=1.0;
% np=500;
% t=linspace(0,tf,np);
% % função original
% for it=1:length(t)
%   if t(it) <= 0.2
%     ft(it)=5000;
%   elseif t(it) < 0.5
%     ft(it)=0;
%   elseif t(it) <= 0.7
%     ft(it)=5000;
%   else
%     ft(it)=0;
%   end
% end
% % série de Fourier
%   % termo constante
% A0=2000;
% F0=A0;
% f0=F0*ones(size(t));
%   % termos harmônicos
% nth=10; % número de termos harmônicos (oscilatórios)
% if nth > 0
%   for j=1:nth
%     A(j)=(5000/(j*pi))*sin(j*0.8*pi);
%     B(j)=(5000/(j*pi))*(1-cos(j*0.8*pi));
%     F(j)=sqrt(A(j)^2+B(j)^2);
%     if F(j) < 1e-10
%       tetaf(j)=0;
%     else
%       tetaf(j)=atan2(A(j),B(j));
%     end
%     f(j,:)=F(j)*sin(j*om1*t+tetaf(j));
%   end
% end
%   % visualização da força e de sua série de Fourier
% figure(1)
% plot(t,ft,'g.') % função original
% axis([0 1 -1500 6500])
% xlabel('t (s)')
% ylabel('f(t) (N)')
% grid on
% pause(4)
% hold on
% plot(t,f0,'r-') % termo constante
% legend('original','termo constante')
% hold off
% pause(4)
% fs=f0;
% for ip=1:nth
%   fs=fs+f(ip,:);
%   plot(t,ft,'g.',t,fs,'r-') % função original e termos da série
%   axis([0 1 -1500 6500])
%   xlabel('t (s)')
%   ylabel('f(t) (N)')
%   grid on
%   th=['harmônicos: ',num2str(ip)];
%   legend('original',th)
%   pause(3) % para nth > 20, reduzir para 2, ou 1 s
% end

% % 2o. passo
% % resposta permanente à série de Fourier da força
% 
% % montagem da resposta
%   % termo constante
% X0=F0/k;
% xp0=X0*ones(size(t));
%   % termos harmônicos
% for j=1:nth
%   om(j)=j*om1;
%   e(j)=om(j)/omn;
%   H(j)=1/(sqrt((1-e(j)^2)^2+(2*zeta*e(j))^2));
%   X(j)=(F(j)/k)*H(j);
%   if X(j) < 1e-10
%     tetae(j)=0;
%   else
%     tetae(j)=atan2((2*zeta*e(j)),(1-e(j)^2));
%   end
%   tetax(j)=tetae(j)-tetaf(j); 
%   xp(j,:)=X(j)*sin(j*om1*t-tetax(j));
% end
% % visualização da força e da resposta permanente
% figure(2)
% plot(t,f0/1000,'r-',t,xp0*1000,'b-') % termos constantes
% axis([0 1 -4.5 8.5]) 
% xlabel('t (s)')
% ylabel('f(t) (kN), xp(t) (mm)')
% grid on
% legend('f(t) (cte.)','xp(t) (cte.)')
% pause(4)
% fs=f0;
% xps=xp0;
% for ip=1:nth
%   fs=fs+f(ip,:);
%   xps=xps+xp(ip,:);
%   plot(t,fs/1000,'r-',t,xps*1000,'b-') % constantes + harmônicos
%   axis([0 1 -4.5 8.5]) 
%   xlabel('t (s)')
%   ylabel('f(t) (kN), xp(t) (mm)')
%   grid on
%   hf=['f(t) (har: ',num2str(ip),')'];
%   hxp=['xp(t) (har: ',num2str(ip),')'];
%   legend(hf,hxp)
%   pause(3)
% end

% % 3o. passo
% % espectros de amplitude e fase da força
% % (atribuir nth = 10 na linha 55)
% 
% figure(3)
% subplot(3,2,1)
% ome=[0 om]; % vetor de frequências
% Fe=[F0 F]; % vetor de amplitudes de força
% stem(ome,Fe,'r.')
% axis([0 max(ome)+10 0 5000])
% xlabel('\omega (rad/s)');
% ylabel('amplitude f(t) (N)');
% grid on
% subplot(3,2,2)
% stem(om,tetaf,'r.')
% axis([0 max(ome)+10 -1.5*pi +pi])
% xlabel('\omega (rad/s)');
% ylabel('fase f(t) (rad)');
% grid on
 
% % 4o. passo
% % espectros de amplitude e fase da resposta permanente
% % (comentar linhas 69 a 93 e 116 a 139)
% 
% subplot(3,2,5)
% Xe=[X0 X]; % vetor de amplitudes de vibração
% stem(ome,Xe*1000,'b.')
% axis([0 max(ome)+10 0 5])
% xlabel('\omega (rad/s)');
% ylabel('amplitude xp(t) (10^{-3}m)');
% grid on
% subplot(3,2,6)
% stem(om,tetax,'b.')
% axis([0 max(ome)+10 -1.5*pi +pi])
% xlabel('\omega (rad/s)');
% ylabel('fase xp(t) (rad)');
% grid on
 
% % 5o. passo
% % gráficos (espectros) do fator de amplitude e do fator de fase
% 
% subplot(3,2,3)
% omv=linspace(0,max(ome),100);
% % fator de amplitude (deslocamento)
% Hx=(1./sqrt((1-(omv/omn).^2).^2+(2*zeta*(omv/omn)).^2));
% plot(omv,Hx,'m-')
% xlim([0 max(ome)+10])
% xlabel('\omega (rad/s)');
% ylabel('fator de amplitude');
% grid on
% subplot(3,2,4)
% % fator de fase correspondente
% tetav=atan2((2*zeta*(omv/omn)),(1-(omv/omn).^2));
% plot(omv,tetav,'m-')
% axis([0 max(ome)+10 -1.5*pi +pi])
% xlabel('\omega (rad/s)');
% ylabel('fator de fase (rad)');
% grid on