def newton(v, alpha): if alpha <= 0 or alpha >= pi/2: print("Bad angle") return .0 def fx(t, v): return -(beta*v**2) def fy(t, v): return -copysign(1, v)*beta*v**2-g dvx, dvy = RK4(fx), RK4(fy) dt = 0.0001 vx, vy = v*cos(alpha), v*sin(alpha) t = dt x_old = y_old = 0 x = vx*dt y = vy*dt while (y > 0): t += dt vx_new = vx + dvx(t, vx, dt) vy_new = vy + dvy(t, vy, dt) x_old = x x += (vx_new + vx) * dt / 2 y_old = y y += (vy_new + vy) * dt / 2 vx = vx_new vy = vy_new if y_old == 0: return x_old return (x_old*(y_old - y) - y_old*(x_old-x))/(y_old - y)