% 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
%
f=3;                         % excitation frequency
%
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]
%
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*omeg;

FF=zeros(length(ome),1);
for ii=1:length(ome)
    if abs(omeg-ome(ii))<=Dome 
           FF(ii)=FFF;
   end    
end


figure(2)
subplot(2,1,1)
plot(ome/2/pi,2*abs(FF))
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])