clear
m =  7.26 ; %%initial mass
k = 0.47; %% drag coefficient of sphere
Vx = 15;
Vx2 = Vx
Vz = 10;
Vz2 = Vz
V= sqrt(Vx^2 + Vz^2);
area = 0.04
airden = 1.2; %% air density at sea level
fdrag = (airden*k*area)/2;
nodragx(1) =0;
nodragz(1) = 0;
x(1) = 0;
z(1)= 0;
t(1) =0
dt = 0.01
i=1
while (min(z)> -0.01)
       t = t + dt;
       i = i + 1;
       Ax = - ( fdrag / m ) * V * Vx;
       Az = -(9.8) - ( fdrag / m ) * V * Vz;
       Vx_new = Vx + Ax * dt;
       Vz_new = Vz + Az * dt;
       x(i) = x(i-1) + Vx_new * dt + 0.5 * Ax * dt^2;
       z(i) = z(i-1) + Vz_new * dt + 0.5 * Az * dt^2;
       Vx = Vx_new;
       Vz = Vz_new;
end
plot(x,z,'b'), hold on;
Vx2 = 10;
Vz2 = 15;
nodragx(1) =0;
nodragz(1) = 0;
t(1) =0
dt = 0.01
i=1
while (min(nodragz)> -0.01)
       t = t + dt;
       i = i + 1;
       Ax = 0;
       Az = -(9.8);
       Vx2_new = Vx2 + Ax * dt;
       Vz2_new = Vz2 + Az * dt;
       nodragx(i) = nodragx(i-1) + Vx2_new * dt + 0.5 * Ax * dt^2;
       nodragz(i) = nodragz(i-1) + Vz2_new * dt + 0.5 * Az * dt^2;
       Vx2 = Vx2_new;
       Vz2 = Vz2_new;
end
              %plots the Projectile Motion with Drag
plot(nodragx,nodragz,'r'), hold off;             %plots the Projectile Motion without Drag
xlabel('Horizontal Distance (m)');    %labels the x axis "Horizontal Distance (m)"
ylabel('Vertical Distance (m)');      %Labels the y axis "Vertical Distance (m)"
title('Projectile Motion Paths (Blue is with drag)');