clear all
close all
clc

% initial conditions
%----------------------------
gx    = 0.0;
gy    = -9.81;   % kg m/s^2

b     = 0.05;     % friction

t0    = 0.;      % s
x0    = 3.;      % m
y0    = 2.;      % m
v0    = 10.;     % m/s (NOTE: total velocity!)
alpha = 23.;     % degrees

vx0   = v0.*cos(alpha.*pi/180.);
vy0   = v0.*sin(alpha.*pi/180.);

% generate 2D vectors
g0 = [gx;  gy];
v0 = [vx0; vy0];
r0 = [x0;  y0];

%=============
% EXERCISE #2: get numerical and analytical solutions
%=============
g  = @(r,v,t) ( v );
h  = @(r,v,t) ( -b.*v.^2 .* v./norm(v) + g0 );

t0 = 0;
tN = -(2*vy0)/gy;
N  = 1000;

[rnum,vnum,t] = rk22_2D(g,h,r0,v0,t0,tN,N);

% get analytical cannonball solution, too (without friction)
rana = r0 + v0.*t+0.5.*g0.*t.^2;

% show them both
figure(1)
hold on
plot(rnum(1,:),rnum(2,:),'.','LineWidth',20)
plot(rana(1,:),rana(2,:))

%--------------------------------------------------------------------------
% EXERCISE #3: find peak position of the analytical solution
%--------------------------------------------------------------------------
vy       = @(t)(vy0+gy.*t);
[tmax,n] = bisec(vy, 0, tN);
rpeak    = r0 + v0.*tmax+0.5.*g0.*tmax.^2;

plot(rpeak(1),rpeak(2),'gx','MarkerSize',30)

%--------------------------------------------------------------------------
% EXERCISE #4: find peak position of the numerical solution
%--------------------------------------------------------------------------
dy    = diff(rnum(2,:));
roots = veczeros(dy);
plot(rnum(1,roots(1)),rnum(2,roots(1)),'cx','MarkerSize',30)


legend('numerical solution (w/ friction)','analytical solution (w/o friction)','analytical peak position','numerical peak position')
xlabel('x [m]')
ylabel('y [m]')
xlim([3,12])
