% Rotina Matlab para determinação da resposta de um 
% sistema de 01 GDL com amortecimento viscoso
% sujeito à várias excitações transientes
% com a hipótese de condições iniciais nulas
%
%   criado em 18/09/2017 - SEM504 - Dinâmica Estrutural
%
clear all
%   parametros do sistema
m = 2 ;   % massa do sistema (kg)
fn = 8; % frequencia natual (hz)
k = (2*pi*fn)^2*m;   % constante elástica (N/m)
cc = 2*sqrt(k*m);   % amortecimento critico
c1 = 0.1*10^(-2)*cc;  % primeira constante de amortecimento
c2 = 5*10^(-2)*cc;  % segunda constante de amortecimento
%
%    definicao da primeira entrada
%
%    pulso retangular
p0 = 0.05*k ;     %amplitude do pulso 
%
t = linspace(0,3,3000)';     % vetor de tempos
%
t0 = 0.5;          % define o instante de inicio do pulso
t1 = 1.5 ;      % define a duracao do pulso a partir da origem
p1 = p0*heaviside(t-t0) - p0*heaviside(t-t1); 
%plot(t,p1)
%
%   Definicao da FT do sistema
%
num = [1];
den1 = [m c1 k] ;
den2 = [m c2 k] ;
H1 = tf(num,den1) ;
H2 = tf(num,den2) ;
%
%    Solução para a primeira entrada
%
u1 = lsim(H1,p1,t);
u2 = lsim(H2,p1,t);
%
%plot(t,p1/k,t,u2)
%
% definir a segunda forma de excitacao
%
w0 = 0.5*2*pi*fn;
p2 = p0*sin(w0*t).*(heaviside(t) - heaviside(t-pi/w0));
%
u3 = lsim(H2,p2,t);
%
plot(t,p2/k,t,u3)
%