density = 1.29;
mass = 0.3;
radius = 0.05;
c = 0.5;
vo = 1.61*125/3.6;
deltat = 0.01;
g = 9.8;
k = 0.5*c*density*Pi*radius^2/mass;
For[theta = 1, 
      theta < 90, 
angle = 
        theta/90*Pi/2; 
vx[1] = 
        vo*Cos[angle]; 
vy[1] = 
        vo*Sin[angle]; 
x[1] = 0; 
y[1] = 
        0; 
t[1] = 0; 
i = 
        1; 
While[
          y[i] ≥ 0, 
t[i + 1] = 
            t[i] + deltat; 
x[i + 1] = 
            x[i] + vx[i]*deltat; 
vx[i + 1] = 
            vx[i] - k*Sqrt[vx[i]^2 + vy[i]^2]*vx[i]*
                deltat; 
y[i + 1] = 
            y[i] + vy[i]*deltat; 
vy[i + 1] = 
            vy[i] - k*Sqrt[vx[i]^2 + vy[i]^2]*vy[i]*deltat - 
              g*deltat; 
i++;
]
If[
          y[i] == 0, 
Print[
            theta]; 
Print[
            x[i]]
]
If[
          y[i] < 0, 
r = -y[i - 1]/
              y[i]; 
Print[
            theta]; 
Print[(x[i - 1] + r*x[i])/(r + 
                  1)];
]
theta++;
]