% Exercise 3: Response of the one degree of freedom system when 
%             it is hamonically solicited
clear all
clc
%
% System data
%
massa=10;                    % equivalent mass of the system
ksi=0.005;                  % damping ratio    
k=1e4;                      % equivalent stiffness
omegan=sqrt(k/massa);		% no damped natual frequency [rad/s]
c=2*massa*ksi*omegan;       % equivalent viscous damping
fn=omegan/(2*pi)            % no damped natural frequency [Hz]
%
% Excitation data
A=0;
B=1;
%
% Time domain data
N=3000;                     % number of points
detaT=0.00431;               % increment time
T=detaT*N;                  % final time
%
t=0:detaT:T;                % time vector
omega=30*2*pi;           % excitation frequency in [rad/s]
%
F=(A-1i*B)/2;               % Complex amplitude
H=1./(-massa*omega.^2+1i*c*omega+k);    % FRF - receptance 
%                                        
X=H*F;                                  % Response for diferent excitation frequencies 
%                                         
for it=1:length(t)
    f1(it)=F*exp(1i*omega*t(it))+conj(F)*exp(-1i*omega*t(it));
    
    x1(it)=X*exp(1i*omega*t(it))+conj(X)*exp(-1i*omega*t(it));
end
%
% Results
%
subplot(2,1,1)
plot(t,f1)
%
ylabel('Excitation [EU]')
subplot(2,1,2)
plot(t,x1)
ylabel('Response [EU]')
%
%
% Response in the frequency domain for a given excitation frequency
%
Dome=1/max(t)*2*pi;
omega_max=1/(2*(t(2)-t(1)))*2*pi;
ome=0:Dome:omega_max-Dome;          % vector frequency (variable)
%
HH=1./(-massa*ome.^2+1i*c*ome+k);   % FRF - receptance 
figure(3)
plot(ome/2/pi,abs(HH(1:length(ome))))
%
FF=fft(f1)/length(t);

ff1=ifft(FF)*length(t);
figure(7)
plot(t,ff1)

XF=fft(x1)/length(t);
xx1=ifft(XF)*length(t);
figure(8)
plot(t,xx1)


figure(4)
plot(ome/2/pi,abs(FF(1:length(ome))))
XX = FF(1:length(ome)).*HH(1:length(ome));

figure(5)
plot(ome(1:length(ome))/2/pi,abs(XX))

XXX=[XX(1:length(XX)-1) 0 conj(fliplr(XX(2:length(XX)-1)))];

xx=ifft(XXX)*length(XXX);
%novo t
DTnovo=max(t)/length(xx);
tnovo=0:DTnovo:max(t)

figure(6)
plot(tnovo(1:length(tnovo)-1),xx)

% figure(2)
% subplot(3,1,1)
% plot(ome,abs(HH))
% ylabel('|Receptance| [m/N]')
% %
% omegaext=omega(3);                  % Chose one excitation frequency 
% %                                    
% FFF=(A-1i*B)/2;                     % Complex amplitude of the excitation
% %
% HHH=1./(-massa*omegaext.^2+1i*c*omegaext+k); % Value of the FRF for the 
% %                                              excitation frequency    
% XXX=HHH*FFF;                        % Response
% %
% FF=zeros(1,length(ome));            % bulding a vector force 
% XX=zeros(1,length(ome));            % bulding a vector displacement
% %
% for ii=1:length(ome)
%     if abs(omegaext-ome(ii))<=Dome 
%            FF(ii)=FFF;
%            XX(ii)=XXX; 
%    end    
% end
% 
% 
% subplot(3,1,2)
% plot(ome,abs(FF))
% ylabel('|Excitation| [N]')
% 
% subplot(3,1,3)
% plot(ome,abs(XX))
% ylabel('|Response| [m]')
% xlabel('frequency [rad/s]')
% 
% 
