% Sistema com 1 gdl em vibração livre

clc;
clearvars;
close all;

% dados do sistema
m1=1000; % kg
m2=200; % kg
k=1e7; % N/m
c=4e5; % kg/s

% dados gerais
g=9.81; % m/s^2
h=1; % m

% passo 1 - condições de impacto
 
% % velocidade da bigorna (1)
% % antes do impacto
% v1a=0;
% % velocidade do martelo (2)
% % antes do impacto
% v2a=sqrt(2*g*h);
% % coeficiente de restituição
% e=1;
% % formulação matricial do impacto
%   % conservação da quantidade de movimento
%     % m1.v1d + m2.v2d = m1.v1a + m2.v2a
%   % coeficiente de restituição
%     % e = (v2d - v1d)/(v1a - v2a)
% A=[m1 m2; -1 1];
% b=[m1*v1a+m2*v2a;e*(v1a-v2a)];
% % velocidades da bigorna (1)
% % e do martelo (2) após o impacto
% v=A\b; % resolução de sistema linear
% disp(['As velocidades da bigorna e do ',...
%       'martelo, em m/s, são, respectivamente, '])
% disp(' ')
% disp(['v(1) = ',num2str(v(1)),'; v(2) = ',...
%       num2str(v(2)),'.']);
% disp(' ')

% passo 2 - parâmetros e condições iniciais
 
% % parâmetros do sistema:
% % massa vibrante
% m=m1;
% % 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; 
% % condições iniciais:
% x0=0;
% v0=v(1); % velocidade da bigorna após impacto
% % visualização de parâmetros
% % e condições iniciais:
% 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 constante de amortecimento ',...
%       'crítico é igual a ',num2str(cc),' kg/s.'])
% disp(' ')
% disp(['A razão de amortecimento é igual a ',...
%        num2str(zeta),'.'])
% disp(' ')
% disp(['O deslocamento inicial é igual a ',...
%        num2str(x0),' m.'])
% disp(' ')
% disp(['A velocidade inicial é igual a ',...
%        num2str(v0),' m/s.'])
% disp(' ')

% passo 3 - cálculo do deslocamento
%           (vibração livre)
 
% % geração de vetor de tempo
% t=linspace(0,0.6,1000);
% % parâmetro delta
% delta=zeta*omn;
%   % amortecimento subcrítico
% if zeta < 1
% % parâmetro ni
% ni=omn*sqrt(1-zeta^2);
% % constante C e ângulo fi
% A=x0;
% B=(x0*delta+v0)/ni;
% C=sqrt(A^2+B^2);
% fi=atan2(A,B);
% % vibração correspondente
% x=C*exp(-delta*t).*sin(ni*t+fi);
%   % amortecimento crítico
% elseif zeta == 1
% % constantes auxiliares  
% A=x0;
% B=x0*omn+v0;
% % vibração correspondente
% x=A*exp(-omn*t)+B*t.*exp(-omn*t);
% % % formulação alternativa
% % % vibração correspondente
% % x=(A+B*t).*exp(-omn*t);
%   % amortecimento supercrítico
% else
% % parâmetro gama
% gama=omn*sqrt(zeta^2-1);
% % constantes auxiliares
% A=x0/gama;
% B=v0/gama;
% % vibração correspondente
% x=A*exp(-delta*t).*...
%   (delta*sinh(gama*t)+gama*cosh(gama*t))+...
%   B*exp(-delta*t).*sinh(gama*t);
% % % formulação alternativa
% % % fatores s
% % s1=-delta+gama;
% % s2=-delta-gama;
% % % constantes auxiliares
% % A=(v0-s2*x0)/(s1-s2);
% % B=(s1*x0-v0)/(s1-s2);
% % % vibração correspondente
% % x=exp(-delta*t).*...
% %   (A*exp(gama*t)+B*exp(-gama*t)); 
% end
%  % visualização do gráfico:
% figure(1);
% plot(t,x,'r-')
% axis([0 0.6 -0.015 0.015])
% hold on
% plot([0 0.6],[0 0],'k-')
% hold off
% xlabel('t (s)')
% ylabel('x(t) (m)')
% grid on

% passo 4 - alteração de amortecimento
%           (c em 0, 2e4 e 6e4 na linha 11)

% passo 5 - efeito da razão de amortecimento
% 1) atribuir c = 2e4 na linha 11;
% 2) retirar comentários da linha 152 até a linha 162;
% 3) rodar;
% 4) atribuir c = 6e4 na linha 11;
% 5) rodar de novo;
% 6) ao final, voltar a comentar o passo 5.

% figure(2);
% xe=C*exp(-delta*t);
% xs=C*sin(ni*t+fi);
% plot(t,x,'r-',t,xs,'m:',t,xe,'m--',t,-xe,'m--')
% axis([0 0.6 -0.015 0.015])
% hold on
% plot([0 0.6],[0 0],'k-')
% hold off
% xlabel('t (s)')
% ylabel('x(t) (m)')
% grid on

% passo 6 - alteração de amortecimento
%           (c em 2e5 e 4e5 na linha 11)

% passo 7 - distinção entre os casos
%           crítico e supercrítico
% 1) comentar clearvars na linha 4;
% 2) atribuir c = 2e5 na linha 11;
% 3) retirar comentário da linha 178 e rodar;
% 4) comentar linha 178;
% 5) atribuir c = 4e5 na linha 11;
% 6) retirar comentários da linha 179 até o final;
% 7) rodar de novo.

% xc=x; % caso crítico
% xsc=x; % caso supercrítico
% % visualização de casos
% pause(4)
% figure(3);
% plot(t,xc,'r-')
% axis([0 0.3 -0.01 0.01])
% xlabel('t (s)')
% ylabel('x(t) (m)')
% grid on
% hold on
% plot([0 0.6],[0 0],'k-')
% pause(3)
% plot(t,xsc,'b-')
% hold off
% legend('crítico','supercrítico')