% -- Sanity check trajectories computed with and without drag --

theta = pi/3;
v0 = s*[cos(theta); sin(theta)];

% Compute the reference trajectory (absent air resistance)
% x(y) = s*cos(theta)*t
% y(t) = -g*t^2/2 + s*sin(theta)*t
% up to time tfinal = 2*s*sin(theta)/g
%
tfinal = 2*v0(2)/g;
tref = linspace(0,tfinal)';
xref = v0(1)*tref;
yref = (v0(2)-g/2*tref).*tref;

% Compute the same reference trajectory with ode45
y0 = [0; 0; v0];
refopt = opt;
refopt.c = 0;
[tout,yout] = ode45(@(t,y) fball(t,y,refopt), tref, y0);

% Compute a similar trajectory with air drag on (no wind)
dopt = opt;
dopt.w = 0;
[toutd,youtd] = ode45(@(t,y) fball(t,y,dopt), tref, y0);

% Do a comparison between analytical and numerical solutions (no drag)
fprintf('Max x error: %g\n', norm( xref-yout(:,1), inf));
fprintf('Max y error: %g\n', norm( yref-yout(:,2), inf));

% Visually compare solutions for drag and no drag
plot(xref, yref, 'r:', youtd(:,1), youtd(:,2), 'b-');
legend('No drag', 'Drag');