% Unit 1 (MATLAB Course @ UAM)
% Harmonic Oscillator exercise
% (c) Alexander Knebe

%======================================================
% the harmonic oscillator
%======================================================
clear all
close all
clc


disp('======================================')
disp('      the harmonic oscillator')
disp('======================================')

%------------------------------------
% the physical problem
%------------------------------------

% initial conditions (to be varied)
t0 = 0
x0 = 1.2
v0 = 3.7
k  = 1.5
m  = 1.4

w0  = sqrt(k/m);
A   = sqrt(x0.^2+(v0/w0).^2);
phi = atan(w0.*x0./v0);

%------------------------------------
% time vector
%------------------------------------
N  = 20;
tN = 2*pi;
t  = linspace(0,tN,N);


%------------------------------------
% plot velocity
%------------------------------------
figure(1)

% analytically
v = A.*w0.*cos(w0.*t+phi);
plot(t,v)
xlabel('t')
ylabel('v')

%------------------------------------
% plot position into same figure
%------------------------------------
figure(1), hold on

% numerically
x_num = x0 + cumsum((v(1:N-1)+v(2:N))./2 .* diff(t));
% Note: we are not plotting at the mid-points as cumsum() gives the
%       cumulative sum up to the right end of each interval; hence we
%       also use that right end of the interval for plotting
plot(t(2:N),x_num,'o')

% analytically
x = A.*sin(w0.*t+phi);
plot(t,x,'r')

%------------------------------------
% plot acceleration into same figure
%------------------------------------
figure(1), hold on

% numerically
a_num = diff(v)./diff(t);
plot(t(1:N-1)+diff(t)./2,a_num,'x')

% analytically
a = -A.*w0.^2.*sin(w0.*t+phi);
plot(t,a,'y')

legend('velocity','numerical position','analytical position','numerical acceleration','analytical acceleration')