% Unit 4 (MATLAB Course @ UAM)
% (c) Alexander Knebe

%============
% CANNONBALL
%============
clear all, clc

% constants
gEarth     = -9.79757; % m/s^2
gMoon      = -1.622;   % m/s^2
gy         = gEarth;
gx         = 0.0;

% input values
N     = 100;     % points
x0    = 3.;      % m
y0    = 2.;      % m
v0    = 10.;     % m/s (NOTE: total velocity!)
alpha = 23.;     % degrees

% inital velocity transformed to [vxy, vy0]
vx0   = v0*cos(alpha * 2*pi/360.);
vy0   = v0*sin(alpha * 2*pi/360.);

% time arrays (Note: Tend is the time until return to y=y0)
Tend  = -(2*vy0)/gy;
t     = linspace(0,Tend,N);

% analytical cannonball solution
x  = x0  + vx0.*t + 0.5.*gx.*t.^2;
y  = y0  + vy0.*t + 0.5.*gy.*t.^2;

%vx = vx0 .* ones(1,N); % making explicit use of the knowledge that gx=0
vx = vx0 +  gx.*t;
vy = vy0 +  gy.*t;

% energy conservation
v = sqrt(vx.^2 + vy.^2);
E = 0.5*v.^2-gy.*y;

% numerical velocities
dx     = diff(x);
dy     = diff(y);
dt     = diff(t);
vx_num = dx./dt;
vy_num = dy./dt;

% ADVANCED: numerical derivatives of numerical velocities
dvx  = diff(vx_num); % NOTE: dvx will have length N-2
dvy  = diff(vy_num);
il   = (1:1:N-2);
ir   = (2:1:N-1);
d2t  = t(ir)-t(il);
gx_num = dvx./d2t;
gy_num = dvy./d2t;

% save some of the data to file
save('cannonball','x','y','vx','vy','t')
