% Projectile with air drag in SI units
clear all
% get input values
%m=input('Enter mass of particle');
m = 10;
g=9.81;

C=input('Enter drag coefficient:');

%%Initial position
%xold=input('Enter initial xposition');
%yold=input('Enter initial yposition');
xold =0;	%initial position at origin
yold =0;

%%Initial velocity
%vxold=input('Enter initial xvelocity');
%vyold=input('Enter initial yvelocity');
vxold = 40;
vyold = 40;

%%timestep
timestep=input('Enter time step');

%
%create array of projectile positions
i=1;
x(i)=xold;
y(i)=yold;
%
% Compute new positions
ynew = yold+timestep*vyold;
while ynew >= 0 % as long as it does not return to ground
	xnew = xold+timestep*vxold;
	ynew = yold+timestep*vyold;
	if ynew <=0 break;
	end
	vmag = sqrt(vxold*vxold+vyold*vyold);
	fx = -C*vxold*vmag;
	fy = -m*g-C*vyold*vmag;
	vxnew = vxold+(timestep*fx/m);
	vynew = vyold+(timestep*fy/m);
	xold=xnew;
	yold=ynew;
	vxold=vxnew;
	vyold=vynew;
	i=i+1;
	x(i)=xold;
	y(i)=yold;
end

plot(x,y)
axis([0 350 0 100])
hold on