% Exercise 1: Hamonic excitation models
%
clear all
clc
%
% The excitation is the hamonic function  f(t) = A cos(ome*t)+ B sin(ome*t)
%                                         f(t) = P cos(ome*t + fi)
%                                         f(t)=F(ome)*exp(i*ome*t)+conjug(F
%                                         (ome))*exp(i*ome*t)
% then F(ome)=(A-i*B)/2
%
f1=3.5;                         % excitation frequency
f2=5;
%
A1=10;                        % Value of A 
B1=0;                         % value of B
A2=0;
B2=20;

P1=sqrt(A1^2+B1^2);             % P amplitude
fi1=atan(-B1/A1);               % phase angle
FFF1=(A1-1i*B1)/2;              % complex amplitude

P2=sqrt(A2^2+B2^2);             % P amplitude
fi2=atan(-B2/A2);               % phase angle
FFF2=(A2-1i*B2)/2;              % complex amplitude
%
N=2000;                      % Number of point in the time domain
detaT=0.001;                 % time increment
T=detaT*N;                   % final time
t=0:detaT:T;                 % time vector

omeg1=2*pi*f1;                 % Excitation frequency in [rad/s]
omeg2=2*pi*f2;

y=A1*cos(omeg1*t)+B1*sin(omeg1*t)+A2*cos(omeg2*t)+B2*sin(omeg2*t);

plot(t,y)
%
% y=P*cos(omeg*t+fi);             % signal with model 1
% y1=A*cos(omeg*t)+B*sin(omeg*t); % signal with model 2
% y2= FFF*exp(1i*omeg*t)+ conj(FFF)*exp(-1i*omeg*t); % signal with model 3
% %
% subplot(3,1,1)
% plot(t,y)

% subplot(3,1,2)
% plot(t,y1)
% ylabel('Harmonic function [N]')
% 
% subplot(3,1,3)
% plot(t,y2)
% xlabel('Time [s]')



Dome=0.013;
ome=0:Dome:3*omeg2;

FF1=zeros(length(ome),1);
FF2=zeros(length(ome),1);
for ii=1:length(ome)
    if abs(omeg1-ome(ii))<=Dome 
           FF1(ii)=FFF1;
    end    
    if abs(omeg2-ome(ii))<=Dome 
           FF2(ii)=FFF2;
    end   
end


figure(2)
%subplot(2,1,1)
plot(ome/2/pi,2*abs(FF1),ome/2/pi,2*abs(FF2))
xlabel('frequency [Hz]')
ylabel('Spectrum [N]')
%axis([0 3*omeg 0 50])

% subplot(2,1,2)
% plot(ome/2/pi,angle(FF))
% xlabel('frequency [Hz]')
% ylabel('Phase angle [rad]')
% axis([0 3*omeg -pi 0])