% Exercise 4: Response at periodic excitation
% Excitation data
%
A=1.0;                  % Amplitude
T0=0.1;                 % Breadhth of the signal
m=2;                    % Integer parameter tha define the period
T1=m*T0;                % Period of the periodic signal
NS=400;   	            % Numbers of terms of the series
%
N=2000;                 % Number of points in the time domain
detaT=0.00131;          % Increment in the time domain
T=detaT*N;              % Final time
f1=1/T1;                % Fundamental frequency, [Hz], of the excitation signal
t=0:detaT:T;            % time vector [s]
omeg1=2*pi*f1;          % Fundamental frequency, [rad/s], of the excitation signal
%
% Sistem data
%    
massa=20;               % mass
ksi=0.05;             % damping coeficiente
k=(2.*omeg1)^2*massa;  % stiffness parameters
%k=k+3690;
omegan=sqrt(k/massa);	% Non damped natural frequency [rad/s]
C=2*massa*ksi*omegan; 	% Damping coeficient
fn=omegan/(2*pi)        % Non damped natural frequency [Hz]
F0=A/m;                 % Half value of f(t)
X0=(1/k)*F0;            % Half value of x(t)

for it=1:length(t)
   soma=0.0;
   soma2=0.0;
  for ii=1:NS
    omeg(ii)=ii*omeg1;
%    
    F(ii)=A/m*sin(omeg(ii)*T0/2)/(omeg(ii)*T0/2);
    H(ii)=1/(-massa*omeg(ii)^2+1i*omeg(ii)*C+k); 
    X(ii)=H(ii)*F(ii);
%    
    soma=soma+conj(F(ii))*exp(-1i*omeg(ii)*t(it))+F(ii)*exp(1i*omeg(ii)*t(it));
    soma2=soma2+conj(X(ii))*exp(-1i*omeg(ii)*t(it))+X(ii)*exp(1i*omeg(ii)*t(it));
  end
  f(it)=soma+F0;
  x(it)=soma2+X0;
end

figure(1)
subplot(2,1,1)
plot(t,f)
ylabel('Excitation [EU]')

subplot(2,1,2)
plot(t,x)
xlabel('Time [s]')
ylabel('Response [m]')

%
%
HH=1./(-massa*omeg.^2+1i*C*omeg+k);   % FRF - receptance 
%
figure(3)
subplot(3,1,1)
plot(omeg(1:50)/2/pi,abs(HH(1:50)))
ylabel('|Receptance| [m/N]')

subplot(3,1,2)
stem(omeg(1:50)/2/pi,abs(F(1:50)))
ylabel('Spectrum [N]')

subplot(3,1,3)
stem(omeg(1:50)/2/pi,abs(X(1:50)))
ylabel('Spectrum [m]')
xlabel('FRequency [Hz]')

