%% projectile motion using textbook solution
v0 = 30;
theta = 60;
g = 9.81;
t = linspace(0, 5);
vx0 = v0 * cosd(theta);
vy0 = v0 * sind(theta);
xtrue = vx0 * t;
ytrue = vy0 * t - 0.5 * g * t .* t;
plot(xtrue, ytrue, 'r-');
xlabel('distance (meters)');
ylabel('height (meters)');
hold on

%% projectile motion using Euler's method
x0 = 0;
y0 = 0;
dt = 1;
[x, y, vx, vy, t] = solveproj(5, dt, x0, y0, vx0, vy0);
plot(x, y, 'k-s')

dt = 0.5;
[x, y, vx, vy, t] = solveproj(5, dt, x0, y0, vx0, vy0);
plot(x, y, 'k-*')

dt = 0.25;
[x, y, vx, vy, t] = solveproj(5, dt, x0, y0, vx0, vy0);
plot(x, y, 'k-+')

legend('textbook solution', ...
    'dt = 1', 'dt = 0.5', 'dt = 0.25', ...
    'Location', 'south');
title('projectile motion solution using Euler''s method')

figure
plot(xtrue, ytrue, 'r-');
xlabel('distance (meters)');
ylabel('height (meters)');
hold on

dt = 0.01;
[x, y, vx, vy, t] = solveproj(5, dt, x0, y0, vx0, vy0);
plot(x(1:10:end), y(1:10:end), 'k.')

legend('textbook solution', ...
    'dt = 0.01', ...
    'Location', 'south');
title('projectile motion solution using Euler''s method')

%% projectile motion with drag
figure
plot(xtrue, ytrue, 'r-');
xlabel('distance (meters)');
ylabel('height (meters)');
hold on
[x, y, vx, vy, t] = solveproj2(5, dt, x0, y0, vx0, vy0);
plot(x(1:10:end), y(1:10:end), 'k.')

legend('textbook solution (no drag)', ...
    'dt = 0.01', ...
    'Location', 'south');
title('projectile motion with drag using Euler''s method')

