%XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX
%X             Exercise 1: Harmonic excitation models                      X
%XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX

close all
clear all
clc

scrsz = get(0,'ScreenSize');
figure('OuterPosition',[scrsz(3)/10 scrsz(4)/6 scrsz(3)/3 2*scrsz(4)/3])

% The excitation is the hamonic function:
%        1- f(t) = A cos(ome*t)+ B sin(ome*t)
%        2- f(t) = P cos(ome*t + fi)
%        3- f(t)=F(ome)*exp(i*ome*t)+conjug(F(ome))*exp(i*ome*t)
%               then F(ome)=(A-i*B)/2

f=10;                         % excitation frequency [Hz]
%
A=10;                        % Value of A 
B=20;                        % value of B
P=sqrt(A^2+B^2);             % P amplitude
fi=atan(-B/A);               % phase angle
FFF=(A-1i*B)/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
omeg=2*pi*f;                 % Excitation frequency in [rad/s]

%--------------------------------------------------------------------------
%                     Graphs using the 3 Models
%--------------------------------------------------------------------------
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,'c-')
title('Model 1')
axis([0 5/f -2*P 2*P])

subplot(3,1,2)
plot(t,y1,'b-')
ylabel('Harmonic function [N]')
title('Model 2')
axis([0 5/f -2*P 2*P])

subplot(3,1,3)
plot(t,y2,'g-')
xlabel('Time [s]')
title('Model 3')
axis([0 5/f -2*P 2*P])
%--------------------------------------------------------------------------
%                     Spectrum[N] and Phase angle
%--------------------------------------------------------------------------

Dome=0.013;                 %frequency increment
ome=0:Dome:3*omeg;          %frequency vector


FF=zeros(length(ome),1);         % creating the vector of the spectrum
for ii=1:length(ome)
    if abs(omeg-ome(ii))<=Dome 
           FF(ii)=FFF;
   end    
end

figure('OuterPosition',[scrsz(3)/2 scrsz(4)/6 scrsz(3)/3 2*scrsz(4)/3])

figure(2)
subplot(2,1,1)
stem(ome/2/pi,2*abs(FF))
xlabel('frequency [Hz]')
ylabel('Spectrum [N]')
axis([0 3*f 0 1.5*P])

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

