% yp = fball(t,y,opt)
% Compute the right hand side of an ODE for projectile motion
% in the presence of wind and air drag. The opt structure should
% describe a scaled air drag coefficient (c), the gravitational field (g),
% and the horizontal wind speed (w).
%
function yp = fball(t,y,opt)

g = opt.g; % Gravitational field
c = opt.c; % Scaled drag coefficient
w = opt.w; % Horizontal wind speed

% Unpack position and velocity
x = y(1:2);
v = y(3:4);

% Velocity relative to wind
vv = v;
vv(1) = vv(1)-w;

% Compute acceleration
a = -c*norm(vv)*vv;
a(2) = a(2)-g;

% Return
yp = [v; a];